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 (for example, 150 bootstrapped predictions) 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	Function
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).	calculate_coverage
Interval width	Sharpness	Average width of the prediction intervals. For a given coverage, narrower is better.	Computed directly with pandas
Interval area	Sharpness	Total area spanned by the intervals over the forecast horizon; a global measure of interval size.	Computed directly with pandas
Winkler score (interval score)	Calibration + sharpness	Interval width plus a penalty, proportional to $2/\alpha$, for each observation that falls outside the interval, where $1 - \alpha$ is the nominal coverage. Rewards narrow intervals that still capture the true value; the lower the better.	winkler_score
Weighted Interval Score (WIS)	Calibration + sharpness	Generalizes the Winkler score to several intervals plus the median forecast. A discrete approximation of the CRPS.	weighted_interval_score
CRPS	Full distribution	Distance between the predicted cumulative distribution function and the step function at the observed value (the empirical CDF of a single observation). Evaluates the entire predictive distribution.	crps_from_quantiles, crps_from_predictions

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:

• Probabilistic forecasting with machine learning
• Probabilistic forecasting: prediction intervals for multi-step time series forecasting
• Continuous Ranked Probability Score (CRPS) in probabilistic forecasting

Libraries and data¶

In [1]:

# 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')

The data used in this guide is the ETTm2 dataset (Electricity Transformer Temperature), available through fetch_dataset. It contains measurements recorded every 15 minutes at an electricity transformer station between July 2016 and July 2018. The target variable is the oil temperature (OT), and six power load variables (HUFL, HULL, MUFL, MULL, LUFL and LULL) are used as exogenous variables. To work at an hourly resolution, the data is aggregated to hourly means.

In [2]:

# 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[2]:

	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

The data is split chronologically, never randomly, into three partitions, so that each one only contains observations later than those of the previous one:

• Train (2016-07-01 to 2017-10-01, 10,991 hours): used to fit the forecaster.

• Validation (2017-10-02 to 2018-04-03, 4,416 hours): used to obtain the out-of-sample residuals from which the prediction intervals are built.

• Test (2018-04-04 to 2018-06-26, 2,013 hours): used only to evaluate the probabilistic forecasts with the metrics described in this guide. Before predicting it, the forecaster is retrained with the training and validation data.

In [3]:

# 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 [4]:

# 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();

The level of the oil temperature changes over time, so the series is not stationary. After differencing (subtracting the previous value from each observation), the series fluctuates around zero with a stable level in all three partitions. For this reason, the forecaster is created later with differentiation=1: it models the differenced series and reverts the transformation automatically when predicting. For more details, see the time series differentiation user guide.

In [5]:

# 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();

Three groups of exogenous features are created:

• Calendar features: year, month, week, day of the week and hour, extracted from the datetime index with CalendarFeatures. The cyclical ones are encoded with sine and cosine transformations, so that, for example, hour 23 is as close to hour 0 as hour 1 is. The year is not cyclical and is kept without encoding.

• Lags of the exogenous variables: past values of the six load variables, created with LagFeatures from the feature-engine library.

• Rolling statistics of the exogenous variables: mean, maximum and minimum of the load variables over the previous day and week, created with WindowFeatures from feature-engine.

Lags and rolling statistics only use values prior to each timestamp, so they do not introduce future information of these variables. The first rows, which do not have enough history to compute them, contain missing values and are removed.

feature-engine is not a dependency of skforecast. If it is not available, install it with pip install feature-engine.

In [6]:

# Calendar features (cyclical encoding)
# ==============================================================================
calendar_transformer = CalendarFeatures(
    features              = ['year', 'month', 'week', 'day_of_week', 'hour'],
    encoding              = 'cyclical',  # 'year' is not cyclical, so it is kept without encoding
    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 and window 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 not cyclical, so it is kept without encoding 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 and window 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 [7]:

# 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' ]

⚠ Warning

Exogenous variables are assumed to be known in advance

The forecaster predicts 24 hours ahead using, as exogenous features, the current value of the load variables (HUFL, HULL, MUFL, MULL, LUFL and LULL) and their lags of 1 to 23 hours. When the forecast is made, these values have not yet been observed for most of the forecast horizon, so this setup implicitly assumes that they are known (or perfectly forecast) for the next 24 hours. This is acceptable to illustrate how the metrics work, but in a real deployment these variables would have to be forecast as well, and the errors and interval widths obtained here would be optimistic.

✏️ 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¶

A ForecasterRecursive is created with a Ridge regressor. As predictors, it uses 14 lags of the oil temperature (between 1 and 42 hours), the mean, minimum and maximum of the last 24 hours (RollingFeatures) and the exogenous features selected above. Since the series is not stationary, differentiation=1 is used, and binner_kwargs={'n_bins': 10} groups the residuals into 10 bins according to the predicted value (binned 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 (backtesting_forecaster) on the validation partition, which was not used during training. The resulting residuals are then stored inside the forecaster with set_out_sample_residuals() and later resampled to generate the intervals on the test set. The test set must never be used to obtain these residuals: if it were, the coverage measured on that same test set would be optimistically biased.

💡 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 [8]:

# 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 [9]:

# 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,   # must match the differentiation of the forecaster
     )

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, # must match the differentiation of the forecaster ) 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();

100%184/184 [00:00<00:00, 1948.72it/s]

	mean_absolute_error
0	2.385951

In [10]:

# 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

There are more positive than negative residuals (2,461 versus 1,955). Since the residual is the observed value minus the predicted value ($e_t = y_t - \hat{y}_t$), this means that the model tends to underestimate the oil temperature in the validation period. The bootstrapping process adds these residuals to the predictions, so this bias is carried over to the intervals, which are slightly shifted upwards with respect to the point forecast.

In [11]:

# 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 with backtesting_forecaster. The forecaster is retrained with the training and validation data, while the out-of-sample residuals stored in it are kept. 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 [12]:

# Backtesting with prediction intervals in test data using out-of-sample residuals
# ==============================================================================
y_test = data.loc[end_validation:, 'OT']  # Observed values of the test partition

cv = TimeSeriesFold(
         initial_train_size = len(data.loc[:end_validation, :]),
         steps              = 24,  # all hours of next day
         differentiation    = 1,   # must match the differentiation of the forecaster
     )

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 # ============================================================================== y_test = data.loc[end_validation:, 'OT'] # Observed values of the test partition cv = TimeSeriesFold( initial_train_size = len(data.loc[:end_validation, :]), steps = 24, # all hours of next day differentiation = 1, # must match the differentiation of the forecaster ) 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)

100%84/84 [00:00<00:00, 638.65it/s]

	mean_absolute_error
0	2.88077

Out[12]:

	fold	pred	lower_bound	upper_bound
2018-04-04 00:00:00	0	32.750911	32.310777	33.372130
2018-04-04 01:00:00	0	31.954627	31.100322	33.217196
2018-04-04 02:00:00	0	31.065174	29.847008	33.013255
2018-04-04 03:00:00	0	30.126560	28.505588	32.764462
2018-04-04 04:00:00	0	29.252528	27.937659	32.588203

In [13]:

# Plot intervals
# ==============================================================================
plt.rcParams['lines.linewidth'] = 1
fig, ax = plt.subplots(figsize=(9, 4))
plot_prediction_intervals(
    predictions         = predictions,
    y_true              = y_test,
    target_variable     = "OT",
    initial_x_zoom      = None,
    title               = "Prediction interval in test data",
    xaxis_title         = "Date time",
    yaxis_title         = "OT",
    ax                  = ax,
    kwargs_fill_between = {'color': 'white', 'alpha': 0.3, 'zorder': 1}
)

# Plot intervals # ============================================================================== plt.rcParams['lines.linewidth'] = 1 fig, ax = plt.subplots(figsize=(9, 4)) plot_prediction_intervals( predictions = predictions, y_true = y_test, target_variable = "OT", initial_x_zoom = None, title = "Prediction interval in test data", xaxis_title = "Date time", yaxis_title = "OT", ax = ax, kwargs_fill_between = {'color': 'white', 'alpha': 0.3, 'zorder': 1} )

Coverage, interval width and area¶

Coverage (calculate_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 [14]:

# Empirical interval coverage (on test data)
# ==============================================================================
coverage = calculate_coverage(
               y_true      = y_test,
               lower_bound = predictions["lower_bound"],
               upper_bound = predictions["upper_bound"]
           )
print(f"Empirical coverage of the interval: {100 * coverage:.2f} %")

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

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

Empirical coverage of the interval: 82.31 %
Area of the interval: 19996.07
Mean width of the interval: 9.93

The empirical coverage (82.31%) is slightly above the nominal 80%: the interval is a little conservative, capturing the true value somewhat more often than promised. On average, the interval is 9.93 units wide (in the units of OT). The area (19,996.07) is the sum of the widths of the 2,013 hourly intervals of the test set, that is, the mean width multiplied by the number of predictions. Because it grows with the length of the test set, the area is only useful to compare models evaluated on the same test set, whereas the mean width can also be compared across test sets of different length.

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, winkler_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 $\alpha$, whenever the true value falls outside the interval:

$$W_{\alpha}(l, u, y) = \begin{cases} (u - l) + \dfrac{2}{\alpha}(l - y) & \text{if } y < l \\[4pt] (u - l) & \text{if } l \le y \le u \\[4pt] (u - l) + \dfrac{2}{\alpha}(y - u) & \text{if } y > u \end{cases}$$

where $l$ and $u$ are the lower and upper bounds, $y$ is the observed value, and $\alpha$ is the significance level of the interval ($\alpha = 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 [15]:

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

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

Winkler score (80% interval): 15.2163

The Winkler score (15.22) is higher than the mean width of the interval (9.93). The difference, about 5.3 units, is the average penalty contributed by the observations that fall outside the interval (about 17.7% of them, according to the coverage). With $\alpha = 0.2$, each of these observations adds $2/\alpha = 10$ times its distance to the nearest bound, so a few large misses have a strong impact on the score. The value has no absolute meaning by itself: it is used to compare models or configurations evaluated on the same data (the lower, the better).

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: when interval is a list of quantiles, all of them are computed from the same bootstrapped predictions and returned in one column per quantile, named q_<level> (for example, q_0.1). In the following example, 21 quantiles between 0.025 and 0.975 are estimated. Pairing the symmetric quantiles yields 10 central intervals with nominal coverage levels of 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90% and 95% (for example, q_0.1 and q_0.9 bound the 80% interval), while the 0.5 quantile is the median forecast.

In [16]:

# Prediction of multiple quantiles
# ==============================================================================
quantiles = [
    0.025, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60,
    0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 0.975
]

metric_q, predictions_q = 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_q.head()

# Prediction of multiple quantiles # ============================================================================== quantiles = [ 0.025, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 0.975 ] metric_q, predictions_q = 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_q.head()

100%84/84 [00:00<00:00, 603.47it/s]

Out[16]:

	fold	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	0	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	0	31.954627	29.694196	30.781977	31.100322	31.380779	31.591285	31.755317	31.865466	31.995785	...	32.451212	32.564607	32.683035	32.813394	32.901922	33.009584	33.126910	33.217196	33.331078	33.416068
2018-04-04 02:00:00	0	31.065174	27.947510	29.222745	29.847008	30.151470	30.518780	30.768850	31.040609	31.232733	...	31.981856	32.195147	32.342294	32.512578	32.659128	32.806228	32.915224	33.013255	33.346519	33.480023
2018-04-04 03:00:00	0	30.126560	27.372049	27.961066	28.505588	29.168214	29.549113	29.797316	30.080699	30.230310	...	31.566522	31.729480	31.957810	32.142916	32.325219	32.467567	32.581610	32.764462	33.250424	36.074330
2018-04-04 04:00:00	0	29.252528	26.003375	26.735354	27.937659	28.527867	28.793625	29.072254	29.460944	29.745509	...	30.980779	31.187488	31.355086	31.541803	31.774576	32.044779	32.362650	32.588203	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 [17]:

# Calculate coverage and area for each interval
# ==============================================================================
intervals = [
    [0.025, 0.975], [0.05, 0.95], [0.10, 0.90], [0.15, 0.85], [0.20, 0.80],
    [0.25, 0.75], [0.30, 0.70], [0.35, 0.65], [0.40, 0.60], [0.45, 0.55]
]
nominal_coverages = [100 * (upper_q - lower_q) for lower_q, upper_q in intervals]
observed_coverages = []
observed_areas = []
for lower_q, upper_q in intervals:
    lower_bound = predictions_q[f'q_{lower_q}']
    upper_bound = predictions_q[f'q_{upper_q}']
    observed_coverage = calculate_coverage(
                            y_true      = y_test,
                            lower_bound = lower_bound,
                            upper_bound = upper_bound
                        )
    observed_coverages.append(100 * observed_coverage)
    observed_areas.append((upper_bound - lower_bound).sum())

results = pd.DataFrame({
              'Interval': intervals,
              'Nominal coverage (%)': nominal_coverages,
              'Observed coverage (%)': observed_coverages,
              'Area': observed_areas
          })
results.round(2)

# Calculate coverage and area for each interval # ============================================================================== intervals = [ [0.025, 0.975], [0.05, 0.95], [0.10, 0.90], [0.15, 0.85], [0.20, 0.80], [0.25, 0.75], [0.30, 0.70], [0.35, 0.65], [0.40, 0.60], [0.45, 0.55] ] nominal_coverages = [100 * (upper_q - lower_q) for lower_q, upper_q in intervals] observed_coverages = [] observed_areas = [] for lower_q, upper_q in intervals: lower_bound = predictions_q[f'q_{lower_q}'] upper_bound = predictions_q[f'q_{upper_q}'] observed_coverage = calculate_coverage( y_true = y_test, lower_bound = lower_bound, upper_bound = upper_bound ) observed_coverages.append(100 * observed_coverage) observed_areas.append((upper_bound - lower_bound).sum()) results = pd.DataFrame({ 'Interval': intervals, 'Nominal coverage (%)': nominal_coverages, 'Observed coverage (%)': observed_coverages, 'Area': observed_areas }) results.round(2)

Out[17]:

	Interval	Nominal coverage (%)	Observed coverage (%)	Area
0	[0.025, 0.975]	95.0	93.64	32488.96
1	[0.05, 0.95]	90.0	89.87	26119.78
2	[0.1, 0.9]	80.0	82.31	19996.07
3	[0.15, 0.85]	70.0	72.83	16048.37
4	[0.2, 0.8]	60.0	64.13	13012.68
5	[0.25, 0.75]	50.0	54.20	10430.50
6	[0.3, 0.7]	40.0	42.87	8114.41
7	[0.35, 0.65]	30.0	31.74	5981.91
8	[0.4, 0.6]	20.0	21.21	3913.18
9	[0.45, 0.55]	10.0	10.73	1935.21

The observed coverage is close to the nominal level for all the intervals: slightly below it for the two widest intervals (93.6% for the 95% interval and 89.9% for the 90% interval) and slightly above it for the rest (for example, 82.3% for the 80% interval and 54.2% for the 50% interval). This indicates that the intervals are well calibrated across the whole distribution, not only for the 80% interval evaluated above. As expected, the area grows with the nominal coverage: wider intervals are the price of capturing a higher proportion of the observations.

The deviations are not symmetric, though: the central intervals (from 10% to 80%) over-cover by between 0.7 and 4.2 percentage points, so they are slightly wider than necessary, while the two widest intervals under-cover slightly, meaning that extreme observations fall beyond the tails of the predicted distribution a little more often than expected.

Weighted Interval Score¶

The Winkler score evaluates a single interval. The Weighted Interval Score (WIS) (weighted_interval_score) 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, when the quantile levels are evenly spaced, it approximates the CRPS (the approximation improves as the number of intervals grows):

$$\text{WIS} = \frac{1}{K + 0.5}\left(0.5\,|y - m| + \sum_{k=1}^{K} \frac{\alpha_k}{2}\, W_{\alpha_k}(l_k, u_k, y)\right)$$

where $m$ is the median forecast, $l_k$ and $u_k$ are the bounds of the $k$-th interval, $\alpha_k$ its significance level, and $W_{\alpha_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 10 nominal intervals as the $K$ intervals ($K = 10$).

In [18]:

# 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       = y_test,
          y_pred       = predictions_q['q_0.5'],  # median forecast
          lower_bounds = predictions_q[[f"q_{lower}" for lower, _ in intervals]],
          upper_bounds = predictions_q[[f"q_{upper}" for _, upper in intervals]],
          alphas       = alphas,
      )
print(f"Weighted Interval Score: {wis:.4f}")

# 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 = y_test, y_pred = predictions_q['q_0.5'], # median forecast lower_bounds = predictions_q[[f"q_{lower}" for lower, _ in intervals]], upper_bounds = predictions_q[[f"q_{upper}" for _, upper in intervals]], alphas = alphas, ) print(f"Weighted Interval Score: {wis:.4f}")

Weighted Interval Score: 2.2780

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. When the predictive distribution is available as a sample instead of quantiles (for example, the bootstrapped predictions returned by predict_bootstrapping()), crps_from_predictions can be used instead.

In [19]:

# Average CRPS
# ==============================================================================
quantile_levels = np.array(quantiles)
pred_quantiles = predictions_q[[f"q_{q}" for q in quantiles]].to_numpy()

# CRPS of each prediction (one row of predicted quantiles per observation)
crps_values = pd.Series(
    [
        crps_from_quantiles(
            y_true=y, pred_quantiles=row_quantiles, quantile_levels=quantile_levels
        )
        for y, row_quantiles in zip(y_test, pred_quantiles)
    ],
    index = y_test.index,
    name  = 'crps'
)

crps = crps_values.mean()
print(f"Average CRPS: {crps:.4f}")
pd.concat([y_test, crps_values], axis=1).head(3)

# Average CRPS # ============================================================================== quantile_levels = np.array(quantiles) pred_quantiles = predictions_q[[f"q_{q}" for q in quantiles]].to_numpy() # CRPS of each prediction (one row of predicted quantiles per observation) crps_values = pd.Series( [ crps_from_quantiles( y_true=y, pred_quantiles=row_quantiles, quantile_levels=quantile_levels ) for y, row_quantiles in zip(y_test, pred_quantiles) ], index = y_test.index, name = 'crps' ) crps = crps_values.mean() print(f"Average CRPS: {crps:.4f}") pd.concat([y_test, crps_values], axis=1).head(3)

Average CRPS: 2.3486

Out[19]:

	OT	crps
date
2018-04-04 00:00:00	32.674500	0.166018
2018-04-04 01:00:00	31.575750	0.455047
2018-04-04 02:00:00	29.763125	1.260602

The average CRPS (2.35) is close to the WIS (2.28), as expected, since the WIS is a discrete approximation of the CRPS computed from the same quantiles. Both are expressed in the units of the target variable, which allows a direct comparison with the mean absolute error of the point forecast (2.88): for a forecast that consists of a single value, the CRPS reduces to the absolute error. Since the CRPS of the probabilistic forecast is lower than the MAE of the point forecast, describing the uncertainty with a distribution scores better than treating the point forecast as certain.

Summary¶

Metric	Value	What it tells
Coverage (80% interval)	82.31%	Calibration: slightly above the nominal level, the interval is a little conservative.
Mean width (80% interval)	9.93	Sharpness: average size of the interval, in the units of OT.
Area (80% interval)	19,996.07	Sharpness: total size of the intervals, only comparable on the same test set.
Winkler score (80% interval)	15.2163	Width plus penalties for misses; the lower the better.
Weighted Interval Score (10 intervals + median)	2.2780	Calibration and sharpness over all the intervals; the lower the better.
CRPS (21 quantiles)	2.3486	Quality of the whole predictive distribution; the lower the better.

Coverage and width describe calibration and sharpness separately and are easy to interpret, but neither of them can rank models on its own: a wider interval always achieves a higher coverage. To compare models or tune hyperparameters, use a proper scoring rule that combines both aspects: the Winkler score when a single interval is of interest, and the WIS or the CRPS when the whole predictive distribution matters.