Probabilistic Forecasting with aeon¶

Overview of this notebook¶

• quick start - probabilistic forecasting
• disambiguation - types of probabilistic forecasts
• details: probabilistic forecasting interfaces
• metrics for, and evaluation of probabilistic forecasts
• advanced composition: pipelines, tuning, reduction
• extender guide

Quick Start - Probabilistic Forecasting with aeon¶

... works exactly like the basic forecasting workflow, replace predict by a probabilistic method!

probabilistic forecasting methods in aeon:

• forecast intervals - predict_interval(fh=None, X=None, coverage=0.90)
• forecast quantiles - predict_quantiles(fh=None, X=None, alpha=[0.05, 0.95])
• forecast variance - predict_var(fh=None, X=None, cov=False)
• distribution forecast - predict_proba(fh=None, X=None, marginal=True)

To check which forecasters in aeon support probabilistic forecasting, use the registry.all_estimators utility and search for estimators which have the capability:pred_int tag (value True).

For composites such as pipelines, a positive tag means that logic is implemented if (some or all) components support it.

What is probabilistic forecasting?¶

Intuition¶

• produce low/high scenarios of forecasts
• quantify uncertainty around forecasts
• produce expected range of variation of forecasts

Interface view¶

Want to produce "distribution" or "range" of forecast values,

at time stamps defined by forecasting horizon fh

given past data y (series), and possibly exogeneous data X

Input, to fit or predict: fh, y, X

Output, from predict_probabilistic: some "distribution" or "range" object

Big caveat: there are multiple possible ways to model "distribution" or "range"!

Used in practice and easily confused! (and often, practically, confused!)

Formal view (endogeneous, one forecast time stamp)¶

Let $y(t_1), \dots, y(t_n)$ be observations at fixed time stamps $t_1, \dots, t_n$.

(we consider $y$ as an $\mathbb{R}^n$-valued random variable)

Let $y'$ be a (true) value, which will be observed at a future time stamp $\tau$.

(we consider $y'$ as an $\mathbb{R}$-valued random variable)

We have the following "types of forecasts" of $y'$:

Notes:

• different forecasters have different capabilities!
• metrics, evaluation & tuning are different by "type of forecast"
• compositors can "add" type of forecast! Example: bootstrap

More formal details & intuition:¶

• a "point forecast" is a prediction/estimate of the conditional expectation $\mathbb{E}[y'|y]$.\

Intuition: "out of many repetitions/worlds, this value is the arithmetic average of all observations".

• a "variance forecast" is a prediction/estimate of the conditional expectation $Var[y'|y]$.\

Intuition: "out of many repetitions/worlds, this value is the average squared distance of the observation to the perfect point forecast".

• a "quantile forecast", at quantile point $\alpha\in (0,1)$ is a prediction/estimate of the $\alpha$-quantile of $y'|y$, i.e., of $F^{-1}_{y'|y}(\alpha)$, where $F^{-1}$ is the (generalized) inverse cdf = quantile function of the random variable $y'|y$.\

Intuition: "out of many repetitions/worlds, a fraction of exactly $\alpha$ will have equal or smaller than this value."

• an "interval forecast" or "predictive interval" with (symmetric) coverage $c\in (0,1)$ is a prediction/estimate pair of lower bound $a$ and upper bound $b$ such that $P(a\le y' \le b| y) = c$ and $P(y' \gneq b| y) = P(y' \lneq a| y) = (1 - c) /2$.\

Intuition: "out of many repetitions/worlds, a fraction of exactly $c$ will be contained in the interval $[a,b]$, and being above is equally likely as being below".

• a "distribution forecast" or "full probabilistic forecast" is a prediction/estimate of the distribution of $y'|y$, e.g., "it's a normal distribution with mean 42 and variance 1".\

Intuition: exhaustive description of the generating mechanism of many repetitions/worlds.

Notes:

• lower/upper of interval forecasts are quantile forecasts at quantile points $0.5 - c/2$ and $0.5 + c/2$ (as long as forecast distributions are absolutely continuous).
• all other forecasts can be obtained from a full probabilistic forecasts; a full probabilistic forecast can be obtained from all quantile forecasts or all interval forecasts.
• there is no exact relation between the other types of forecasts (point or variance vs quantile)
• in particular, point forecast does not need to be median forecast aka 0.5-quantile forecast. Can be $\alpha$-quantile for any $\alpha$!

Frequent confusion in literature & python packages:

• coverage c vs quantile \alpha
• coverage c vs significance p = 1-c
• quantile of lower interval bound, p/2, vs p
• interval forecasts vs related, but substantially different concepts: confidence interval on predictive mean; Bayesian posterior or credibility interval of the predictive mean
• all forecasts above can be Bayesian, confusion: "posteriors are different" or "have to be evaluted differently"

Probabilistic forecasting interfaces in aeon¶

This section:

• walkthrough of probabilistic predict methods
• use in update/predict workflow
• multivariate and hierarchical data

All forecasters with tag capability:pred_int provide the following:

• forecast intervals - predict_interval(fh=None, X=None, coverage=0.90)
• forecast quantiles - predict_quantiles(fh=None, X=None, alpha=[0.05, 0.95])
• forecast variance - predict_var(fh=None, X=None, cov=False)
• distribution forecast - predict_proba(fh=None, X=None, marginal=True)

Generalities:

• methods do not change state, multiple can be called
• fh is optional, if passed late
• exogeneous data X can be passed

predict_interval - interval predictions¶

Inputs:
fh - forecasting horizon (not necessary if seen in fit)
coverage, float or list of floats, default=0.9
nominal coverage(s) of the prediction interval(s) queried

Output: pandas.DataFrame
Row index is fh
Column has multi-index:
1st level = variable name from y in fit
2nd level = coverage fractions in coverage
3rd level = string "lower" or "upper"\

Entries = forecasts of lower/upper interval at nominal coverage in 2nd lvl, for var in 1st lvl, for time in row

pretty-plotting the predictive interval forecasts:

multiple coverages:

predict_quantiles - quantile forecasts¶

Inputs:
fh - forecasting horizon (not necessary if seen in fit)
alpha, float or list of floats, default = [0.1, 0.9]
quantile points at which quantiles are queried

Output: pandas.DataFrame
Row index is fh
Column has multi-index:
1st level = variable name from y in fit
2nd level = quantile points in alpha\

Entries = forecasts of quantiles at quantile point in 2nd lvl, for var in 1st lvl, for time in row

pretty-plotting the quantile interval forecasts:

predict_var - variance forecasts¶

Inputs:
fh - forecasting horizon (not necessary if seen in fit)
cov, boolean, default=False
whether covariance forecasts should also be returned (not all estimators support this)

Output: pandas.DataFrame, for cov=False:
Row index is fh
Column is equal to column index of y (variables)

Entries = variance forecast for variable in col, for time in row

with covariance, using a forecaster which can return covariance forecasts:

return is pandas.DataFrame with fh indexing rows and columns;
entries are forecast covariance between row and column time
(diagonal = forecast variances)

Using probabilistic forecasts with update/predict stream workflow¶

Example:

• data observed monthly
• make probabilistic forecasts for an entire year ahead
• update forecasts every month
• start in Dec 1950

repeat the same commands with further data batches:

... and so on.

predict_proba - distribution forecasts aka "full" probabilistic forecasts#¶

Inputs:
fh - forecasting horizon (not necessary if seen in fit)
marginal, bool, optional, default=True
whether returned distribution is marginal over time points (True), or joint over time points (False)
(not all forecasters support marginal=False)

Output: scipy Distribution object. if marginal=True: batch shape 1D, len(fh) (time); event shape 1D, len(y.columns) (variables)
if marginal=False: event shape 2D, [len(fh), len(y.columns)]

a note on consistence of methods¶

Outputs of predict_interval, predict_quantiles, predict_var, predict_proba are typically but not necessarily consistent with each other!

Consistency is a weak interface requirement but not strictly enforced.

Probabilistic forecasting for multivariate and hierarchical data¶

multivariate data: first column index for different variables

hierarchical data: probabilistic forecasts per level are row-concatenated with a row hierarchy index

(more about this in the hierarchical forecasting notebook)

Metrics for probabilistic forecasts and evaluation¶

overview - theory¶

Predicted y has different form from true y, so metrics have form

metric(y_true: series, y_pred: proba_prediction) -> float

where proba_prediction is the type of the specific "probabilistic prediction type".

I.e., we have the following function signature for a loss/metric $L$:

metrics: general signature and averaging¶

intro using the example of the quantile loss aka interval loss aka pinball loss, in the univariate case.

For one quantile value $\alpha$, the (per-sample) pinball loss function is defined as
$L_{\alpha}(\widehat{y}, y) := \alpha \cdot \Theta (y - \widehat{y}) + (1-\alpha) \cdot \Theta (\widehat{y} - y)$,
where $\Theta (x) := [1$ if $x\ge 0$ and $0$ otherwise $]$ is the Heaviside function.

This can be used to evaluate:

• multiple quantile forecasts $\widehat{\bf y}:=\widehat{y}_1, \dots, \widehat{y}_k$ for quantiles ${\bf \alpha} = \alpha_1,\dots, \alpha_k$ via\

$L_{{\bf \alpha}}(\widehat{\bf y}, y) := \frac{1}{k}\sum_{i=1}^k L_{\alpha_i}(\widehat{y}_i, y)$

• interval forecasts $[\widehat{a}, \widehat{b}]$ at symmetric coverage $c$ via\

$L_c([\widehat{a},\widehat{b}], y) := \frac{1}{2} L_{\alpha_{low}}(\widehat{a}, y) + \frac{1}{2}L_{\alpha_{high}}(\widehat{b}, y)$ where $\alpha_{low} = \frac{1-c}{2}, \alpha_{high} = \frac{1+c}{2}$

(all are known to be strictly proper losses for their respective prediction object)

There are three things we can choose to average over:

• quantile values, if multiple are predicted - elements of alpha in predict_interval(fh, alpha)
• time stamps in the forecasting horizon fh - elements of fh in fit(fh) resp predict_interval(fh, alpha)
• variables in y, in case of multivariate (later, first we look at univariate)

We will show quantile values and time stamps first:

1. averaging by fh time stamps only -> one number per quantile value in alpha

2. averaging over nothing -> one number per quantile value in alpha and fh time stamp

3. averaging over both fh and quantile values in alpha -> one number

first, generating some quantile predictions. pred_quantiles now contains quantile forecasts
formally, forecasts $\widehat{y}_j(t_i)$ where $\widehat{y_j}$ are forecasts at quantile $\alpha_j$, with range $i=1\dots N, j=1\dots k$
$\alpha_j$ are the elements of alpha, and $t_i$ are the future time stamps indexed by fh

1. computing the loss by quantile point or interval end, averaged over fh time stamps\

i.e., $\frac{1}{N} \sum_{i=1}^N L_{\alpha}(\widehat{y}(t_i), y(t_i))$ for $t_i$ in the fh, and every alpha, this is one number per quantile value in alpha

2. computing the the individual loss values, by sample, no averaging,\

i.e., $L_{\alpha}(\widehat{y}(t_i), y(t_i))$ for every $t_i$ in fh and every $\alpha$ in alpha
this is one number per quantile value $\alpha$ in alpha and time point $t_i$ in fh

3. computing the loss for a multiple quantile forecast, averaged over fh time stamps and quantile values alpha\

i.e., $\frac{1}{Nk} \sum_{j=1}^k\sum_{i=1}^N L_{\alpha_j}(\widehat{y_j}(t_i), y(t_i))$ for $t_i$ in fh, and quantile values $\alpha_j$,
this is a single number that can be used in tuning (e.g., grid search) or evaluation overall

4. computing the loss for a multiple quantile forecast, averaged quantile values alpha, for individual time stamps\

i.e., $\frac{1}{k} \sum_{j=1}^k L_{\alpha_j}(\widehat{y_j}(t_i), y(t_i))$ for $t_i$ in fh, and quantile values $\alpha_j$,
this is a univariate time series at fh times $t_i$, it can be used for tuning or evaluation by horizon index

Question: why is score_average a constructor flag, and evaluate_by_index a method?

• not all losses are "by index", so evaluate_by_index logic can vary (e.g., pseudo-samples)
• constructor args define "mathematical object" of scientific signature: series -> non-temporal object\

methods define action or "way to apply", e.g., as used in tuning or reporting

Compare score_average to multioutput arg in scikit-learn metrics and aeon.

metrics: interval vs quantile metrics¶

Interval and quantile metrics can be used interchangeably:

internally, these are easily convertible to each other
recall: lower/upper interval = quantiles at $\frac{1}{2} \pm \frac{1}2$ coverage

loss object recognizes an input type automatically and computes the corresponding interval loss

evaluation by backtesting¶

• each row is one train/test split in the walkforward setting
• first col is the loss on the test fold
• last two columns summarise length of training window, cutoff between train/test

roadmap items:

implementing further metrics

• distribution prediction metrics - may need tfp extension
• advanced evaluation set-ups
• variance loss

contributions are appreciated!

Advanced composition: pipelines, tuning, reduction, adding proba forecasts to any estimator¶

composition = constructing "composite" estimators out of multiple "component" estimators

• reduction = building estimator type A using estimator type B
  • special case: adding proba forecasting capability to non-proba forecaster
  • special case: using proba supervised learner for proba forecasting
• pipelining = chaining estimators, here: transformers to a forecaster
• tuning = automated hyperparameter fitting, usually via internal evaluation loop
  • special case: grid parameter search and random parameter search tuning
  • special case: "Auto-ML", optimising not just estimator hyperparameter but also choice of estimator

Adding probabilistic forecasts to non-probabilistic forecasters¶

start with a forecaster that does not produce probabilistic predictions:

adding probabilistic predictions is possible via reduction wrappers:

the composite can now be used like any probabilistic forecaster:

roadmap items:

more compositors to enable probabilistic forecasting

• bootstrap forecast intervals
• reduction to probabilistic supervised learning
• popular "add probabilistic capability" wrappers

contributions are appreciated!

Tuning and AutoML¶

tuning and autoML with probabilistic forecasters works exactly like with "ordinary" forecasters
via ForecastingGridSearchCV or ForecastingRandomSearchCV

• change metric to tune to a probabilistic metric
• use a corresponding probabilistic metric or loss function

Internally, evaluation will be done using probabilistic metric, via backtesting evaluation.

important: to evaluate the tuned estimator, use it in evaluate or a separate benchmarking workflow.
Using internal metric/loss values amounts to in-sample evaluation, which is over-optimistic.

inspect hyperparameter fit obtained by tuning:

obtain predictions:

for AutoML, use the MultiplexForecaster to select among multiple forecasters:

Pipelines with probabilistic forecasters¶

aeon pipelines are compatible with probabilistic forecasters:

• ForecastingPipeline applies transformers to the exogeneous X argument before passing them to the forecaster
• TransformedTargetForecaster transforms y and back-transforms forecasts, including interval or quantile forecasts

ForecastingPipeline takes a chain of transformers and forecasters, say,

[t1, t2, ..., tn, f],

where t[i] are forecasters that pre-process, and tp[i] are forecasters that postprocess

fit(y, X, fh) does:¶

X1 = t1.fit_transform(X)
X2 = t2.fit_transform(X1)
etc
X[n] = t3.fit_transform(X[n-1])\

f.fit(y=y, x=X[n])

predict_[sth](X, fh) does:¶

X1 = t1.transform(X)
X2 = t2.transform(X1)
etc
X[n] = t3.transform(X[n-1])

f.predict_[sth](X=X[n], fh)

TransformedTargetForecaster takes a chain of transformers and forecasters, say,

[t1, t2, ..., tn, f, tp1, tp2, ..., tk],

where t[i] are forecasters that pre-process, and tp[i] are forecasters that postprocess

fit(y, X, fh) does:\¶

y1 = t1.fit_transform(y)
y2 = t2.fit_transform(y1)
y3 = t3.fit_transform(y2)
etc
y[n] = t3.fit_transform(y[n-1])

f.fit(y[n])

yp1 = tp1.fit_transform(yn)
yp2 = tp2.fit_transform(yp1)
yp3 = tp3.fit_transform(yp2)
etc

predict_quantiles(y, X, fh) does:¶

y1 = t1.transform(y)
y2 = t2.transform(y1)
etc
y[n] = t3.transform(y[n-1])

y_pred = f.predict_quantiles(y[n])

y_pred = t[n].inverse_transform(y_pred)
y_pred = t[n-1].inverse_transform(y_pred)
etc
y_pred = t1.inverse_transform(y_pred)
y_pred = tp1.transform(y_pred)
y_pred = tp2.transform(y_pred)
etc
y_pred = tp[n].transform(y_pred)\

Note: the remaining proba predictions are inferred from predict_quantiles.

quick creation is also possible via the * dunder method, same pipeline:

Building your own probabilistic forecaster¶

Getting started:

• follow the "implementing estimator" developer guide
• use the advanced forecasting extension template

Extension template = python "fill-in" template with to-do blocks that allow you to implement your own, aeon-compatible forecasting algorithm.

Check estimators using check_estimator

For probabilistic forecasting:

• implement at least one of predict_quantiles, predict_interval, predict_var, predict_proba
• optimally, implement all, unless identical with defaulting behaviour as below
• if only one is implemented, others use the following defaults (in this sequence, dependent availability):
  • predict_interval uses quantiles from predict_quantiles and vice versa
  • predict_var uses variance from predict_proba, or variance of normal with IQR as obtained from predict_quantiles
  • predict_interval or predict_quantiles uses quantiles from predict_proba distribution
  • predict_proba returns normal with mean predict and variance predict_var
• so if predictive residuals not normal, implement predict_proba or predict_quantiles
• if interfacing, implement the ones where least "conversion" is necessary
• ensure to set the capability:pred_int tag to True

Summary¶

• unified API for probabilistic forecasting and probabilistic metrics
• integrating other packages (e.g. scikit-learn, statsmodels, pmdarima, prophet)
• interface for composite model building is same, proba or not (pipelining, ensembling, tuning, reduction)
• easily extensible with custom estimators

Useful resources¶

• For a good introduction to forecasting, see Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2018.
• For comparative benchmarking studies/forecasting competitions, see the M4 competition and the M5 competition.