Prerequisites This tutorial assumes basic familiarity with StatsForecast. For a minimal example visit the Quick Start

​
Introduction

When we generate a forecast, we usually produce a single value known as the point forecast. This value, however, doesn’t tell us anything about the uncertainty associated with the forecast. To have a measure of this uncertainty, we need prediction intervals. A prediction interval is a range of values that the forecast can take with a given probability. Hence, a 95% prediction interval should contain a range of values that include the actual future value with probability 95%. Probabilistic forecasting aims to generate the full forecast distribution. Point forecasting, on the other hand, usually returns the mean or the median or said distribution. However, in real-world scenarios, it is better to forecast not only the most probable future outcome, but many alternative outcomes as well. The problem is that some timeseries models provide forecast distributions, but some other ones only provide point forecasts. How can we then estimate the uncertainty of predictions?
Prediction Intervals For models that already provide the forecast distribution, check Prediction Intervals.

​
Conformal Prediction

For a video introduction, see the PyData Seattle presentation. Multi-quantile losses and statistical models can provide provide prediction intervals, but the problem is that these are uncalibrated, meaning that the actual frequency of observations falling within the interval does not align with the confidence level associated with it. For example, a calibrated 95% prediction interval should contain the true value 95% of the time in repeated sampling. An uncalibrated 95% prediction interval, on the other hand, might contain the true value only 80% of the time, or perhaps 99% of the time. In the first case, the interval is too narrow and underestimates the uncertainty, while in the second case, it is too wide and overestimates the uncertainty. Statistical methods also assume normality. Here, we talk about another method called conformal prediction that doesn’t require any distributional assumptions. More information on the approach can be found in this repo owned by Valery Manokhin. Conformal prediction intervals use cross-validation on a point forecaster model to generate the intervals. This means that no prior probabilities are needed, and the output is well-calibrated. No additional training is needed, and the model is treated as a black box. The approach is compatible with any model. Statsforecast now supports Conformal Prediction on all available models.

​
Install libraries

We assume that you have StatsForecast already installed. If not, check this guide for instructions on how to install StatsForecast Install the necessary packages using pip install statsforecast 
pip install statsforecast -U

​
Load and explore the data

For this example, we’ll use the hourly dataset from the M4 Competition. We first need to download the data from a URL and then load it as a pandas dataframe. Notice that we’ll load the train and the test data separately. We’ll also rename the y column of the test data as y_test. 
import pandas as pd

train = pd.read_csv('https://auto-arima-results.s3.amazonaws.com/M4-Hourly.csv')
test = pd.read_csv('https://auto-arima-results.s3.amazonaws.com/M4-Hourly-test.csv').rename(columns={'y': 'y_test'})
train.head()
unique_iddsy
0H11605.0
1H12586.0
2H13586.0
3H14559.0
4H15511.0
Since the goal of this notebook is to generate prediction intervals, we’ll only use the first 8 series of the dataset to reduce the total computational time. 
n_series = 8 
uids = train['unique_id'].unique()[:n_series] # select first n_series of the dataset
train = train.query('unique_id in @uids')
test = test.query('unique_id in @uids')
We can plot these series using the plot_series function from the utilsforecast library. Thisfunctionmethod has multiple parameters, and the required ones to generate the plots in this notebook are explained below.
  • df: A pandas dataframe with columns [unique_id, ds, y].
  • forecasts_df: A pandas dataframe with columns [unique_id, ds] and models.
  • plot_random: bool = True. Plots the time series randomly.
  • models: List[str]. A list with the models we want to plot.
  • level: List[float]. A list with the prediction intervals we want to plot.
  • engine: str = matplotlib. It can also be plotly. plotly generates interactive plots, while matplotlib generates static plots.
from utilsforecast.plotting import plot_series

plot_series(train, test, plot_random=False)

​
Train models

StatsForecast can train multiple models on different time series efficiently. Most of these models can generate a probabilistic forecast, which means that they can produce both point forecasts and prediction intervals. For this example, we’ll use SimpleExponentialSmoothing and ADIDA which do not provide a prediction interval natively. Thus, it makes sense to use Conformal Prediction to generate the prediction interval. We’ll also show using it with ARIMA to provide prediction intervals that don’t assume normality. To use these models, we first need to import them from statsforecast.models and then we need to instantiate them. 
from statsforecast.models import SeasonalExponentialSmoothing, ADIDA, ARIMA
from statsforecast.utils import ConformalIntervals

# Create a list of models and instantiation parameters 
intervals = ConformalIntervals(h=24, n_windows=2)
# P.S. n_windows*h should be less than the count of data elements in your time series sequence.
# P.S. Also value of n_windows should be atleast 2 or more.

models = [
    SeasonalExponentialSmoothing(season_length=24, alpha=0.1, prediction_intervals=intervals),
    ADIDA(prediction_intervals=intervals),
    ARIMA(order=(24,0,12), prediction_intervals=intervals),
]
To instantiate a new StatsForecast object, we need the following parameters:
  • df: The dataframe with the training data.
  • models: The list of models defined in the previous step.
  • freq: A string indicating the frequency of the data. See pandas’ available frequencies.
  • n_jobs: An integer that indicates the number of jobs used in parallel processing. Use -1 to select all cores.
sf = StatsForecast(models=models, freq=1, n_jobs=-1)
Now we’re ready to generate the forecasts and the prediction intervals. To do this, we’ll use the forecast method, which takes two arguments:
  • h: An integer that represent the forecasting horizon. In this case, we’ll forecast the next 24 hours.
  • level: A list of floats with the confidence levels of the prediction intervals. For example, level=[95] means that the range of values should include the actual future value with probability 95%.
levels = [80, 90] # confidence levels of the prediction intervals 

forecasts = sf.forecast(df=train, h=24, level=levels)
forecasts.head()
unique_iddsSeasonalESSeasonalES-lo-90SeasonalES-lo-80SeasonalES-hi-80SeasonalES-hi-90ADIDAADIDA-lo-90ADIDA-lo-80ADIDA-hi-80ADIDA-hi-90ARIMAARIMA-lo-90ARIMA-lo-80ARIMA-hi-80ARIMA-hi-90
0H1701624.132703553.097423556.359139691.906266695.167983747.292568599.519220600.030467894.554670895.065916618.078274609.440076610.583304625.573243626.716472
1H1702555.698193496.653559506.833156604.563231614.742827747.292568491.669220498.330467996.2546701002.915916549.789291510.464070515.232352584.346231589.114513
2H1703514.403029462.673117464.939840563.866218566.132941747.292568475.105038475.7937911018.7913461019.480099508.099925496.574844496.990264519.209587519.625007
3H1704482.057899433.030711436.161413527.954385531.085087747.292568440.069220440.1304671054.4546701054.515916486.376622471.141813471.516997501.236246501.611431
4H1705460.222522414.270186416.959492503.485552506.174858747.292568415.805038416.1937911078.3913461078.780099470.159478445.162316446.808608493.510348495.156640

​
Plot prediction intervals

Here we’ll plot the different intervals for one timeseries. The prediction interval with the SeasonalExponentialSmoothing seen below. Even if the model generates a point forecast, we are able to get a prediction interval. The 80% prediction interval does not cross the 90% prediction interval, which is a sign that the intervals are calibrated. 
plot_series(train, forecasts, level=levels, ids=['H105'], models=['SeasonalES'])
For weaker fitting models, the conformal prediction interval can be larger. A better model corresponds to a narrower interval. 
plot_series(train, forecasts, level=levels, ids=['H105'], models=['ADIDA'])
ARIMA is an example of a model that provides a forecast distribution, but we can still use conformal prediction to generate the prediction interval. As mentioned earlier, this method has the benefit of not assuming normality. 
plot_series(train, forecasts, level=levels, ids=['H105'], models=['ARIMA'])

​
StatsForecast Object

Alternatively, the prediction interval can be defined on the StatsForecast object. This will apply to all models that don’t have the prediction_intervals defined. 
from statsforecast.models import SimpleExponentialSmoothing, ADIDA
from statsforecast.utils import ConformalIntervals
from statsforecast import StatsForecast

models = [
    SimpleExponentialSmoothing(alpha=0.1),
    ADIDA()
]

res = StatsForecast(
    models=models, 
    freq=1,
).forecast(df=train, h=24, prediction_intervals=ConformalIntervals(h=24, n_windows=2), level=[80]) 
res.head()
unique_iddsSESSES-lo-80SES-hi-80ADIDAADIDA-lo-80ADIDA-hi-80
0H1701742.669064649.221405836.116722747.292568600.030467894.554670
1H1702742.669064550.551324934.786804747.292568498.330467996.254670
2H1703742.669064523.621405961.716722747.292568475.7937911018.791346
3H1704742.669064488.121405997.216722747.292568440.1304671054.454670
4H1705742.669064464.0214051021.316722747.292568416.1937911078.391346

​
Future work

Conformal prediction has become a powerful framework for uncertainty quantification, providing well-calibrated prediction intervals without making any distributional assumptions. Its use has surged in both academia and industry over the past few years. We’ll continue working on it, and future tutorials may include:
  • Exploring larger datasets
  • Incorporating industry-specific examples
  • Investigating specialized methods like the jackknife+ that are closely related to conformal prediction (for details on the jackknife+ see here).
If you’re interested in any of these, or in any other related topic, please let us know by opening an issue on GitHub

​
Acknowledgements

We would like to thank Kevin Kho for writing this tutorial, and Valeriy Manokhin for his expertise on conformal prediction, as well as for promoting this work.

​
References

Manokhin, Valery. (2022). Machine Learning for Probabilistic Prediction. 10.5281/zenodo.6727505.