Probabilistic Forecasting: Global Models¶
Skforecast allows applying all its implemented probabilistic forecasting methods (bootstrapping, conformal prediction and quantile regression) to global models. A global model, such as ForecasterRecursiveMultiSeries, is a single model trained with the observations of all the available time series, which can then forecast any of them, and even series not seen during training. This guide shows how to estimate prediction intervals for multiple series with conformal prediction and bootstrapped residuals.
For detailed information about the available probabilistic forecasting methods, see the following user guides:
💡 Tip
When predicting multiple series and scaling up to hundreds or thousands of series, computational time can become a bottleneck. In such cases, the conformal framework is a good option, since it is usually faster than the other methods: it does not need to simulate n_boot prediction paths (bootstrapping), nor to train one model per quantile (quantile regression).
Libraries and data¶
# Data processing
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset
# Plots
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme, plot_prediction_intervals
from statsmodels.graphics.tsaplots import plot_pacf
# Modelling and Forecasting
# ==============================================================================
from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from lightgbm import LGBMRegressor
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import TimeSeriesFold, backtesting_forecaster_multiseries
from skforecast.metrics import calculate_coverage
from skforecast.stats import calculate_lag_autocorrelation
# Configuration
# ==============================================================================
import warnings
from pprint import pprint
warnings.filterwarnings('once')
# Data
# ==============================================================================
data = fetch_dataset(name="bdg2_hourly_sample")
data = data.loc['2016-01-01 00:00:00' : '2016-09-30 00:00:00', :]
series = ['building_1', 'building_2']
data.head(2)
╭────────────────────────────── bdg2_hourly_sample ───────────────────────────────╮ │ Description: │ │ Hourly energy consumption data of two buildings sampled from the The Building │ │ Data Genome Project 2. https://github.com/buds-lab/building-data-genome- │ │ project-2 │ │ │ │ Source: │ │ Miller, C., Kathirgamanathan, A., Picchetti, B. et al. The Building Data Genome │ │ Project 2, energy meter data from the ASHRAE Great Energy Predictor III │ │ competition. Sci Data 7, 368 (2020). https://doi.org/10.1038/s41597-020-00712-x │ │ │ │ URL: │ │ https://raw.githubusercontent.com/skforecast/skforecast- │ │ datasets/main/data/bdg2_hourly_sample.csv │ │ │ │ Shape: 6553 rows x 2 columns │ ╰─────────────────────────────────────────────────────────────────────────────────╯
| building_1 | building_2 | |
|---|---|---|
| timestamp | ||
| 2016-01-01 00:00:00 | 186.532 | 219.27 |
| 2016-01-01 01:00:00 | 186.532 | 219.27 |
The day of the week and the hour of the day are used as exogenous variables to capture the weekly and daily patterns of energy consumption. Both are one-hot encoded, so that each category becomes a binary column. Since calendar features are deterministic and known in advance for any date, encoding them with the whole dataset does not leak information from the future.
# Calendar features
# ==============================================================================
data['day_of_week'] = data.index.dayofweek
data['hour'] = data.index.hour
transformer = ColumnTransformer(
transformers=[
('one_hot_encoder', OneHotEncoder(sparse_output=False), ['day_of_week', 'hour'])
],
remainder='passthrough',
verbose_feature_names_out=False
).set_output(transform='pandas')
data = transformer.fit_transform(data)
data.head()
| day_of_week_0 | day_of_week_1 | day_of_week_2 | day_of_week_3 | day_of_week_4 | day_of_week_5 | day_of_week_6 | hour_0 | hour_1 | hour_2 | ... | hour_16 | hour_17 | hour_18 | hour_19 | hour_20 | hour_21 | hour_22 | hour_23 | building_1 | building_2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| timestamp | |||||||||||||||||||||
| 2016-01-01 00:00:00 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 186.532 | 219.270 |
| 2016-01-01 01:00:00 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 186.532 | 219.270 |
| 2016-01-01 02:00:00 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 183.479 | 218.348 |
| 2016-01-01 03:00:00 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 184.303 | 220.003 |
| 2016-01-01 04:00:00 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 183.657 | 217.731 |
5 rows × 33 columns
# Plot time series
# ==============================================================================
set_dark_theme()
plt.rcParams['lines.linewidth'] = 0.5
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(2, 1, figsize=(7, 4), sharex=True)
for i, col in enumerate(series):
axs[i].plot(data[col], label=col, color=colors[i])
axs[i].legend(loc='lower left', fontsize=8)
axs[i].tick_params(axis='both', labelsize=8)
plt.tight_layout()
Lag selection¶
The partial autocorrelation function (PACF) measures the correlation between a series and its lag after removing the effect of the intermediate lags ( to ). Lags with a high absolute partial autocorrelation carry information that the previous lags do not, so they are good candidates for predictors. The partial autocorrelation of each building is calculated up to one week ( lags) with calculate_lag_autocorrelation, and the 10 lags with the highest absolute value are selected for each series. Since a global model uses the same lags for all the series, the union of both sets is used as predictors.
# Partial autocorrelation values and plots
# ==============================================================================
n_lags = 7 * 24
pacf_dfs = []
fig, axs = plt.subplots(2, 1, figsize=(6, 3))
for i, col in enumerate(series):
pacf_values = calculate_lag_autocorrelation(data[col], n_lags=n_lags)
pacf_values['variable'] = col
pacf_dfs.append(pacf_values)
plot_pacf(data[col], lags=n_lags, ax=axs[i])
axs[i].set_title(col, fontsize=10)
axs[i].set_ylim(-0.5, 1.1)
plt.tight_layout()
# Top n lags with highest absolute partial autocorrelation per variable
# ==============================================================================
n = 10
top_lags = set()
for pacf_values in pacf_dfs:
variable = pacf_values['variable'].iloc[0]
lags = pacf_values.nlargest(n, 'partial_autocorrelation_abs')['lag'].sort_values().tolist()
top_lags.update(lags)
print(f"{variable}: {lags}")
top_lags = sorted(top_lags)
print(f"\nAll lags: {top_lags}")
building_1: [1, 2, 3, 15, 16, 19, 20, 24, 25, 145] building_2: [1, 2, 3, 20, 21, 22, 25, 26, 141, 145] All lags: [1, 2, 3, 15, 16, 19, 20, 21, 22, 24, 25, 26, 141, 145]
The target series exhibit similar dynamics, and the selected lags can be grouped into three blocks:
Lags 1 to 3: the consumption of the last few hours (short-term inertia).
Lags 15 to 26: values from the previous 15 to 26 hours, mostly around 24 hours, which capture the daily cycle.
Lags 141 and 145: values from around six days before, related to the weekly pattern.
The union of the lags selected for both series contains 14 lags, which are used as predictors by the global model.
For simplicity, the partial autocorrelation is calculated with the whole dataset. In a strict evaluation, lag selection should use only the training data, since otherwise information from the test period influences the design of the model.
Train, validation and test partitions¶
The data is split into three consecutive partitions, each with a different role:
Train: used to fit the forecaster.
Validation (calibration): the forecaster, trained only with the training data, predicts this period through backtesting. The errors made (out-of-sample residuals) are used later to build the prediction intervals.
Test: used to evaluate the coverage and width of the prediction intervals on data not involved in any of the previous steps.
Label-based slicing with .loc includes both ends. The end dates are set to 23:59, a time that does not exist in the hourly index, so that consecutive partitions do not share any observation.
# Split train-validation-test
# ==============================================================================
end_train = '2016-07-01 23:59:00'
end_validation = '2016-09-01 23:59:00'
data_train = data.loc[: end_train, :]
data_val = data.loc[end_train:end_validation, :]
data_test = data.loc[end_validation:, :]
print(f"Dates train : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Dates validation : {data_val.index.min()} --- {data_val.index.max()} (n={len(data_val)})")
print(f"Dates test : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
Dates train : 2016-01-01 00:00:00 --- 2016-07-01 23:00:00 (n=4392) Dates validation : 2016-07-02 00:00:00 --- 2016-09-01 23:00:00 (n=1488) Dates test : 2016-09-02 00:00:00 --- 2016-09-30 00:00:00 (n=673)
# Plot partitions
# ==============================================================================
fig, axs = plt.subplots(2, 1, figsize=(7, 5), sharex=True)
for i, col in enumerate(series):
axs[i].plot(data[col], label=col, color=colors[i])
axs[i].legend(loc='lower left', fontsize=8)
axs[i].tick_params(axis='both', labelsize=8)
axs[i].axvline(pd.to_datetime(end_train), color='white', linestyle='--', linewidth=1) # End train
axs[i].axvline(pd.to_datetime(end_validation), color='white', linestyle='--', linewidth=1) # End validation
plt.tight_layout()
Forecaster¶
A ForecasterRecursiveMultiSeries is created with a LightGBM regressor. It is a global model: a single estimator is trained with the observations of both buildings, using the 14 lags selected above as predictors. With encoding='ordinal', the identifier of each series is added as an extra feature (encoded as an integer), so that the model can learn patterns specific to each building. The one-hot encoded calendar features are passed as exogenous variables.
# Create forecaster
# ==============================================================================
exog_features = data.columns[data.columns.str.contains('day_of_|hour_')].tolist()
params = {
"n_estimators": 300,
"learning_rate": 0.05,
"max_depth": 5,
}
forecaster = ForecasterRecursiveMultiSeries(
estimator = LGBMRegressor(random_state=15926, verbose=-1, **params),
lags = top_lags,
encoding = 'ordinal',
)
Out-of-sample residuals¶
To estimate prediction intervals that reflect the errors expected on new data, the residuals must be out-of-sample, that is, generated by a model that did not see those observations during training. They are obtained by running a backtesting over the validation partition with backtesting_forecaster_multiseries. The folds are defined with TimeSeriesFold: the forecaster is trained with the training data only (initial_train_size = len(data_train)), and the validation period is predicted in folds of 24 steps (one day) without refitting. As a result, the residuals contain the errors of steps 1 to 24 of the forecast horizon, the same horizon used later on the test set.
# Backtesting on validation data to obtain out-of-sample residuals
# ==============================================================================
cv = TimeSeriesFold(
initial_train_size = len(data_train),
steps = 24,
)
metric_val, predictions_val = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data.loc[:end_validation, series],
exog = data.loc[:end_validation, exog_features],
cv = cv,
metric = "mean_absolute_error"
)
metric_val
╭────────────────────────────────── InputTypeWarning ──────────────────────────────────╮ │ Passing a DataFrame (either wide or long format) as `series` requires additional │ │ internal transformations, which can increase computational time. It is recommended │ │ to use a dictionary of pandas Series instead. For more details, see: │ │ https://skforecast.org/latest/user_guides/independent-multi-time-series-forecasting. │ │ html#input-data │ │ │ │ Category : skforecast.exceptions.InputTypeWarning │ │ Location : │ │ /Users/javier.escobar/code/github/skforecast/skforecast/utils/utils.py:4204 │ │ Suppress : warnings.simplefilter('ignore', category=InputTypeWarning) │ ╰──────────────────────────────────────────────────────────────────────────────────────╯
| levels | mean_absolute_error | |
|---|---|---|
| 0 | building_1 | 6.564278 |
| 1 | building_2 | 9.104094 |
| 2 | average | 7.834186 |
| 3 | weighted_average | 7.834186 |
| 4 | pooling | 7.834186 |
# Plot predictions on validation data
# ==============================================================================
fig, axs = plt.subplots(2, 1, figsize=(7, 5.5), sharex=True, sharey=False)
for i, level in enumerate(predictions_val['level'].unique()):
predictions_val_level = predictions_val.loc[predictions_val['level'] == level, 'pred']
data.loc[end_train:end_validation, level].plot(ax=axs[i], label='Real value')
predictions_val_level.plot(ax=axs[i], label='Prediction')
axs[i].set_xlabel("")
axs[i].set_title(level)
axs[i].legend(loc='lower right', fontsize=7)
plt.tight_layout()
# Out-of-sample residuals over time
# ==============================================================================
fig, axs = plt.subplots(2, 1, figsize=(7, 4), sharex=True, sharey=False)
for i, level in enumerate(predictions_val["level"].unique()):
residuals = (
data.loc[predictions_val.index, level]
- predictions_val.loc[predictions_val["level"] == level, "pred"]
)
residuals.plot(ax=axs[i], label=f"Residuals {level}")
axs[i].legend(loc="upper right", fontsize=7)
fig.tight_layout()
Next, the forecaster is trained with the same training data used in the backtesting, and the validation residuals are stored in it with the set_out_sample_residuals() method. Both y_true and y_pred are dictionaries whose keys are the names of the series, so that each series keeps its own residuals. In addition, the forecaster stores a random sample of the residuals of all the series under the key _unknown_level, which is used to estimate intervals for series not seen during training.
# Store out-of-sample residuals in the forecaster
# ==============================================================================
forecaster.fit(
series = data_train[series],
exog = data_train[exog_features],
suppress_warnings = True
)
forecaster.set_out_sample_residuals(
y_true = {k: data_val[k] for k in series},
y_pred = {k: v for k, v in predictions_val.groupby('level')['pred']}
)
✏️ Note
When the out-of-sample residuals are stored in the forecaster, they are binned according to the predicted value to which they correspond (10 bins by default, as shown below). Later, when the intervals are estimated with use_binned_residuals=True, only the residuals of the bin where the prediction falls are used: in bootstrapping, the residuals are sampled from that bin; in conformal prediction, the correction factor is calculated from that bin. In this way, the width of the intervals adapts to the magnitude of the prediction, which is useful when the size of the errors depends on the predicted value. For more information about how residuals are used in interval estimation, visit Probabilistic forecasting: Bootstrapped residuals and Probabilistic forecasting: Conformal prediction.
# Intervals of the residual bins (conditioned on predicted values) for each level
# ==============================================================================
for k, v in forecaster.binner_intervals_.items():
print(k)
pprint(v)
print("")
building_1
{0: (165.20499932654647, 188.29862772750283),
1: (188.29862772750283, 192.71291989855519),
2: (192.71291989855519, 196.50159701655218),
3: (196.50159701655218, 200.36282013339397),
4: (200.36282013339397, 205.78088677076306),
5: (205.78088677076306, 213.05048123203005),
6: (213.05048123203005, 223.81891494270874),
7: (223.81891494270874, 247.23223476967067),
8: (247.23223476967067, 257.59956347695794),
9: (257.59956347695794, 292.5578128598947)}
building_2
{0: (185.31199428441397, 221.82899062428837),
1: (221.82899062428837, 227.77775215345284),
2: (227.77775215345284, 231.22953122462755),
3: (231.22953122462755, 234.1632458534549),
4: (234.1632458534549, 237.5259183049393),
5: (237.5259183049393, 241.3815688800525),
6: (241.3815688800525, 260.5481777111192),
7: (260.5481777111192, 288.13779883487916),
8: (288.13779883487916, 299.9804948765519),
9: (299.9804948765519, 316.86744258452006)}
_unknown_level
{0: (165.20499932654647, 192.70687277925953),
1: (192.70687277925953, 200.33743695492512),
2: (200.33743695492512, 211.78881040116786),
3: (211.78881040116786, 222.55028533988022),
4: (222.55028533988022, 230.37899103921075),
5: (230.37899103921075, 235.9675135656299),
6: (235.9675135656299, 243.06523999048886),
7: (243.06523999048886, 258.36155747914506),
8: (258.36155747914506, 288.367317637904),
9: (288.367317637904, 316.86744258452006)}
# Distribution of the residuals by bin and level
# ==============================================================================
flierprops = dict(marker='o', markerfacecolor='gray', markersize=5, linestyle='none')
levels = list(forecaster.out_sample_residuals_by_bin_)
fig, axs = plt.subplots(len(levels), 1, figsize=(7, 7), sharex=True, sharey=True)
for i, level in enumerate(levels):
residuals_by_bin = forecaster.out_sample_residuals_by_bin_[level]
residuals_by_bin_df = pd.DataFrame(
{k: pd.Series(v) for k, v in residuals_by_bin.items()}
)
residuals_by_bin_df.boxplot(ax=axs[i], flierprops=flierprops)
axs[i].set_title(level, fontsize=10)
axs[i].set_ylabel("Residual", fontsize=8)
axs[-1].set_xlabel("Bin", fontsize=8)
fig.suptitle("Distribution of residuals by bin", fontsize=12)
fig.tight_layout()
✏️ Note
Two arguments control the use of residuals in predict_interval() and backtesting_forecaster_multiseries():
use_in_sample_residuals: IfTrue, the in-sample residuals are used to compute the prediction intervals. Since these residuals are obtained from the training set, they are always available, but usually lead to overoptimistic intervals. IfFalse, the out-of-sample residuals (calibration) are used to calculate the prediction intervals. These residuals are obtained from the validation/calibration set and are only available if theset_out_sample_residuals()method has been called. It is recommended to use out-of-sample residuals to achieve the desired coverage.use_binned_residuals: IfFalse, all the residuals are used regardless of the predicted value: a single correction factor is applied to all the predictions (conformal), or the residuals are sampled from the whole set (bootstrapping). IfTrue, only the residuals of the bin where the prediction falls are used, so the correction factor (conformal) or the sampled residuals (bootstrapping) depend on the predicted value. This produces intervals whose width adapts to the magnitude of the prediction.
Intervals using conformal prediction¶
The forecaster is now evaluated on the test set with a backtesting of 24-step folds. The argument interval = [0.1, 0.9] requests the 0.1 and 0.9 quantiles, that is, an interval with a nominal coverage of 80%. The intervals are calculated with interval_method = "conformal", using the out-of-sample residuals conditioned on the predicted value.
Two metrics are used to evaluate the intervals on the test set:
Empirical coverage: proportion of observed values that fall within the interval, calculated with
calculate_coverage. Well-calibrated intervals have a coverage close to the nominal one (80%).Area: sum of the widths of the interval over all the predicted steps. For a similar coverage, a smaller area means narrower, more informative intervals.
The intervals are plotted with plot_prediction_intervals.
✏️ Note
Since initial_train_size covers the training and validation partitions, the backtesting retrains the forecaster with both of them before predicting the test set, so the final model uses all the data available before the test period. The out-of-sample residuals and their bins, stored with set_out_sample_residuals(), are kept after this refit, so the intervals are still calibrated with the validation residuals.
# Backtesting with conformal prediction intervals
# ==============================================================================
cv = TimeSeriesFold(
initial_train_size = len(data.loc[:end_validation, :]),
steps = 24,
)
metric, predictions_conformal = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data[series],
exog = data[exog_features],
cv = cv,
metric = "mean_absolute_error",
interval = [0.1, 0.9],
interval_method = "conformal",
use_in_sample_residuals = False, # Use out-of-sample residuals
use_binned_residuals = True, # Adaptive intervals
suppress_warnings = True
)
predictions_conformal
| level | fold | pred | lower_bound | upper_bound | |
|---|---|---|---|---|---|
| 2016-09-02 00:00:00 | building_1 | 0 | 198.821665 | 189.903492 | 207.739838 |
| 2016-09-02 00:00:00 | building_2 | 0 | 204.475048 | 194.299621 | 214.650475 |
| 2016-09-02 01:00:00 | building_1 | 0 | 195.874267 | 188.227810 | 203.520723 |
| 2016-09-02 01:00:00 | building_2 | 0 | 202.250982 | 192.075555 | 212.426409 |
| 2016-09-02 02:00:00 | building_1 | 0 | 192.720768 | 185.074312 | 200.367225 |
| ... | ... | ... | ... | ... | ... |
| 2016-09-29 22:00:00 | building_2 | 27 | 212.441112 | 202.265685 | 222.616539 |
| 2016-09-29 23:00:00 | building_1 | 27 | 205.382013 | 194.846273 | 215.917753 |
| 2016-09-29 23:00:00 | building_2 | 27 | 206.619372 | 196.443945 | 216.794799 |
| 2016-09-30 00:00:00 | building_1 | 28 | 200.266870 | 191.348697 | 209.185043 |
| 2016-09-30 00:00:00 | building_2 | 28 | 203.849044 | 193.673617 | 214.024471 |
1346 rows × 5 columns
# Function to plot intervals and calculate coverage for each level
# ==============================================================================
def plot_intervals_and_coverage(predictions, y_true):
"""
Plot the prediction intervals of each level and print their empirical
coverage and area.
"""
plt.rcParams["lines.linewidth"] = 1
for level in predictions["level"].unique():
print(f"\nLevel: {level}")
predictions_level = (
predictions[predictions["level"] == level].drop(columns="level")
)
# Plot intervals
fig, ax = plt.subplots(figsize=(7, 3))
plot_prediction_intervals(
predictions = predictions_level,
y_true = y_true[[level]],
target_variable = level,
initial_x_zoom = None,
title = "Prediction intervals",
xaxis_title = "",
yaxis_title = level,
ax = ax
)
ax.legend(loc="upper left", fontsize=7)
fill_between_obj = ax.collections[0]
fill_between_obj.set_facecolor("white")
fill_between_obj.set_alpha(0.3)
# Empirical coverage of the interval
coverage = calculate_coverage(
y_true = y_true[level],
lower_bound = predictions_level["lower_bound"],
upper_bound = predictions_level["upper_bound"]
)
print(f"Empirical coverage of the interval: {round(100 * coverage, 2)} %")
# Area of the interval
area = (predictions_level["upper_bound"] - predictions_level["lower_bound"]).sum()
print(f"Area of the interval: {round(area, 2)}")
plot_intervals_and_coverage(predictions_conformal, y_true=data_test)
Level: building_1 Empirical coverage of the interval: 83.95 % Area of the interval: 13485.99 Level: building_2 Empirical coverage of the interval: 83.06 % Area of the interval: 17790.12
With conformal prediction, the empirical coverage on the test set is 83.95% for building_1 and 83.06% for building_2, slightly above the nominal coverage of 80%. The intervals are slightly conservative, but reasonably well calibrated. Their areas are 13,485.99 and 17,790.12, respectively. The larger area of building_2 is consistent with its larger error in the validation backtesting (mean absolute error of 9.10 versus 6.56 for building_1).
Intervals using bootstrapped residuals¶
The same process is repeated, but this time the bootstrapping method is used instead of the conformal method. With n_boot = 150, 150 alternative prediction paths are simulated by adding to the predictions residuals sampled from the bin of each prediction. The bounds of the interval are the 0.1 and 0.9 quantiles of these paths at each step.
# Backtesting with prediction intervals using bootstrapping
# ==============================================================================
# The cv object of the conformal backtesting is reused
metric, predictions_bootstrapping = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data[series],
exog = data[exog_features],
cv = cv,
metric = "mean_absolute_error",
interval = [0.1, 0.9],
interval_method = "bootstrapping",
n_boot = 150,
use_in_sample_residuals = False, # Use out-of-sample residuals
use_binned_residuals = True, # Residuals conditioned on predicted values
suppress_warnings = True
)
predictions_bootstrapping
| level | fold | pred | lower_bound | upper_bound | |
|---|---|---|---|---|---|
| 2016-09-02 00:00:00 | building_1 | 0 | 198.821665 | 192.420494 | 209.905264 |
| 2016-09-02 00:00:00 | building_2 | 0 | 204.475048 | 194.847599 | 215.465116 |
| 2016-09-02 01:00:00 | building_1 | 0 | 195.874267 | 186.090498 | 213.503387 |
| 2016-09-02 01:00:00 | building_2 | 0 | 202.250982 | 188.650391 | 218.721453 |
| 2016-09-02 02:00:00 | building_1 | 0 | 192.720768 | 183.154797 | 218.639839 |
| ... | ... | ... | ... | ... | ... |
| 2016-09-29 22:00:00 | building_2 | 27 | 212.441112 | 192.717301 | 234.768066 |
| 2016-09-29 23:00:00 | building_1 | 27 | 205.382013 | 193.015801 | 252.393190 |
| 2016-09-29 23:00:00 | building_2 | 27 | 206.619372 | 183.959878 | 230.196653 |
| 2016-09-30 00:00:00 | building_1 | 28 | 200.266870 | 193.766870 | 209.104257 |
| 2016-09-30 00:00:00 | building_2 | 28 | 203.849044 | 193.687843 | 215.279858 |
1346 rows × 5 columns
# Plot intervals and calculate coverage for each level
# ==============================================================================
plot_intervals_and_coverage(predictions_bootstrapping, y_true=data_test)
Level: building_1 Empirical coverage of the interval: 96.29 % Area of the interval: 31297.34 Level: building_2 Empirical coverage of the interval: 97.92 % Area of the interval: 37029.52
With bootstrapping, the empirical coverage is 96.29% for building_1 and 97.92% for building_2, well above the nominal 80%. The intervals are overly conservative: their areas (31,297.34 and 37,029.52) are about 2.3 and 2.1 times larger than those obtained with conformal prediction.
A likely reason is the way the residuals are used. The out-of-sample residuals come from a backtesting with folds of 24 steps, so each residual already contains the error accumulated over up to 24 predicted steps. In recursive bootstrapping, a residual is added at every step and, since each simulated value is used as a lag to predict the following steps, the noise propagates along the horizon. Part of the accumulated error is therefore counted twice. Conformal prediction, instead, adds the correction factor once per step, without propagating it. In addition, as discussed in the bootstrapped residuals user guide, the residuals are obtained with a model trained only with the training partition, while the final model is also trained with the validation data and is expected to make smaller errors.
In this example, conformal prediction provides intervals that are much closer to the nominal coverage and considerably narrower.
Forecasting intervals for unknown series¶
The ForecasterRecursiveMultiSeries class allows forecasting series not seen during training, with the only requirement that at least window_size observations of the new series are available (145 in this example, the largest lag). These observations are passed through the last_window argument of predict_interval().
In this example, the new series is simulated by adding a small amount of Gaussian noise to the last observations of building_1 in the validation partition. Two warnings are raised:
The identifier of the new series was not seen during training, so it is encoded as
NaN. This works because LightGBM handles missing values natively; with an estimator that does not acceptNaN, the prediction would fail.There are no residuals for the new series, so the intervals are calculated with the residuals stored under the key
_unknown_level, a random sample of the residuals of all the known series.
# Predictions for an unknown series
# ==============================================================================
# Simulate last_window of the new series
rng = np.random.default_rng(123)
len_last_window = forecaster.window_size
last_window = pd.DataFrame({
'new_series': (
data_val['building_1'].iloc[-len_last_window:]
+ rng.normal(0, 0.1, len_last_window)
)
})
forecaster.predict_interval(
steps = 24,
last_window = last_window,
exog = data_test[exog_features],
method = 'conformal',
interval = [0.1, 0.9],
use_in_sample_residuals = False
)
╭──────────────────────────────── UnknownLevelWarning ─────────────────────────────────╮ │ `levels` {'new_series'} were not included in training. Unknown levels are encoded as │ │ NaN, which may cause the prediction to fail if the estimator does not accept NaN │ │ values. │ │ │ │ Category : skforecast.exceptions.UnknownLevelWarning │ │ Location : │ │ /Users/javier.escobar/code/github/skforecast/skforecast/utils/utils.py:1491 │ │ Suppress : warnings.simplefilter('ignore', category=UnknownLevelWarning) │ ╰──────────────────────────────────────────────────────────────────────────────────────╯
╭──────────────────────────────── UnknownLevelWarning ─────────────────────────────────╮ │ `levels` {'new_series'} are not present in │ │ `forecaster.out_sample_residuals_by_bin_`. A random sample of the residuals from │ │ other levels will be used. This can lead to inaccurate intervals for the unknown │ │ levels. Otherwise, Use the `set_out_sample_residuals()` method before predicting to │ │ set the residuals for these levels. │ │ │ │ Category : skforecast.exceptions.UnknownLevelWarning │ │ Location : │ │ /Users/javier.escobar/code/github/skforecast/skforecast/utils/utils.py:1925 │ │ Suppress : warnings.simplefilter('ignore', category=UnknownLevelWarning) │ ╰──────────────────────────────────────────────────────────────────────────────────────╯
| level | pred | lower_bound | upper_bound | |
|---|---|---|---|---|
| 2016-09-02 00:00:00 | new_series | 199.294313 | 190.929439 | 207.659187 |
| 2016-09-02 01:00:00 | new_series | 196.150552 | 187.785678 | 204.515426 |
| 2016-09-02 02:00:00 | new_series | 192.871179 | 184.506305 | 201.236053 |
| 2016-09-02 03:00:00 | new_series | 191.774456 | 185.001310 | 198.547603 |
| 2016-09-02 04:00:00 | new_series | 191.367051 | 184.593904 | 198.140198 |
| 2016-09-02 05:00:00 | new_series | 191.546486 | 184.773339 | 198.319632 |
| 2016-09-02 06:00:00 | new_series | 192.314902 | 185.541755 | 199.088049 |
| 2016-09-02 07:00:00 | new_series | 198.286753 | 189.921879 | 206.651627 |
| 2016-09-02 08:00:00 | new_series | 215.400580 | 203.870512 | 226.930649 |
| 2016-09-02 09:00:00 | new_series | 236.501373 | 224.016532 | 248.986214 |
| 2016-09-02 10:00:00 | new_series | 249.815882 | 234.829805 | 264.801959 |
| 2016-09-02 11:00:00 | new_series | 252.806961 | 237.820885 | 267.793038 |
| 2016-09-02 12:00:00 | new_series | 255.773568 | 240.787491 | 270.759644 |
| 2016-09-02 13:00:00 | new_series | 255.838280 | 240.852204 | 270.824357 |
| 2016-09-02 14:00:00 | new_series | 255.514821 | 240.528744 | 270.500898 |
| 2016-09-02 15:00:00 | new_series | 253.333388 | 238.347311 | 268.319464 |
| 2016-09-02 16:00:00 | new_series | 249.001096 | 234.015020 | 263.987173 |
| 2016-09-02 17:00:00 | new_series | 238.398423 | 225.913582 | 250.883264 |
| 2016-09-02 18:00:00 | new_series | 227.150036 | 215.931363 | 238.368710 |
| 2016-09-02 19:00:00 | new_series | 215.810711 | 204.280643 | 227.340780 |
| 2016-09-02 20:00:00 | new_series | 210.391208 | 199.368534 | 221.413881 |
| 2016-09-02 21:00:00 | new_series | 209.709003 | 198.686329 | 220.731676 |
| 2016-09-02 22:00:00 | new_series | 208.672495 | 197.649822 | 219.695169 |
| 2016-09-02 23:00:00 | new_series | 204.863147 | 193.840474 | 215.885821 |
Probabilistic forecasting in production¶
Once the forecaster is deployed, the out-of-sample residuals should be refreshed periodically with the most recent errors of the model. A simple workflow is to keep the predictions made by the deployed forecaster (with the same horizon used in production) and, once the real values are observed, store the new residuals with set_out_sample_residuals(). With append=False (default), the stored residuals are replaced; with append=True, the new residuals are added to the existing ones.
# Illustrative example: y_pred_recent contains the predictions made by the
# deployed forecaster over the last weeks and y_true_recent the values observed
# afterwards, both as dictionaries {series_name: pandas Series}
forecaster.set_out_sample_residuals(
y_true = y_true_recent,
y_pred = y_pred_recent,
append = False # True to add the new residuals to the existing ones
)
Calling fit() removes the stored out-of-sample residuals, so they must be set again every time the forecaster is retrained.
⚠ Warning
The correct estimation of prediction intervals with conformal methods depends on the residuals being representative of future errors. For this reason, calibration residuals should be used. However, the dynamics of the series and models can change over time, so it is important to monitor and regularly update the residuals. It can be done easily using the set_out_sample_residuals() method.
Key takeaways¶
The probabilistic forecasting methods of skforecast can be applied to global models: a single
ForecasterRecursiveMultiSeriesestimates intervals for every series, and each series keeps its own residuals and bins.The intervals must be built with out-of-sample residuals, obtained by backtesting on a validation (calibration) partition with the same horizon used later. In-sample residuals lead to overconfident intervals.
With
use_binned_residuals=True, the residuals are conditioned on the predicted value, so the width of the intervals adapts to the magnitude of the prediction.In this example, conformal prediction achieves a coverage close to the nominal 80% (83.95% and 83.06%), while bootstrapping produces intervals about twice as wide and overly conservative (96.29% and 97.92%), since the multi-step residuals are propagated recursively along the horizon.
Series not seen during training can also be forecast with intervals, using a random sample of the residuals of all the known series (
_unknown_level).The empirical coverage must always be validated with backtesting, and the residuals should be refreshed periodically in production and after every refit.