<jupyter_start><jupyter_text>First, we specify several configuration hyperparameters and we store them in a dictionary.Not all of them are used at the same time. For example, if we decide to use an LSTM model, the parameters that specifies the TCN and RF model are not used.<jupyter_code>target_idx = [0]            # target variables to predict
B = 3                       # number of ensembles
alpha = 0.1                 # confidence level            
quantiles = [alpha/2,       # quantiles to predict
             0.5,
             1-(alpha/2)] 

# rf only
n_trees = 20                # number of trees in each rf model

# lstm and tcn only
regression = 'quantile'     # options: {'quantile', 'linear'}. If 'linear', just set one quantile
l2_lambda = 1e-4            # weight of l2 regularization in the lstm and tcn models
batch_size = 16             # size of batches using to train the lstm and tcn models

# lstm only
units = 128                 # number of units in each lstm layer
n_layers = 3                # number of lstm layers in the model

# tcn only
dilations = [1,2,4,8]       # dilation rate of the Conv1D layers
n_filters = 128             # filters in each Conv1D layer 
kernel_size = 7             # kernel size in each ConvID layer

# Store the configuration in a dictionary
P = {'B':B, 'alpha':alpha, 'quantiles':quantiles,
     'n_trees':n_trees,  
     'regression':regression,'l2':l2_lambda, 'batch_size':batch_size,
     'units':units,'n_layers':n_layers,
     'dilations':dilations, 'n_filters':n_filters, 'kernel_size':kernel_size}<jupyter_output><empty_output><jupyter_text>Data loadingFor this example, we will use 3 years of data relative to solar power production.We will use the first year for training, the second for validation, and the last year as test set.*Note:* To use your own data, you must write a data loader which returns a single DataFrame or, like in this case, 3 DataFrames (one for training, one for validation, and one for test). You can find two examples of data loaders in [data_loaders.py](https://github.com/FilippoMB/Ensemble-Conformalized-Quantile-Regression/blob/main/data_loaders.py).<jupyter_code>train_df, val_df, test_df = data_loaders.get_solar_data()

train_df.head()<jupyter_output><empty_output><jupyter_text>Data preprocessingThe ``data_windowing()`` function transforms each DataFrame into 3-dimensional arrays of shape \[*number of samples*, *time steps*, *number of variables* \].The input data, X, might have a different number of time steps and a different number of variables than the output data, Y. In this case, we want to predict the energy production for the next day given the measurements of the past week. Therefore, the second dimension of X is ``time_steps_in=168`` (hours in the past week) and the second dimension of Y is ``time_steps_out=24`` (hours of the next day). The input variables are the historical energy production plus 5 exogenous variables, so the last dimension of X is ``n_vars=6``. Since we want to predict the future energy production, we specify the target variable to predict: ``label_columns=['MWH']``. Note that in Y ``n_vars=1``.``data_windowing()`` also rescales each variable in \[0,1\] and return the scaler, which is used to invert the transformation.In addition, it also splits training data in *B* disjoint sets, used to train the ensemble model. In this case, ``B=3``.![data_shape.drawio.png](attachment:data_shape.drawio.png)<jupyter_code>train_data, val_x, val_y, test_x, test_y, Scaler = data_preprocessing.data_windowing(df=train_df, 
                                                                                     val_data=val_df,
                                                                                     test_data=test_df,
                                                                                     B=3, 
                                                                                     time_steps_in=168, 
                                                                                     time_steps_out=24, 
                                                                                     label_columns=['MWH'])

print("-- Training data --")
for i in range(len(train_data)):
    print(f"Set {i} - x: {train_data[i][0].shape}, y: {train_data[i][1].shape}")
print("-- Validation data --")
print(f"x: {val_x.shape}, y: {val_y.shape}")
print("-- Test data --")
print(f"x: {test_x.shape}, y: {test_y.shape}")

# Update configuration dict
P['time_steps_in'] = test_x.shape[1]
P['n_vars'] = test_x.shape[2] 
P['time_steps_out'] = test_y.shape[1]<jupyter_output>-- Training data --
Set 0 - x: (119, 168, 6), y: (119, 24)
Set 1 - x: (119, 168, 6), y: (119, 24)
Set 2 - x: (119, 168, 6), y: (119, 24)
-- Validation data --
x: (357, 168, 6), y: (357, 24)
-- Test data --
x: (357, 168, 6), y: (357, 24)<jupyter_text>Training the quantile regression modelsBefore looking into the conformalization of the PI, let's see how we can train different models that perform quantile regression.In the paper we considered three models:- a random forest (rf)- a recurrent neural network with LSTM cells- a feedforward neural network with 1-dimensional convolutional cells (TCN).In principle, any other model performing quantile regression can be used. Each model must implement a ``fit()`` function with is used to train the model parameters and a ``transform()`` function used to predict new data.The ``fit()`` function uses ``val_x`` and ``val_y`` to perform early stopping.Let's start with the **TCN** model.<jupyter_code>P['model_type'] = 'tcn' 

# Train 
model = regression_model(P)
hist = model.fit(train_data[0][0], train_data[0][1], val_x, val_y)
utils.plot_history(hist)

# Test
PI = model.transform(test_x)
utils.plot_PIs(test_y, PI[:,:,1],
                PI[:,:,0], PI[:,:,2],
                x_lims=[0,168], scaler=Scaler, title='TCN model')<jupyter_output><empty_output><jupyter_text>The function ``plot_hist()`` plots how the loss, coverage, and PI length evolve during training on the train and validation set.Note that here we trained the model only on the first subset of the training set.Next we train the **LSTM** model. To do that, we just change ``model_type`` in the hyperparameters dictionary.<jupyter_code>P['model_type'] = 'lstm'

# Train
model = regression_model(P)
hist = model.fit(train_data[0][0], train_data[0][1], val_x, val_y)
utils.plot_history(hist)

# Test
PI = model.transform(test_x)
utils.plot_PIs(test_y, PI[:,:,1],
                PI[:,:,0], PI[:,:,2],
                x_lims=[0,168], scaler=Scaler, title='LSTM model')<jupyter_output><empty_output><jupyter_text>Finally, we train the **RF** model. As before, we change  ``model_type`` in the hyperparameters dictionary. Contrairly to the previous two neural network model, the ``fit()`` function does not use ``val_x`` and ``val_y`` since there is no early stopping.<jupyter_code># Train
P['model_type'] = 'rf'
model = regression_model(P)
model.fit(train_data[0][0], train_data[0][1])

# Test
PI = model.transform(test_x)
utils.plot_PIs(test_y, PI[:,:,1],
                PI[:,:,0], PI[:,:,2],
                x_lims=[0,168], scaler=Scaler, title='RF model')<jupyter_output><empty_output><jupyter_text>EnCQRFinally, we compute the intervals with the EnCQR method.This is done by calling the function ``conformalized_PI()``, which returns two intervals:- the PI computed by the ensemble of QR models- the conformalized PIIn this example, we consider an ensemble of TCN models and show that after conformalization the coverage of the PI gets much closer to the desired confidence level.<jupyter_code>P['model_type'] = 'tcn'

# compute the conformalized PI with EnCQR
PI, conf_PI = EnCQR(train_data, val_x, val_y, test_x, test_y, P)

# Plot original and conformalized PI
utils.plot_PIs(test_y, PI[:,:,1],
               PI[:,:,0], PI[:,:,2],
               conf_PI[:,:,0], conf_PI[:,:,2],
               x_lims=[0,168], scaler=Scaler)

# Compute PI coverage and length before and after conformalization
print("Before conformalization:")
utils.compute_coverage_len(test_y.flatten(), PI[:,:,0].flatten(), PI[:,:,2].flatten(), verbose=True)
print("After conformalization:")
utils.compute_coverage_len(test_y.flatten(), conf_PI[:,:,0].flatten(), conf_PI[:,:,2].flatten(), verbose=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11517696460707859, 'fraction_numerical': 0.020155968806238753, 'mean_word_length': 3.231472081218274, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '507443', 'license': 'MIT', 'n_tokens_mistral': 2616, 'n_tokens_neox': 2466, 'n_words': 873, 'path': 'example.ipynb', 'repo_name': 'FilippoMB/Ensemble-Conformalized-Quantile-Regression-forProbabilistic-Time-Series-Forecasting', 'size': '8335'}