Metrics in probabilistic forecasting¶

In point estimate forecasting, the model outputs a single value that ideally represents the most likely value of the time series at future steps. In this scenario, the quality of the predictions can be assessed by comparing the predicted value with the true value of the series. Examples of metrics used for this purpose include Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE).

In probabilistic forecasting, however, the model does not produce a single value, but rather a representation of the entire distribution of possible predicted values. In practice, this is often represented by a sample of the underlying distribution (e.g. 50 possible predicted values) or by specific quantiles that capture most of the information in the distribution. This approach provides richer insights by allowing the creation of prediction intervals (ranges within which the true value is likely to fall).

In this context, the quality of the predictions cannot be assessed with the same metrics used for point forecasts. Instead, specific metrics are needed to evaluate different aspects of a probabilistic forecast: whether the intervals are well calibrated (they contain the true value as often as their nominal level promises), how sharp they are (narrower intervals are more informative), and how well the whole predictive distribution matches the observations.

The metrics covered in this notebook are summarized below:

Metric	Aspect evaluated	Description
Coverage	Calibration	Proportion of true values that fall within the prediction interval. Should be close to the nominal level (e.g. ~80% for an 80% interval).
Interval width	Sharpness	Average width of the prediction intervals. For a given coverage, narrower is better.
Interval area	Sharpness	Total area spanned by the intervals over the forecast horizon; a global measure of interval size.
Winkler score (interval score)	Calibration + sharpness	Interval width plus a penalty for observations that fall outside the interval. Rewards narrow intervals that still capture the true value.
Weighted Interval Score (WIS)	Calibration + sharpness	Generalizes the Winkler score to several intervals plus the median forecast. A discrete approximation of the CRPS.
CRPS	Full distribution	Distance between the predicted and the empirical cumulative distribution functions. Evaluates the entire predictive distribution.

In general, the goal is to produce prediction intervals that are as narrow as possible while still capturing the true values with the desired probability. This is a trade-off between the sharpness (width) of the intervals and their calibration (coverage of the true values). Proper scoring rules such as the Winkler score, the WIS, and the CRPS combine both aspects into a single number, which makes them especially convenient for model comparison and hyperparameter tuning.

💡 Tip

For more examples on how to use probabilistic forecasting, check out the following articles:

Libraries and data¶

In [2]:

Copied!





# 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, plot_residuals

# Modelling and Forecasting
# ==============================================================================
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from feature_engine.timeseries.forecasting import LagFeatures
from feature_engine.timeseries.forecasting import WindowFeatures
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold, backtesting_forecaster
from skforecast.preprocessing import CalendarFeatures, RollingFeatures
from skforecast.metrics import (
    calculate_coverage,
    crps_from_quantiles,
    winkler_score,
    weighted_interval_score,
)

# Warnings configuration
# ==============================================================================
import warnings
warnings.filterwarnings('once')
# 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, plot_residuals

# Modelling and Forecasting
# ==============================================================================
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from feature_engine.timeseries.forecasting import LagFeatures
from feature_engine.timeseries.forecasting import WindowFeatures
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold, backtesting_forecaster
from skforecast.preprocessing import CalendarFeatures, RollingFeatures
from skforecast.metrics import (
    calculate_coverage,
    crps_from_quantiles,
    winkler_score,
    weighted_interval_score,
)

# Warnings configuration
# ==============================================================================
import warnings
warnings.filterwarnings('once')

In [3]:

Copied!





# Load data
# ==============================================================================
data = fetch_dataset('ett_m2')
data = data.resample(rule="1h", closed="left", label="right").mean()
data.head(3)
# Load data
# ==============================================================================
data = fetch_dataset('ett_m2')
data = data.resample(rule="1h", closed="left", label="right").mean()
data.head(3)

╭───────────────────────────────────── ett_m2 ─────────────────────────────────────╮
│ Description:                                                                     │
│ Data from an electricity transformer station was collected between July 2016 and │
│ July 2018 (2 years x 365 days x 24 hours x 4 intervals per hour = 70,080 data    │
│ points). Each data point consists of 8 features, including the date of the       │
│ point, the predictive value "Oil Temperature (OT)", and 6 different types of     │
│ external power load features: High UseFul Load (HUFL), High UseLess Load (HULL), │
│ Middle UseFul Load (MUFL), Middle UseLess Load (MULL), Low UseFul Load (LUFL),   │
│ Low UseLess Load (LULL).                                                         │
│                                                                                  │
│ Source:                                                                          │
│ Zhou, Haoyi & Zhang, Shanghang & Peng, Jieqi & Zhang, Shuai & Li, Jianxin &      │
│ Xiong, Hui & Zhang, Wancai. (2020). Informer: Beyond Efficient Transformer for   │
│ Long Sequence Time-Series Forecasting.                                           │
│ [10.48550/arXiv.2012.07436](https://arxiv.org/abs/2012.07436).                   │
│ https://github.com/zhouhaoyi/ETDataset                                           │
│                                                                                  │
│ URL:                                                                             │
│ https://raw.githubusercontent.com/skforecast/skforecast-                         │
│ datasets/main/data/ETTm2.csv                                                     │
│                                                                                  │
│ Shape: 69680 rows x 7 columns                                                    │
╰──────────────────────────────────────────────────────────────────────────────────╯

Out[3]:

	HUFL	HULL	MUFL	MULL	LUFL	LULL	OT
date
2016-07-01 01:00:00	38.784501	10.88975	34.753500	8.55100	4.12575	1.26050	37.83825
2016-07-01 02:00:00	36.041249	9.44475	32.696001	7.13700	3.59025	0.62900	36.84925
2016-07-01 03:00:00	38.240000	11.41350	35.343501	9.10725	3.06000	0.31175	35.91575

In [4]:

Copied!





# Split train-validation-test
# ==============================================================================
end_train = '2017-10-01 23:59:00'
end_validation = '2018-04-03 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)})")
# Split train-validation-test
# ==============================================================================
end_train = '2017-10-01 23:59:00'
end_validation = '2018-04-03 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-07-01 01:00:00 --- 2017-10-01 23:00:00  (n=10991)
Dates validation : 2017-10-02 00:00:00 --- 2018-04-03 23:00:00  (n=4416)
Dates test       : 2018-04-04 00:00:00 --- 2018-06-26 20:00:00  (n=2013)

In [5]:

Copied!





# Plot partitions of the target series
# ==============================================================================
set_dark_theme()
plt.rcParams['lines.linewidth'] = 0.5
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(data_train['OT'], label='Train')
ax.plot(data_val['OT'], label='Validation')
ax.plot(data_test['OT'], label='Test')
ax.set_title('Oil Temperature')
ax.legend();
# Plot partitions of the target series
# ==============================================================================
set_dark_theme()
plt.rcParams['lines.linewidth'] = 0.5
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(data_train['OT'], label='Train')
ax.plot(data_val['OT'], label='Validation')
ax.plot(data_test['OT'], label='Test')
ax.set_title('Oil Temperature')
ax.legend();

No description has been provided for this image

In [6]:

Copied!





# Plot partitions after differencing
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(data_train['OT'].diff(1), label='Train')
ax.plot(data_val['OT'].diff(1), label='Validation')
ax.plot(data_test['OT'].diff(1), label='Test')
ax.set_title('Differenced Oil Temperature')
ax.legend();
# Plot partitions after differencing
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(data_train['OT'].diff(1), label='Train')
ax.plot(data_val['OT'].diff(1), label='Validation')
ax.plot(data_test['OT'].diff(1), label='Test')
ax.set_title('Differenced Oil Temperature')
ax.legend();

In [7]:

Copied!





# Calendar features (cyclical encoding)
# ==============================================================================
calendar_transformer = CalendarFeatures(
    features              = ['year', 'month', 'week', 'day_of_week', 'hour'],
    encoding              = 'cyclical',  # year is automatically ignored because it is not cyclical
    keep_original_columns = True,
)

# Lags of exogenous variables
# ==============================================================================
lag_transformer = LagFeatures(
    variables = ["HUFL", "MUFL", "MULL", "HULL", "LUFL", "LULL"],
    periods   = [1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 23, 24, 42]
)

# Rolling features for exogenous variables
# ==============================================================================
wf_transformer = WindowFeatures(
    variables = ["HUFL", "MUFL", "MULL", "HULL", "LUFL", "LULL"],
    window    = ["1D", "7D"],
    functions = ["mean", "max", "min"],
    freq      = "1h",
)

exog_transformer = make_pipeline(
                       calendar_transformer,
                       lag_transformer,
                       wf_transformer
                   )
display(exog_transformer)

data = exog_transformer.fit_transform(data)
# Remove rows with NaNs created by lag features
data = data.dropna()
display(data.head(3))
# Calendar features (cyclical encoding)
# ==============================================================================
calendar_transformer = CalendarFeatures(
    features              = ['year', 'month', 'week', 'day_of_week', 'hour'],
    encoding              = 'cyclical',  # year is automatically ignored because it is not cyclical
    keep_original_columns = True,
)

# Lags of exogenous variables
# ==============================================================================
lag_transformer = LagFeatures(
    variables = ["HUFL", "MUFL", "MULL", "HULL", "LUFL", "LULL"],
    periods   = [1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 23, 24, 42]
)

# Rolling features for exogenous variables
# ==============================================================================
wf_transformer = WindowFeatures(
    variables = ["HUFL", "MUFL", "MULL", "HULL", "LUFL", "LULL"],
    window    = ["1D", "7D"],
    functions = ["mean", "max", "min"],
    freq      = "1h",
)

exog_transformer = make_pipeline(
                       calendar_transformer,
                       lag_transformer,
                       wf_transformer
                   )
display(exog_transformer)

data = exog_transformer.fit_transform(data)
# Remove rows with NaNs created by lag features
data = data.dropna()
display(data.head(3))

Pipeline(steps=[('calendarfeatures',
                 CalendarFeatures(features=['year', 'month', 'week',
                                            'day_of_week', 'hour'])),
                ('lagfeatures',
                 LagFeatures(periods=[1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14,
                                      15, 16, 17, 18, 19, 20, 21, 23, 24, 42],
                             variables=['HUFL', 'MUFL', 'MULL', 'HULL', 'LUFL',
                                        'LULL'])),
                ('windowfeatures',
                 WindowFeatures(freq='1h', functions=['mean', 'max', 'min'],
                                variables=['HUFL', 'MUFL', 'MULL', 'HULL',
                                           'LUFL', 'LULL'],
                                window=['1D', '7D']))])

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

	HUFL	HULL	MUFL	MULL	LUFL	LULL	OT	year	month_sin	month_cos	...	MULL_window_7D_min	HULL_window_7D_mean	HULL_window_7D_max	HULL_window_7D_min	LUFL_window_7D_mean	LUFL_window_7D_max	LUFL_window_7D_min	LULL_window_7D_mean	LULL_window_7D_max	LULL_window_7D_min
date
2016-07-02 19:00:00	33.6330	9.08900	29.874751	6.95600	3.753	0.61025	28.719500	2016	-0.5	-0.866025	...	5.2945	10.177786	14.115	6.806	2.130976	4.12575	0.30375	0.088702	1.2605	0.0
2016-07-02 20:00:00	31.7065	7.72750	27.884500	5.63600	3.753	0.67425	29.103875	2016	-0.5	-0.866025	...	5.2945	10.152465	14.115	6.806	2.168698	4.12575	0.30375	0.100831	1.2605	0.0
2016-07-02 21:00:00	31.8110	7.81125	27.362000	5.18025	4.435	1.42075	29.598500	2016	-0.5	-0.866025	...	5.2945	10.097352	14.115	6.806	2.204705	4.12575	0.30375	0.113864	1.2605	0.0

3 rows × 184 columns

In [8]:

Copied!





# Lags and exogenous features
# ==============================================================================
lags = [1, 2, 3, 4, 5, 6, 9, 12, 15, 17, 20, 23, 24, 42]
exog_features = [
    'year', 'month_sin', 'month_cos', 'week_sin', 'week_cos',
    'day_of_week_sin', 'day_of_week_cos', 'hour_sin', 'hour_cos',

    'HUFL', 'HUFL_lag_1', 'HUFL_lag_12', 'HUFL_lag_13', 'HUFL_lag_15', 'HUFL_lag_19',
    'HUFL_lag_2', 'HUFL_lag_20', 'HUFL_lag_23', 'HUFL_lag_4', 'HUFL_lag_5', 'HUFL_lag_9',
    'HUFL_window_1D_mean', 

    'HULL', 'HULL_lag_1', 'HULL_lag_12', 'HULL_lag_14', 'HULL_lag_15', 'HULL_lag_2', 
    'HULL_lag_20', 'HULL_lag_21', 'HULL_lag_23', 'HULL_lag_3', 'HULL_lag_4', 
    'HULL_window_1D_mean', 

    'LUFL', 'LUFL_lag_1', 'LUFL_lag_10', 'LUFL_lag_15', 'LUFL_lag_19', 'LUFL_lag_2', 
    'LUFL_lag_20', 'LUFL_lag_23', 'LUFL_lag_3', 'LUFL_lag_4', 'LUFL_lag_5', 'LUFL_window_1D_mean',

    'LULL', 'LULL_lag_1', 'LULL_lag_12', 'LULL_lag_13', 'LULL_lag_14', 'LULL_lag_18',
    'LULL_lag_19', 'LULL_lag_20', 'LULL_lag_21', 'LULL_lag_23', 'LULL_lag_24',
    'LULL_lag_4', 'LULL_lag_5', 'LULL_lag_6', 'LULL_window_1D_max', 'LULL_window_1D_min',

    'MUFL', 'MUFL_lag_1', 'MUFL_lag_11', 'MUFL_lag_12', 'MUFL_lag_13', 'MUFL_lag_15',
    'MUFL_lag_2', 'MUFL_lag_20', 'MUFL_lag_23', 'MUFL_lag_4', 'MUFL_lag_9',
    'MUFL_window_1D_mean', 

    'MULL', 'MULL_lag_1', 'MULL_lag_11', 'MULL_lag_12', 'MULL_lag_13', 'MULL_lag_14', 
    'MULL_lag_15', 'MULL_lag_17', 'MULL_lag_19', 'MULL_lag_2', 'MULL_lag_20', 
    'MULL_lag_3', 'MULL_lag_4', 'MULL_lag_5', 'MULL_lag_9', 'MULL_window_1D_mean', 
    'MULL_window_1D_min', 'MULL_window_7D_mean'
]
# Lags and exogenous features
# ==============================================================================
lags = [1, 2, 3, 4, 5, 6, 9, 12, 15, 17, 20, 23, 24, 42]
exog_features = [
    'year', 'month_sin', 'month_cos', 'week_sin', 'week_cos',
    'day_of_week_sin', 'day_of_week_cos', 'hour_sin', 'hour_cos',

    'HUFL', 'HUFL_lag_1', 'HUFL_lag_12', 'HUFL_lag_13', 'HUFL_lag_15', 'HUFL_lag_19',
    'HUFL_lag_2', 'HUFL_lag_20', 'HUFL_lag_23', 'HUFL_lag_4', 'HUFL_lag_5', 'HUFL_lag_9',
    'HUFL_window_1D_mean', 

    'HULL', 'HULL_lag_1', 'HULL_lag_12', 'HULL_lag_14', 'HULL_lag_15', 'HULL_lag_2', 
    'HULL_lag_20', 'HULL_lag_21', 'HULL_lag_23', 'HULL_lag_3', 'HULL_lag_4', 
    'HULL_window_1D_mean', 

    'LUFL', 'LUFL_lag_1', 'LUFL_lag_10', 'LUFL_lag_15', 'LUFL_lag_19', 'LUFL_lag_2', 
    'LUFL_lag_20', 'LUFL_lag_23', 'LUFL_lag_3', 'LUFL_lag_4', 'LUFL_lag_5', 'LUFL_window_1D_mean',

    'LULL', 'LULL_lag_1', 'LULL_lag_12', 'LULL_lag_13', 'LULL_lag_14', 'LULL_lag_18',
    'LULL_lag_19', 'LULL_lag_20', 'LULL_lag_21', 'LULL_lag_23', 'LULL_lag_24',
    'LULL_lag_4', 'LULL_lag_5', 'LULL_lag_6', 'LULL_window_1D_max', 'LULL_window_1D_min',

    'MUFL', 'MUFL_lag_1', 'MUFL_lag_11', 'MUFL_lag_12', 'MUFL_lag_13', 'MUFL_lag_15',
    'MUFL_lag_2', 'MUFL_lag_20', 'MUFL_lag_23', 'MUFL_lag_4', 'MUFL_lag_9',
    'MUFL_window_1D_mean', 

    'MULL', 'MULL_lag_1', 'MULL_lag_11', 'MULL_lag_12', 'MULL_lag_13', 'MULL_lag_14', 
    'MULL_lag_15', 'MULL_lag_17', 'MULL_lag_19', 'MULL_lag_2', 'MULL_lag_20', 
    'MULL_lag_3', 'MULL_lag_4', 'MULL_lag_5', 'MULL_lag_9', 'MULL_window_1D_mean', 
    'MULL_window_1D_min', 'MULL_window_7D_mean'
]

✏️ Note

The lags, features and hyperparameters used in this document were selected after a hyperparameter optimization and feature selection process. For more details, visit the full document Probabilistic forecasting: prediction intervals for multi-step time series forecasting.

Forecaster and out-of-sample residuals¶

Prediction intervals are built by bootstrapping the residuals of the forecaster. For the intervals to be reliable, these residuals must reflect the model's out-of-sample error, so they are obtained by backtesting the forecaster on a validation partition that was not used during training. The resulting residuals are then stored inside the forecaster and later resampled to generate the intervals on the test set.

💡 Tip

This section only summarizes how the out-of-sample residuals are obtained. For a detailed explanation, including binned residuals and how they are used to build prediction intervals, see the Bootstrapped residuals user guide.

In [9]:

Copied!





# Create forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=24)
forecaster = ForecasterRecursive(
                 estimator       = Ridge(random_state=15926, alpha=1.1075),
                 lags            = lags,
                 window_features = window_features,
                 differentiation = 1,
                 binner_kwargs   = {'n_bins': 10}
             )
# Create forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=24)
forecaster = ForecasterRecursive(
                 estimator       = Ridge(random_state=15926, alpha=1.1075),
                 lags            = lags,
                 window_features = window_features,
                 differentiation = 1,
                 binner_kwargs   = {'n_bins': 10}
             )

In [10]:

Copied!





# Backtesting on validation data to obtain out-of-sample residuals
# ==============================================================================
cv = TimeSeriesFold(
         initial_train_size = len(data.loc[:end_train, :]),
         steps              = 24,  # all hours of next day
         differentiation    = 1,
     )

metric_val, predictions_val = backtesting_forecaster(
                                  forecaster = forecaster,
                                  y          = data.loc[:end_validation, 'OT'],
                                  exog       = data.loc[:end_validation, exog_features],
                                  cv         = cv,
                                  metric     = 'mean_absolute_error'
                              )

display(metric_val)
fig, ax = plt.subplots(figsize=(8, 3))
data.loc[end_train:end_validation, 'OT'].plot(ax=ax, label='real value')
predictions_val['pred'].plot(ax=ax, label='prediction')
ax.set_title("Backtesting on validation data")
ax.legend();
# Backtesting on validation data to obtain out-of-sample residuals
# ==============================================================================
cv = TimeSeriesFold(
         initial_train_size = len(data.loc[:end_train, :]),
         steps              = 24,  # all hours of next day
         differentiation    = 1,
     )

metric_val, predictions_val = backtesting_forecaster(
                                  forecaster = forecaster,
                                  y          = data.loc[:end_validation, 'OT'],
                                  exog       = data.loc[:end_validation, exog_features],
                                  cv         = cv,
                                  metric     = 'mean_absolute_error'
                              )

display(metric_val)
fig, ax = plt.subplots(figsize=(8, 3))
data.loc[end_train:end_validation, 'OT'].plot(ax=ax, label='real value')
predictions_val['pred'].plot(ax=ax, label='prediction')
ax.set_title("Backtesting on validation data")
ax.legend();

  0%|          | 0/184 [00:00<?, ?it/s]

	mean_absolute_error
0	2.385951

In [11]:

Copied!





# Out-of-sample residuals distribution
# ==============================================================================
residuals = data.loc[predictions_val.index, 'OT'] - predictions_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(residuals=residuals, figsize=(7, 4))
# Out-of-sample residuals distribution
# ==============================================================================
residuals = data.loc[predictions_val.index, 'OT'] - predictions_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(residuals=residuals, figsize=(7, 4))

positive    2461
negative    1955
Name: count, dtype: int64

In [12]:

Copied!





# Store out-of-sample residuals in the forecaster
# ==============================================================================
forecaster.fit(y=data.loc[:end_train, 'OT'], exog=data.loc[:end_train, exog_features])
forecaster.set_out_sample_residuals(
    y_true = data.loc[predictions_val.index, 'OT'], 
    y_pred = predictions_val['pred']
)
# Store out-of-sample residuals in the forecaster
# ==============================================================================
forecaster.fit(y=data.loc[:end_train, 'OT'], exog=data.loc[:end_train, exog_features])
forecaster.set_out_sample_residuals(
    y_true = data.loc[predictions_val.index, 'OT'], 
    y_pred = predictions_val['pred']
)

Metrics for a single interval¶

With the out-of-sample residuals stored, an 80% prediction interval (bounded by the 10th and 90th percentiles) is generated on the test set through backtesting. The following metrics evaluate the quality of this interval from two complementary angles: how often it captures the true value (calibration) and how narrow it is (sharpness).

In [13]:

Copied!





# Backtesting with prediction intervals in test data using out-of-sample residuals
# ==============================================================================
cv = TimeSeriesFold(
         initial_train_size = len(data.loc[:end_validation, :]),
         steps              = 24,  # all hours of next day
         differentiation    = 1
     )

metric, predictions = backtesting_forecaster(
                          forecaster              = forecaster,
                          y                       = data['OT'],
                          exog                    = data[exog_features],
                          cv                      = cv,
                          metric                  = 'mean_absolute_error',
                          interval                = [0.1, 0.9],  # 80% prediction interval
                          interval_method         = 'bootstrapping',
                          n_boot                  = 150,
                          use_in_sample_residuals = False,  # Use out-of-sample residuals
                          use_binned_residuals    = True
                      )
display(metric)
predictions.head(5)
# Backtesting with prediction intervals in test data using out-of-sample residuals
# ==============================================================================
cv = TimeSeriesFold(
         initial_train_size = len(data.loc[:end_validation, :]),
         steps              = 24,  # all hours of next day
         differentiation    = 1
     )

metric, predictions = backtesting_forecaster(
                          forecaster              = forecaster,
                          y                       = data['OT'],
                          exog                    = data[exog_features],
                          cv                      = cv,
                          metric                  = 'mean_absolute_error',
                          interval                = [0.1, 0.9],  # 80% prediction interval
                          interval_method         = 'bootstrapping',
                          n_boot                  = 150,
                          use_in_sample_residuals = False,  # Use out-of-sample residuals
                          use_binned_residuals    = True
                      )
display(metric)
predictions.head(5)

  0%|          | 0/84 [00:00<?, ?it/s]

	mean_absolute_error
0	2.88077

Out[13]:

	pred	lower_bound	upper_bound
2018-04-04 00:00:00	32.750911	32.310777	33.372130
2018-04-04 01:00:00	31.954627	31.100322	33.223175
2018-04-04 02:00:00	31.065174	29.903329	33.082871
2018-04-04 03:00:00	30.126560	28.505588	32.792591
2018-04-04 04:00:00	29.252528	27.937659	32.638833

In [14]:

Copied!





# Plot intervals
# ==============================================================================
plt.rcParams['lines.linewidth'] = 1
fig, ax = plt.subplots(figsize=(9, 4))
plot_prediction_intervals(
    predictions     = predictions,
    y_true          = data_test,
    target_variable = "OT",
    initial_x_zoom  = None,
    title           = "Prediction interval in test data",
    xaxis_title     = "Date time",
    yaxis_title     = "OT",
    ax              = ax
)
fill_between_obj = ax.collections[0]
fill_between_obj.set_facecolor('white')
fill_between_obj.set_alpha(0.3)
# Plot intervals
# ==============================================================================
plt.rcParams['lines.linewidth'] = 1
fig, ax = plt.subplots(figsize=(9, 4))
plot_prediction_intervals(
    predictions     = predictions,
    y_true          = data_test,
    target_variable = "OT",
    initial_x_zoom  = None,
    title           = "Prediction interval in test data",
    xaxis_title     = "Date time",
    yaxis_title     = "OT",
    ax              = ax
)
fill_between_obj = ax.collections[0]
fill_between_obj.set_facecolor('white')
fill_between_obj.set_alpha(0.3)

Coverage, interval width and area¶

Coverage is the proportion of true values that fall inside the interval; it should be close to the nominal level (80% here). Interval width and interval area measure sharpness: for a given coverage, narrower intervals (smaller width and area) are more informative.

In [15]:

Copied!





# Empirical interval coverage (on test data)
# ==============================================================================
coverage = calculate_coverage(
               y_true      = data.loc[end_validation:, 'OT'],
               lower_bound = predictions["lower_bound"], 
               upper_bound = predictions["upper_bound"]
           )
print(f"Predicted interval coverage: {round(100 * coverage, 2)} %")

# Mean width and area of the interval
# ==============================================================================
area = (predictions["upper_bound"] - predictions["lower_bound"]).sum()
mean_width = (predictions["upper_bound"] - predictions["lower_bound"]).mean()
print(f"Area of the interval: {round(area, 2)}")
print(f"Mean width of the interval: {round(mean_width, 2)}")
# Empirical interval coverage (on test data)
# ==============================================================================
coverage = calculate_coverage(
               y_true      = data.loc[end_validation:, 'OT'],
               lower_bound = predictions["lower_bound"], 
               upper_bound = predictions["upper_bound"]
           )
print(f"Predicted interval coverage: {round(100 * coverage, 2)} %")

# Mean width and area of the interval
# ==============================================================================
area = (predictions["upper_bound"] - predictions["lower_bound"]).sum()
mean_width = (predictions["upper_bound"] - predictions["lower_bound"]).mean()
print(f"Area of the interval: {round(area, 2)}")
print(f"Mean width of the interval: {round(mean_width, 2)}")

Predicted interval coverage: 82.27 %
Area of the interval: 19950.97
Mean width of the interval: 9.91

Winkler score¶

Coverage, width, and area must be looked at together: an interval can reach perfect coverage simply by being extremely wide. The Winkler score (or interval score) combines both aspects into a single value. For each observation it takes the interval width and adds a penalty, scaled by the significance level $α$ , whenever the true value falls outside the interval:

$W_{α} (l, u, y) = {\begin{cases} (u - l) + \frac{2}{α} (l - y) & if y < l \\ (u - l) & if l \leq y \leq u \\ (u - l) + \frac{2}{α} (y - u) & if y > u \end{cases}$

where $l$ and $u$ are the lower and upper bounds, $y$ is the observed value, and $α$ is the significance level of the interval ( $α = 0.2$ for an 80% interval).

The reported value is the average across all observations. Lower is better, and because it is a proper scoring rule it can be used directly to compare or tune models.

In [ ]:

Copied!





# Winkler score (80% interval -> alpha = 0.2)
# ==============================================================================
winkler = winkler_score(
              y_true      = data.loc[end_validation:, 'OT'],
              lower_bound = predictions['lower_bound'],
              upper_bound = predictions['upper_bound'],
              alpha       = 0.2,  # (80% interval -> alpha = 0.2)
          )
print(f"Winkler score (80% interval): {round(winkler, 4)}")
# Winkler score (80% interval -> alpha = 0.2)
# ==============================================================================
winkler = winkler_score(
              y_true      = data.loc[end_validation:, 'OT'],
              lower_bound = predictions['lower_bound'],
              upper_bound = predictions['upper_bound'],
              alpha       = 0.2,  # (80% interval -> alpha = 0.2)
          )
print(f"Winkler score (80% interval): {round(winkler, 4)}")

Winkler score (80% interval): 15.1168

Metrics for multiple intervals¶

So far the evaluation has focused on a single 80% interval. A more complete assessment looks at several intervals at once: this reveals whether the model is well calibrated across the entire range of probabilities and enables distributional scores such as the Weighted Interval Score and the CRPS.

The backtesting_forecaster function can estimate many quantiles in a single run at almost no additional computational cost compared to a single interval. In the following example, a grid of quantiles is estimated and used to build prediction intervals at nominal coverage levels of 5%, 10%, ..., 95%.

In [17]:

Copied!





# Prediction intervals for different nominal coverages
# ==============================================================================
quantiles = [0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.975]
intervals = [[0.025, 0.975], [0.05, 0.95], [0.1, 0.9], [0.15, 0.85], [0.2, 0.8], [0.3, 0.7], [0.35, 0.65], [0.4, 0.6], [0.45, 0.55]]
observed_coverages = []
observed_areas = []

metric, predictions = backtesting_forecaster(
                          forecaster              = forecaster,
                          y                       = data['OT'],
                          exog                    = data[exog_features],
                          cv                      = cv,
                          metric                  = 'mean_absolute_error',
                          interval                = quantiles,
                          interval_method         = 'bootstrapping',
                          n_boot                  = 150,
                          use_in_sample_residuals = False,  # Use out-of-sample residuals
                          use_binned_residuals    = True
                      )
predictions.head()
# Prediction intervals for different nominal coverages
# ==============================================================================
quantiles = [0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.975]
intervals = [[0.025, 0.975], [0.05, 0.95], [0.1, 0.9], [0.15, 0.85], [0.2, 0.8], [0.3, 0.7], [0.35, 0.65], [0.4, 0.6], [0.45, 0.55]]
observed_coverages = []
observed_areas = []

metric, predictions = backtesting_forecaster(
                          forecaster              = forecaster,
                          y                       = data['OT'],
                          exog                    = data[exog_features],
                          cv                      = cv,
                          metric                  = 'mean_absolute_error',
                          interval                = quantiles,
                          interval_method         = 'bootstrapping',
                          n_boot                  = 150,
                          use_in_sample_residuals = False,  # Use out-of-sample residuals
                          use_binned_residuals    = True
                      )
predictions.head()

  0%|          | 0/84 [00:00<?, ?it/s]

Out[17]:

	pred	q_0.025	q_0.05	q_0.1	q_0.15	q_0.2	q_0.25	q_0.3	q_0.35	...	q_0.55	q_0.6	q_0.65	q_0.7	q_0.75	q_0.8	q_0.85	q_0.9	q_0.95	q_0.975
2018-04-04 00:00:00	32.750911	31.437079	31.902244	32.310777	32.474818	32.539450	32.600715	32.674090	32.731729	...	33.014322	33.106131	33.168963	33.196513	33.237600	33.296089	33.315526	33.372130	33.455466	33.528920
2018-04-04 01:00:00	31.954627	29.694196	30.781977	31.100322	31.380779	31.591285	31.755317	31.865466	32.011015	...	32.453554	32.564607	32.700793	32.829090	32.953904	33.061983	33.173696	33.223175	33.331078	33.416068
2018-04-04 02:00:00	31.065174	27.947510	29.222745	29.903329	30.163986	30.535393	30.813658	31.033502	31.254830	...	32.062752	32.255418	32.402917	32.591945	32.720248	32.852155	32.921553	33.082871	33.346519	33.480023
2018-04-04 03:00:00	30.126560	27.372049	27.961066	28.505588	29.168214	29.580439	29.807495	30.124945	30.401308	...	31.618759	31.823120	31.979337	32.191032	32.372613	32.523199	32.657103	32.792591	33.164071	36.074330
2018-04-04 04:00:00	29.252528	26.003375	26.735354	27.937659	28.498142	28.831877	29.229707	29.488917	29.905211	...	30.980779	31.230206	31.387892	31.588711	31.833475	32.123811	32.400639	32.638833	33.515192	36.151012

5 rows × 23 columns

Coverage and area by interval¶

The empirical coverage and area are computed for each nominal interval. Comparing the observed coverage with the nominal level shows whether the model is well calibrated across the full range of probabilities.

In [18]:

Copied!





# Calculate coverage and area for each interval
# ==============================================================================
for interval in intervals:
    observed_coverage = calculate_coverage(
                            y_true      = data.loc[end_validation:, 'OT'],
                            lower_bound = predictions[f"q_{interval[0]}"], 
                            upper_bound = predictions[f"q_{interval[1]}"]
                        )
    observed_area = (predictions[f"q_{interval[1]}"] - predictions[f"q_{interval[0]}"]).sum()
    observed_coverages.append(100 * observed_coverage)
    observed_areas.append(observed_area)

results = pd.DataFrame({
              'Interval': intervals,
              'Nominal coverage': [100 * (interval[1] - interval[0]) for interval in intervals],
              'Observed coverage': observed_coverages,
              'Area': observed_areas
          })
results.round(1)
# Calculate coverage and area for each interval
# ==============================================================================
for interval in intervals:
    observed_coverage = calculate_coverage(
                            y_true      = data.loc[end_validation:, 'OT'],
                            lower_bound = predictions[f"q_{interval[0]}"], 
                            upper_bound = predictions[f"q_{interval[1]}"]
                        )
    observed_area = (predictions[f"q_{interval[1]}"] - predictions[f"q_{interval[0]}"]).sum()
    observed_coverages.append(100 * observed_coverage)
    observed_areas.append(observed_area)

results = pd.DataFrame({
              'Interval': intervals,
              'Nominal coverage': [100 * (interval[1] - interval[0]) for interval in intervals],
              'Observed coverage': observed_coverages,
              'Area': observed_areas
          })
results.round(1)

Out[18]:

	Interval	Nominal coverage	Observed coverage	Area
0	[0.025, 0.975]	95.0	93.9	32544.9
1	[0.05, 0.95]	90.0	90.5	26129.0
2	[0.1, 0.9]	80.0	82.3	19951.0
3	[0.15, 0.85]	70.0	73.2	15990.2
4	[0.2, 0.8]	60.0	63.4	12942.2
5	[0.3, 0.7]	40.0	43.3	8067.2
6	[0.35, 0.65]	30.0	32.0	5953.4
7	[0.4, 0.6]	20.0	21.2	3893.7
8	[0.45, 0.55]	10.0	10.5	1936.2

Weighted Interval Score¶

The Winkler score evaluates a single interval. The Weighted Interval Score (WIS) extends it to a set of $K$ intervals plus the median forecast, summarizing the whole set of quantiles in one number. It is a weighted average of the individual interval scores and the absolute error of the median, and it converges to the CRPS as the number of intervals grows:

$WIS = \frac{1}{K + 0.5} (0.5 | y - m | + \sum_{k = 1}^{K} \frac{α_{k}}{2} W_{α_{k}} (l_{k}, u_{k}, y))$

where $m$ is the median forecast, $l_{k}$ and $u_{k}$ are the bounds of the $k$ -th interval, $α_{k}$ its significance level, and $W_{α_{k}}$ the Winkler score defined above. Lower is better. It is computed from the same grid of quantiles used above, taking the central 0.5 quantile as the median and the nominal intervals as the $K$ intervals.

In [21]:

Copied!





# Weighted Interval Score across the predicted intervals
# ==============================================================================
# For each interval [lower, upper] the significance level is alpha = 1 - (upper - lower)
alphas = [(1 - (upper - lower)) for lower, upper in intervals]

wis = weighted_interval_score(
          y_true       = data.loc[end_validation:, 'OT'].to_numpy(),
          y_pred       = predictions['q_0.5'].to_numpy(),  # median forecast
          lower_bounds = np.column_stack([predictions[f"q_{lower}"] for lower, upper in intervals]),
          upper_bounds = np.column_stack([predictions[f"q_{upper}"] for lower, upper in intervals]),
          alphas       = alphas,
      )
print(f"Weighted Interval Score: {round(wis, 4)}")
# Weighted Interval Score across the predicted intervals
# ==============================================================================
# For each interval [lower, upper] the significance level is alpha = 1 - (upper - lower)
alphas = [(1 - (upper - lower)) for lower, upper in intervals]

wis = weighted_interval_score(
          y_true       = data.loc[end_validation:, 'OT'].to_numpy(),
          y_pred       = predictions['q_0.5'].to_numpy(),  # median forecast
          lower_bounds = np.column_stack([predictions[f"q_{lower}"] for lower, upper in intervals]),
          upper_bounds = np.column_stack([predictions[f"q_{upper}"] for lower, upper in intervals]),
          alphas       = alphas,
      )
print(f"Weighted Interval Score: {round(wis, 4)}")

Weighted Interval Score: 2.2341

CRPS¶

The Continuous Ranked Probability Score (CRPS) measures the distance between the predicted and the empirical cumulative distribution functions, evaluating the full predictive distribution rather than a single interval. It is computed for each prediction with the function crps_from_quantiles and averaged to obtain a single value that summarizes the quality of the forecast.

In [20]:

Copied!





# Average CRPS
# ==============================================================================

# Collapse predictions to a single column (two first are excluded)
predicted_q = predictions.iloc[:, 2:].apply(
    lambda row: np.array(row), axis=1
).to_frame(name='predicted_q')
predicted_q = pd.concat([data.loc[end_validation:, 'OT'], predicted_q], axis=1)

# Calculate CRPS for each row
predicted_q['crps'] = predicted_q.apply(
    lambda row: crps_from_quantiles(
        y_true=row['OT'], pred_quantiles=row['predicted_q'], quantile_levels=np.array(quantiles)
    ), 
    axis=1
)

crps = predicted_q['crps'].mean()
print(f"Average CRPS: {round(crps, 2)}")
predicted_q.head(3)
# Average CRPS
# ==============================================================================

# Collapse predictions to a single column (two first are excluded)
predicted_q = predictions.iloc[:, 2:].apply(
    lambda row: np.array(row), axis=1
).to_frame(name='predicted_q')
predicted_q = pd.concat([data.loc[end_validation:, 'OT'], predicted_q], axis=1)

# Calculate CRPS for each row
predicted_q['crps'] = predicted_q.apply(
    lambda row: crps_from_quantiles(
        y_true=row['OT'], pred_quantiles=row['predicted_q'], quantile_levels=np.array(quantiles)
    ), 
    axis=1
)

crps = predicted_q['crps'].mean()
print(f"Average CRPS: {round(crps, 2)}")
predicted_q.head(3)

Average CRPS: 2.29

Out[20]:

	OT	predicted_q	crps
2018-04-04 00:00:00	32.674500	[31.43707893897429, 31.902244140023598, 32.310...	0.166012
2018-04-04 01:00:00	31.575750	[29.6941960913904, 30.781976558802896, 31.1003...	0.458671
2018-04-04 02:00:00	29.763125	[27.94750953838116, 29.22274495229103, 29.9033...	1.294523

	steps steps: list of tuples List of (name of step, estimator) tuples that are to be chained in sequential order. To be compatible with the scikit-learn API, all steps must define `fit`. All non-last steps must also define `transform`. See :ref:`Combining Estimators <combining_estimators>` for more details.	[('calendarfeatures', ...), ('lagfeatures', ...), ...]
	transform_input transform_input: list of str, default=None The names of the :term:`metadata` parameters that should be transformed by the pipeline before passing it to the step consuming it. This enables transforming some input arguments to ``fit`` (other than ``X``) to be transformed by the steps of the pipeline up to the step which requires them. Requirement is defined via :ref:`metadata routing <metadata_routing>`. For instance, this can be used to pass a validation set through the pipeline. You can only set this if metadata routing is enabled, which you can enable using ``sklearn.set_config(enable_metadata_routing=True)``. .. versionadded:: 1.6	None
	memory memory: str or object with the joblib.Memory interface, default=None Used to cache the fitted transformers of the pipeline. The last step will never be cached, even if it is a transformer. By default, no caching is performed. If a string is given, it is the path to the caching directory. Enabling caching triggers a clone of the transformers before fitting. Therefore, the transformer instance given to the pipeline cannot be inspected directly. Use the attribute ``named_steps`` or ``steps`` to inspect estimators within the pipeline. Caching the transformers is advantageous when fitting is time consuming. See :ref:`sphx_glr_auto_examples_neighbors_plot_caching_nearest_neighbors.py` for an example on how to enable caching.	None
	verbose verbose: bool, default=False If True, the time elapsed while fitting each step will be printed as it is completed.	False

	features	['year', 'month', ...]
	features_to_encode	None
	encoding	'cyclical'
	max_values	None
	spline_kwargs	None
	keep_original_columns	True
	tol	1e-12

	variables	['HUFL', 'MUFL', ...]
	periods	[1, 2, ...]
	freq	None
	fill_value	None
	sort_index	True
	missing_values	'raise'
	drop_original	False
	drop_na	False

	variables	['HUFL', 'MUFL', ...]
	window	['1D', '7D']
	functions	['mean', 'max', ...]
	freq	'1h'
	min_periods	None
	periods	1
	sort_index	True
	missing_values	'raise'
	drop_original	False
	drop_na	False