Probabilistic Forecasting: Quantile Regression¶

Quantile regression trains a model to predict a given quantile of the target instead of its mean. Unlike bootstrapping or conformal prediction, which build intervals from the residuals of a point forecaster, quantile regression predicts the bounds of the interval directly: one model is trained for the lower quantile and another for the upper one. For example, models trained for the quantiles $\alpha = 0.1$ and $\alpha = 0.9$ produce an interval with a nominal coverage of 80% ($0.9 - 0.1 = 0.8$).

No special forecaster is needed. When the estimator is trained with a quantile loss (for example, LGBMRegressor(objective='quantile', alpha=0.1)), the predict method of the forecaster returns that quantile. Two forecasters, one per bound, are combined to obtain the interval.

This guide shows how to tune each quantile forecaster with the pinball loss, generate the intervals with backtesting, handle quantile crossing, and evaluate the coverage of the intervals against the other methods of this series.

For a continuous distribution, the $\alpha$-quantile $Q_{\alpha}(x)$ is the value below which the target $Y$ falls with probability $\alpha$, given $X = x$. For example, the 0.5-quantile is the median: half of the values are expected to be below it.

Several machine learning algorithms are capable of modeling quantiles. Some of them are:

• LightGBM

• XGBoost

• CatBoost

• Scikit-learn HistGradientBoostingRegressor

• Scikit-learn QuantileRegressor

Just as the squared-error loss function is used to train models that predict the mean value, a specific loss function is needed in order to train models that predict quantiles. The most common metric used for quantile regression is called quantile loss or pinball loss:

$$\text{pinball}(y, \hat{y}) = \frac{1}{n_{\text{samples}}} \sum_{i=0}^{n_{\text{samples}}-1} \left[ \alpha \max(y_i - \hat{y}_i, 0) + (1 - \alpha) \max(\hat{y}_i - y_i, 0) \right]$$

where $\alpha$ is the target quantile, $y$ the real value and $\hat{y}$ the quantile prediction. Note that here $\alpha$ denotes the target quantile (as in the alpha argument of LightGBM and of create_mean_pinball_loss), not the miscoverage level used in the Winkler score.

It can be seen that the loss differs depending on the evaluated quantile. The higher the quantile, the more the loss function penalizes underestimates, and the less it penalizes overestimates. As with MSE and MAE, the goal is to minimize its values (the lower the loss, the better).

Two disadvantages of quantile regression, compared to the bootstrap approach to prediction intervals, are that each quantile requires its own estimator and quantile regression is not available for all types of regression models. However, once the models are trained, the inference is much faster, since no bootstrap resampling of residuals is needed: each bound is obtained with a single prediction.

💡 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

⚠ Warning

Limitations of quantile regression in multi-step forecasting

• Quantile crossing: the two models are trained independently, so nothing prevents the predicted lower bound from being above the predicted upper bound at some steps.
• Recursive predictions: in a recursive forecaster (ForecasterRecursive), each prediction is used as a lag to predict the next step. A model trained for the quantile 0.1 feeds its own low predictions back as inputs, so, beyond the first step, its output is no longer a true quantile of the multi-step distribution. Direct strategies do not have this limitation, since predictions are never used as predictors. For this reason, this guide uses ForecasterDirect.
• Cost of the direct strategy: ForecasterDirect is slower than ForecasterRecursive because it requires training one model per step (with steps = 24, each forecaster of this guide trains 24 models, 48 in total). Although it can achieve better performance, its scalability is an important limitation when many steps need to be predicted.

For these reasons, the empirical coverage should always be validated. More details can be found in Probabilistic forecasting: prediction intervals for multi-step time series forecasting.

Libraries and data¶

The data used in this guide is the hourly number of users of the bike sharing system of Washington D.C., together with weather and holiday information (fetch_dataset with name='bike_sharing'). The goal is to predict the next 24 hours with an 80% prediction interval. The data and the train, validation and test partitions are the same as in the other probabilistic forecasting guides, so that the results can be compared on the same test set.

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

# Modelling and Forecasting
# ==============================================================================
from lightgbm import LGBMRegressor
from skforecast.direct import ForecasterDirect
from skforecast.preprocessing import CalendarFeatures, RollingFeatures
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.metrics import calculate_coverage, winkler_score
from skforecast.metrics import create_mean_pinball_loss

# 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 # Modelling and Forecasting # ============================================================================== from lightgbm import LGBMRegressor from skforecast.direct import ForecasterDirect from skforecast.preprocessing import CalendarFeatures, RollingFeatures from skforecast.model_selection import TimeSeriesFold from skforecast.model_selection import backtesting_forecaster from skforecast.model_selection import bayesian_search_forecaster from skforecast.metrics import calculate_coverage, winkler_score from skforecast.metrics import create_mean_pinball_loss # Configuration # ============================================================================== import warnings warnings.filterwarnings('once')

In [2]:

# Data download
# ==============================================================================
data = fetch_dataset(name='bike_sharing', raw=False)
data = data[['users', 'temp', 'hum', 'windspeed', 'holiday']]
data = data.loc['2011-04-01 00:00:00':'2012-10-20 23:00:00', :].copy()
data.head(3)

# Data download # ============================================================================== data = fetch_dataset(name='bike_sharing', raw=False) data = data[['users', 'temp', 'hum', 'windspeed', 'holiday']] data = data.loc['2011-04-01 00:00:00':'2012-10-20 23:00:00', :].copy() data.head(3)

╭───────────────────────────────── bike_sharing ──────────────────────────────────╮
│ Description:                                                                    │
│ Hourly usage of the bike share system in the city of Washington D.C. during the │
│ years 2011 and 2012. In addition to the number of users per hour, information   │
│ about weather conditions and holidays is available.                             │
│                                                                                 │
│ Source:                                                                         │
│ Fanaee-T,Hadi. (2013). Bike Sharing Dataset. UCI Machine Learning Repository.   │
│ https://doi.org/10.24432/C5W894.                                                │
│                                                                                 │
│ URL:                                                                            │
│ https://raw.githubusercontent.com/skforecast/skforecast-                        │
│ datasets/main/data/bike_sharing_dataset_clean.csv                               │
│                                                                                 │
│ Shape: 17544 rows x 11 columns                                                  │
╰─────────────────────────────────────────────────────────────────────────────────╯

Out[2]:

	users	temp	hum	windspeed	holiday
date_time
2011-04-01 00:00:00	6.0	10.66	100.0	11.0014	0.0
2011-04-01 01:00:00	4.0	10.66	100.0	11.0014	0.0
2011-04-01 02:00:00	7.0	10.66	93.0	12.9980	0.0

Additional features are created based on calendar information: month, week, day of the week and hour. These variables are cyclical (hour 23 is as close to hour 0 as hour 1 is), so they are encoded with sine and cosine transformations that preserve this continuity.

In [3]:

# Split train-validation-test
# ==============================================================================
end_train = '2012-06-30 23:59:00'
end_validation = '2012-10-01 23:59:00'
data_train = data.loc[: end_train, :]
data_val   = data.loc[end_train:end_validation, :]
data_test  = data.loc[end_validation:, :]

print(f"Dates train      : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Dates validation : {data_val.index.min()} --- {data_val.index.max()}  (n={len(data_val)})")
print(f"Dates test       : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")

# Split train-validation-test # ============================================================================== end_train = '2012-06-30 23:59:00' end_validation = '2012-10-01 23:59:00' data_train = data.loc[: end_train, :] data_val = data.loc[end_train:end_validation, :] data_test = data.loc[end_validation:, :] print(f"Dates train : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})") print(f"Dates validation : {data_val.index.min()} --- {data_val.index.max()} (n={len(data_val)})") print(f"Dates test : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")

Dates train      : 2011-04-01 00:00:00 --- 2012-06-30 23:00:00  (n=10968)
Dates validation : 2012-07-01 00:00:00 --- 2012-10-01 23:00:00  (n=2232)
Dates test       : 2012-10-02 00:00:00 --- 2012-10-20 23:00:00  (n=456)

In [4]:

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

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

Quantile forecasters¶

Two quantile forecasters are trained, one for the 10% quantile and another for the 90% quantile. Their predictions can be combined to generate a prediction interval with a nominal coverage of 80%.

Both use the same data, partitions and predictors as the bootstrapped residuals user guide (lags, rolling mean of the last 72 hours created with RollingFeatures, calendar features and exogenous variables), so that the intervals can be compared with those of the other methods on the same test set. There are two differences, though: the forecasters are of type ForecasterDirect instead of ForecasterRecursive (see the warning above), and the hyperparameters of each quantile model are tuned separately.

In [5]:

# Create forecasters: one for each bound of the interval
# ==============================================================================
# The forecasters obtained for alpha=0.1 and alpha=0.9 produce an 80% prediction
# interval (90% - 10% = 80%).
lags = [1, 2, 3, 23, 24, 25, 167, 168, 169]
window_features = RollingFeatures(stats=["mean"], window_sizes=24 * 3)

# Calendar features (cyclical encoding)
calendar_transformer = CalendarFeatures(
                           features = ['month', 'week', 'day_of_week', 'hour'],
                           encoding = 'cyclical'
                       )

# Exogenous variables provided by the user (calendar features are created internally)
exog_features = ['holiday', 'hum', 'temp', 'windspeed']

# Forecaster for quantile 10%
forecaster_q10 = ForecasterDirect(
                     estimator = LGBMRegressor(
                                     objective    = 'quantile',
                                     alpha        = 0.1,
                                     random_state = 15926,
                                     verbose      = -1
                                 ),
                     steps             = 24,
                     lags              = lags,
                     window_features   = window_features,
                     calendar_features = calendar_transformer
                 )

# Forecaster for quantile 90%
forecaster_q90 = ForecasterDirect(
                     estimator = LGBMRegressor(
                                     objective    = 'quantile',
                                     alpha        = 0.9,
                                     random_state = 15926,
                                     verbose      = -1
                                 ),
                     steps             = 24,
                     lags              = lags,
                     window_features   = window_features,
                     calendar_features = calendar_transformer
                 )

# Create forecasters: one for each bound of the interval # ============================================================================== # The forecasters obtained for alpha=0.1 and alpha=0.9 produce an 80% prediction # interval (90% - 10% = 80%). lags = [1, 2, 3, 23, 24, 25, 167, 168, 169] window_features = RollingFeatures(stats=["mean"], window_sizes=24 * 3) # Calendar features (cyclical encoding) calendar_transformer = CalendarFeatures( features = ['month', 'week', 'day_of_week', 'hour'], encoding = 'cyclical' ) # Exogenous variables provided by the user (calendar features are created internally) exog_features = ['holiday', 'hum', 'temp', 'windspeed'] # Forecaster for quantile 10% forecaster_q10 = ForecasterDirect( estimator = LGBMRegressor( objective = 'quantile', alpha = 0.1, random_state = 15926, verbose = -1 ), steps = 24, lags = lags, window_features = window_features, calendar_features = calendar_transformer ) # Forecaster for quantile 90% forecaster_q90 = ForecasterDirect( estimator = LGBMRegressor( objective = 'quantile', alpha = 0.9, random_state = 15926, verbose = -1 ), steps = 24, lags = lags, window_features = window_features, calendar_features = calendar_transformer )

Hyperparameter tuning with the pinball loss¶

Next, a Bayesian search (bayesian_search_forecaster) is performed to find the best hyperparameters for each quantile forecaster. When validating a quantile regression model, it is important to use a metric that is coherent with the quantile being evaluated. In this case, the pinball loss is used. Skforecast provides the function create_mean_pinball_loss to create the pinball loss of a given quantile, so that each forecaster is evaluated with the quantile it was trained for.

The validation strategy is defined with TimeSeriesFold: each candidate model is trained with the training set and evaluated on the validation set in folds of 24 steps. The test set is never used during the search. Since return_best = True by default, once the search is finished, each forecaster is refitted with the best configuration found, using all the data passed to the search (training and validation).

The number of trials (n_trials = 7) is kept small to limit the execution time of this guide. In practice, a larger number of trials is recommended.

In [6]:

# Bayesian search of hyperparameters for each quantile forecaster
# ==============================================================================
def search_space(trial):
    return {
        'n_estimators'  : trial.suggest_int('n_estimators', 100, 500, step=50),
        'max_depth'     : trial.suggest_int('max_depth', 3, 10, step=1),
        'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.1)
    }

cv = TimeSeriesFold(
         steps              = 24,
         initial_train_size = len(data.loc[:end_train]),
         refit              = False
     )

results_q10, _ = bayesian_search_forecaster(
                     forecaster   = forecaster_q10,
                     y            = data.loc[:end_validation, 'users'],
                     exog         = data.loc[:end_validation, exog_features],
                     cv           = cv,
                     metric       = create_mean_pinball_loss(alpha=0.1),
                     search_space = search_space,
                     n_trials     = 7
                 )

results_q90, _ = bayesian_search_forecaster(
                     forecaster   = forecaster_q90,
                     y            = data.loc[:end_validation, 'users'],
                     exog         = data.loc[:end_validation, exog_features],
                     cv           = cv,
                     metric       = create_mean_pinball_loss(alpha=0.9),
                     search_space = search_space,
                     n_trials     = 7
                 )

print('Best results for quantile 0.1')
display(results_q10.head(3))
print('Best results for quantile 0.9')
display(results_q90.head(3))

# Bayesian search of hyperparameters for each quantile forecaster # ============================================================================== def search_space(trial): return { 'n_estimators' : trial.suggest_int('n_estimators', 100, 500, step=50), 'max_depth' : trial.suggest_int('max_depth', 3, 10, step=1), 'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.1) } cv = TimeSeriesFold( steps = 24, initial_train_size = len(data.loc[:end_train]), refit = False ) results_q10, _ = bayesian_search_forecaster( forecaster = forecaster_q10, y = data.loc[:end_validation, 'users'], exog = data.loc[:end_validation, exog_features], cv = cv, metric = create_mean_pinball_loss(alpha=0.1), search_space = search_space, n_trials = 7 ) results_q90, _ = bayesian_search_forecaster( forecaster = forecaster_q90, y = data.loc[:end_validation, 'users'], exog = data.loc[:end_validation, exog_features], cv = cv, metric = create_mean_pinball_loss(alpha=0.9), search_space = search_space, n_trials = 7 ) print('Best results for quantile 0.1') display(results_q10.head(3)) print('Best results for quantile 0.9') display(results_q90.head(3))

Best trial: 2. Best value: 12.2885: 100%7/7 [02:30<00:00, 19.67s/it]

Best trial: 0. Best value: 15.4231: 100%7/7 [02:30<00:00, 19.31s/it]

Best results for quantile 0.1

	trial_number	lags	params	mean_pinball_loss_q	n_estimators	max_depth	learning_rate
0	2	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 500, 'max_depth': 8, 'learnin...	12.288520	500.0	8.0	0.053284
1	6	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 300, 'max_depth': 7, 'learnin...	12.416458	300.0	7.0	0.067096
2	3	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 250, 'max_depth': 5, 'learnin...	12.801295	250.0	5.0	0.075614

Best results for quantile 0.9

	trial_number	lags	params	mean_pinball_loss_q	n_estimators	max_depth	learning_rate
0	0	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 400, 'max_depth': 5, 'learnin...	15.423148	400.0	5.0	0.030417
1	3	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 250, 'max_depth': 5, 'learnin...	15.577917	250.0	5.0	0.075614
2	1	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 300, 'max_depth': 8, 'learnin...	16.209979	300.0	8.0	0.048080

The best configuration is different for each quantile: the 10% quantile forecaster selects 500 trees with a maximum depth of 8 and a learning rate of about 0.053, while the 90% quantile forecaster selects 400 trees with a maximum depth of 5 and a learning rate of about 0.030. This illustrates why each bound benefits from its own tuning. Note that the pinball losses of the two quantiles (12.29 and 15.42 in the validation set) should not be compared with each other: each one weights underestimates and overestimates differently and measures the error of a different quantile.

Predictions on the test set (backtesting)¶

Once the quantile forecasters are trained, they are used to predict each of the bounds of the interval on the test set with backtesting_forecaster. The test set (19 days) is split into 19 folds of 24 hours. Since refit = False, each forecaster is trained only once, with the training and validation data, and then predicts every fold of the test set. The predictions of the two forecasters are combined into a single DataFrame with the lower and upper bounds of the interval.

In [7]:

# Backtesting on test data
# ==============================================================================
cv = TimeSeriesFold(
         steps              = 24,
         initial_train_size = len(data.loc[:end_validation]),
         refit              = False
     )

metric_q10, predictions_q10 = backtesting_forecaster(
                                  forecaster = forecaster_q10,
                                  y          = data['users'],
                                  exog       = data[exog_features],
                                  cv         = cv,
                                  metric     = create_mean_pinball_loss(alpha=0.1)
                              )

metric_q90, predictions_q90 = backtesting_forecaster(
                                  forecaster = forecaster_q90,
                                  y          = data['users'],
                                  exog       = data[exog_features],
                                  cv         = cv,
                                  metric     = create_mean_pinball_loss(alpha=0.9)
                              )

print(f"Pinball loss of the 10% quantile forecaster (test): {metric_q10.iloc[0, 0]:.2f}")
print(f"Pinball loss of the 90% quantile forecaster (test): {metric_q90.iloc[0, 0]:.2f}")

predictions = pd.DataFrame({
    'lower_bound': predictions_q10['pred'],
    'upper_bound': predictions_q90['pred']
})
predictions.head()

# Backtesting on test data # ============================================================================== cv = TimeSeriesFold( steps = 24, initial_train_size = len(data.loc[:end_validation]), refit = False ) metric_q10, predictions_q10 = backtesting_forecaster( forecaster = forecaster_q10, y = data['users'], exog = data[exog_features], cv = cv, metric = create_mean_pinball_loss(alpha=0.1) ) metric_q90, predictions_q90 = backtesting_forecaster( forecaster = forecaster_q90, y = data['users'], exog = data[exog_features], cv = cv, metric = create_mean_pinball_loss(alpha=0.9) ) print(f"Pinball loss of the 10% quantile forecaster (test): {metric_q10.iloc[0, 0]:.2f}") print(f"Pinball loss of the 90% quantile forecaster (test): {metric_q90.iloc[0, 0]:.2f}") predictions = pd.DataFrame({ 'lower_bound': predictions_q10['pred'], 'upper_bound': predictions_q90['pred'] }) predictions.head()

100%19/19 [00:00<00:00, 351.54it/s]

100%19/19 [00:00<00:00, 371.54it/s]

Pinball loss of the 10% quantile forecaster (test): 12.34
Pinball loss of the 90% quantile forecaster (test): 13.84

Out[7]:

	lower_bound	upper_bound
2012-10-02 00:00:00	35.148866	66.333425
2012-10-02 01:00:00	13.008664	29.715711
2012-10-02 02:00:00	6.298987	15.876301
2012-10-02 03:00:00	0.796944	8.647740
2012-10-02 04:00:00	5.321382	11.967840

The pinball losses in the test set (12.34 for the 10% quantile and 13.84 for the 90% quantile) are similar to or lower than those obtained in the validation set (12.29 and 15.42). However, a low pinball loss does not guarantee that the interval reaches its nominal coverage, which is evaluated in the next section.

As anticipated in the warning above, the two forecasters are trained independently, so nothing prevents the predicted 10% quantile from being above the predicted 90% quantile at some steps (quantile crossing). Before the intervals are evaluated, the number of steps affected is counted and, if any, the bounds are swapped at those steps so that the lower bound is always the smaller of the two predictions. Without this correction, functions such as winkler_score raise an error because they require the upper bound to be greater than or equal to the lower bound. In this example no step is affected, but the check should always be part of the workflow.

In [8]:

# Quantile crossing: steps where the lower bound is above the upper bound
# ==============================================================================
crossing = predictions['lower_bound'] > predictions['upper_bound']
print(f"Steps with quantile crossing: {crossing.sum()} of {len(predictions)}")

# Enforce the order of the bounds by swapping them at the steps where they cross
predictions[['lower_bound', 'upper_bound']] = np.sort(
    predictions[['lower_bound', 'upper_bound']].to_numpy(), axis=1
)
predictions.head()

# Quantile crossing: steps where the lower bound is above the upper bound # ============================================================================== crossing = predictions['lower_bound'] > predictions['upper_bound'] print(f"Steps with quantile crossing: {crossing.sum()} of {len(predictions)}") # Enforce the order of the bounds by swapping them at the steps where they cross predictions[['lower_bound', 'upper_bound']] = np.sort( predictions[['lower_bound', 'upper_bound']].to_numpy(), axis=1 ) predictions.head()

Steps with quantile crossing: 0 of 456

Out[8]:

	lower_bound	upper_bound
2012-10-02 00:00:00	35.148866	66.333425
2012-10-02 01:00:00	13.008664	29.715711
2012-10-02 02:00:00	6.298987	15.876301
2012-10-02 03:00:00	0.796944	8.647740
2012-10-02 04:00:00	5.321382	11.967840

In [9]:

# Plot intervals: whole test set and zoom ["2012-10-08", "2012-10-15"]
# ==============================================================================
fig, axs = plt.subplots(2, 1, figsize=(7, 6))
titles = ['Prediction intervals in test data', 'Prediction intervals in test data (zoom in)']
for ax, title in zip(axs, titles):
    data_test['users'].plot(ax=ax, label='Real value', color='orange')
    ax.fill_between(
        predictions.index,
        predictions['lower_bound'],
        predictions['upper_bound'],
        color  = 'gray',
        alpha  = 0.6,
        zorder = 1,
        label  = '80% prediction interval'
    )
    ax.set_xlabel('')
    ax.set_title(title)
axs[0].legend()
axs[1].set_xlim(pd.to_datetime(["2012-10-08 00:00:00", "2012-10-15 00:00:00"]))
fig.tight_layout();

# Plot intervals: whole test set and zoom ["2012-10-08", "2012-10-15"] # ============================================================================== fig, axs = plt.subplots(2, 1, figsize=(7, 6)) titles = ['Prediction intervals in test data', 'Prediction intervals in test data (zoom in)'] for ax, title in zip(axs, titles): data_test['users'].plot(ax=ax, label='Real value', color='orange') ax.fill_between( predictions.index, predictions['lower_bound'], predictions['upper_bound'], color = 'gray', alpha = 0.6, zorder = 1, label = '80% prediction interval' ) ax.set_xlabel('') ax.set_title(title) axs[0].legend() axs[1].set_xlim(pd.to_datetime(["2012-10-08 00:00:00", "2012-10-15 00:00:00"])) fig.tight_layout();

Interval evaluation¶

The intervals are evaluated on the test set with the metrics described in the metrics user guide:

• Empirical coverage (calculate_coverage): proportion of true values that fall within the interval. For an 80% interval, it should be close to 80%.

• Area: sum of the widths of all the intervals. For the same coverage, a smaller area means sharper, more informative intervals.

• Winkler score (winkler_score): the width of each interval plus a penalty of $2/\alpha$ times the distance by which the true value falls outside it, averaged over all predictions (the lower, the better). Here $\alpha$ is the miscoverage level of the interval ($\alpha = 1 - 0.8 = 0.2$ for an 80% interval), not the quantile of the forecasters.

In [10]:

# Coverage, area and Winkler score of the intervals (on test data)
# ==============================================================================
y_test = data_test['users']
coverage = calculate_coverage(
               y_true      = y_test,
               lower_bound = predictions["lower_bound"],
               upper_bound = predictions["upper_bound"]
           )
print(f"Empirical coverage of the interval: {round(100 * coverage, 2)} %")

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

# 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: {round(winkler, 2)}")

# Coverage, area and Winkler score of the intervals (on test data) # ============================================================================== y_test = data_test['users'] coverage = calculate_coverage( y_true = y_test, lower_bound = predictions["lower_bound"], upper_bound = predictions["upper_bound"] ) print(f"Empirical coverage of the interval: {round(100 * coverage, 2)} %") # Area of the interval # ============================================================================== area = (predictions["upper_bound"] - predictions["lower_bound"]).sum() print(f"Area of the interval: {round(area, 2)}") # 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: {round(winkler, 2)}")

Empirical coverage of the interval: 71.93 %
Area of the interval: 71009.82
Winkler score: 261.8

The prediction intervals generated by quantile regression achieve an empirical coverage of 71.9%, below the nominal coverage of 80%. A likely cause is that the quantiles are learned by minimizing the pinball loss on the training data: flexible models such as gradient boosting fit the tails of the training data closely, so the bounds tend to be too narrow on new data. This is the same phenomenon that makes in-sample residuals produce overoptimistic intervals in the bootstrapped residuals guide. Differences between the test period and the data used for training may also contribute.

When the coverage is not satisfactory, the intervals can be calibrated with the ConformalIntervalCalibrator transformer, as shown in the conformal calibration guide. These results correspond to a single time series and a test period of 19 days, so they should not be generalized: the empirical coverage should always be validated for each use case.

Key takeaways¶

• Quantile regression builds an interval with one model per bound, each trained with the pinball loss of its quantile (for example, 0.1 and 0.9 for an 80% interval). No residuals or bootstrapping are needed, but the estimator must support a quantile objective.

• The hyperparameters should be tuned for each quantile with a metric coherent with it, such as the pinball loss created with create_mean_pinball_loss, using validation data and never the test set.

• In multi-step forecasting, the direct strategy (ForecasterDirect) is preferable, since a recursive forecaster feeds its own quantile predictions back as lags. The cost is one model per step.

• The two models are trained independently, so quantile crossing must be checked and fixed (for example, by sorting the bounds) before the intervals are evaluated.

• The coverage of quantile regression intervals is not guaranteed: in this example it was 71.9% for a nominal 80%. The empirical coverage must always be validated with backtesting and, if needed, corrected with ConformalIntervalCalibrator.