Models

PostForecasts.jl currently provides four models for handling postprocessing of point predictions into probabilistic forecasts:

• Normal error distribution
• Conformal Prediction (and Historical Simulation)
• Isotonic Distributional Regression
• Quantile Regression

Every model belongs to the PostModel supertype. Models that work exclusively with a single point forecast as a regressor are of type UniPostModel, while models that support multiple regressors are of type MultiPostModel.

Normal error distribution

The naive model for probabilistic forecasting, which assumes normally distributed errors of point forecasts.

The predictive distribution conditional on point forecast $\hat{y}$ is a Gaussian $\mathcal{N}(\hat{y} + \mu, \sigma)$, where $\mu$ and $\sigma$ are mean and sample standard deviation of errors in the training window.

The $\tau$-th quantile conditional on $\hat{y}$ of such parameterized distribution can be obtained via an analytic expression:

\[\hat{q}_{\tau|\hat{y}} = \hat{y} + \mu + \sigma \sqrt{2} \cdot \text{erf}^{-1} (2\tau - 1).\]

Normal([type::Type{F}=Float64]; zeromean::Bool=false) where {F<:AbstractFloat}

getmean(m::Normal)

getstd(m::Normal)

Conformal prediction and historical simulation

Conformal Prediction (CP) is a machine learning framework for computing prediction intervals based on the outputs of an arbitrary point forecasting model. The implemented method corresponds to inductive conformal prediction (Papadopoulos et al., 2002) and is analogous to the approach of Kath and Ziel (2021).

In the training step, the non-conformity scores $\lambda_i$ are calculated on the training set $(\hat{y}_i, y_i)_{i\in\text{training window}}$ as $\lambda_i := |\hat{y}_i - y_i|$.

In the prediction step, the $\tau$-th quantile conditional on $\hat{y_t}$ is obtained by shifting the prediction by an appropriate sample quantile of non-conformity scores:

\[\hat{q}_{\tau|\hat{y}} = \hat{y} - \mathbf{1}_{\{\tau < 0.5\}} Q_{1 - 2\tau}(\lambda) + \mathbf{1}_{\{\tau > 0.5\}} Q_{2\tau - 1}(\lambda),\]

where $Q_{\alpha}(\lambda)$ is the $\alpha$-th sample quantile of non-conformity scores from the training window. Although the intervals in the form of $[\hat{y} - Q_{\alpha}(\lambda), \hat{y} +Q_{\alpha}(\lambda)]$ are valid $\alpha$ prediction intervals without any requirements on the underlying distribution, translating them into quantiles requires the assumption of symmetrically distributed errors.

However, it is also possible to use conformal prediction to obtain non-symmetric distributions, by using non-absolute errors $\lambda_i := \hat{y}_i - y_i$. Then, in the prediction step $\tau$-th quantile conditional on $\hat{y_t}$ is computed as:

\[\hat{q}_{\tau|\hat{y}} = \hat{y} + Q_{\tau}(\lambda),\]

the method that predates conformal prediction and is widely known as Historical Simulation (HS) (Hendricks, 1996).

CP([type::Type{F}=Float64,] n::Integer[; abs::Bool=true]) where {F<:AbstractFloat}

getscores(m::CP)

Isotonic distributional regression

Isotonic Distributional Regression (IDR; Henzi et al., 2021) has been recently introduced as a novel nonparametric method for estimating distributions that are isotonic in the regressed variable, which means that the quantiles of such distributions are non-decreasing w.r.t the regressor. In the training step, $n$ observations $(\hat{y}_i, y_i)_{i \in \text{training window}}$ are first sorted to be ascending in $\hat{y}_i$. Then, $n$ conditional distributions $\hat{F_i}(z) = \hat{F}(z|x_i)$ are obtained by solving the following min-max problem via abridged pool-adjacent-violators algorithm:

\[\hat{F}_i(z) = \min_{k=1,...,i} \max_{j=k,..,n} \frac{1}{j-k+1}\sum_{l=k}^{j} \mathbb{1}\{y_{l} < z\},\]

where $z \in (y_i)_{i \in \text{training window}}$ To obtain conditional distribution for any $\hat{y}\in\mathbb{R}$, the obtained distribution functions are interpolated

\[\hat{F}(z|\hat{y}) = \frac{\hat{y}-\hat{y}_i}{\hat{y}_{i+1} - \hat{y}_i}\hat{F}_i(z) + \frac{\hat{y}_{i+1} - \hat{y}}{\hat{y}_{i+1} - \hat{y}_i} \hat{F}_{i+1}(z),\]

If $\hat{y} < \hat{y}_1$ or $\hat{y} > \hat{y}_n$, we set $\hat{F}(z|\hat{y})$ to $\hat{F}_1(z)$ or $\hat{F}_n(z)$, respectively. Finally, since ProbcastSeries stores predictive distributions in the form of quantiles, we determine quantiles at specified levels as

\[\hat{q}_{\tau|\hat{y}} = \min\{z : \hat{F}(z|\hat{y}) \geq \tau\}.\]

The multivariate version of IDR is not supported, but ForecastSeries containing multiple forecasts can be used as input for computing ProbcastSeries. In such a case, multiple univariate IDR models are estimated and the resulting distributions functions $\hat{F}(z)$ are averaged. Since $z$ is limited to true values of the timeseries in training window, the distributions resulting from estimated IDRs are defined at the exact same points, which allows to efficiently and precisely compute the average across probability.

The implemented IDR estimation uses abridged pool-adjacent-violators algorithm introduced by Henzi et al. (2022).

IDR([type::Type{F}=Float64,] n::Integer, r::Integer) where {F<:AbstractFloat}

getcdf(m::IDR [, r])

getx(m::IDR [, r])

gety(m::IDR)

Quantile regression

Quantile Regression Averaging (QRA; Nowotarski and Weron, 2014) is a well-established method in for obtaining probabilistic forecasts of electricity prices and load. It learns conditional quantiles as linear combination of $m$ point forecasts:

\[\hat{q}_{\tau|\hat{y}^{(1)}, ..., \hat{y}^{(m)}} = \beta^{(\tau)}_0 + \beta^{(\tau)}_1\hat{y}^{(1)} + ... + \beta^{(\tau)}_m\hat{y}^{(m)}.\]

The coefficients $\beta^{(\tau)}_{0...m}$ are selected to minimize the pinball loss on the training window and estimated by solving a linear programming problem. For this task, Probcasts.jl employs JuMP.jl and HiGHS.jl packages. Different LP solvers compatible with JuMP can be used, but the constructor defaults to an open source HiGHS.

Apart from the standard QRA introduced by Nowotarski and Weron (2014), PostForecasts.jl allows to readily compute Quantile Regression Machine (QRM; Marcjasz et al., 2020) and Quantile Regression with probability or quantile averaging (Uniejewski et al., 2019). See Different flavors of quantile regression for details.

QR([type::Type{F}=Float64,] n::Integer, r::Integer, prob::Union{AbstractFloat, AbstractVector{<:AbstractFloat}}) where {F<:AbstractFloat}

getweights(m::QR)

getquantprob(m::QR)

Isotonic quantile regression

iQR(args...)

Lasso quantile regression

LassoQR([type::Type{F}=Float64,] n::Integer, r::Integer, prob::Union{AbstractFloat, AbstractVector{<:AbstractFloat}}; kwargs...) where {F<:AbstractFloat}

Smoothed quantile regression

SQR([type::Type{F}=Float64,] n::Integer, r::Integer, prob::Union{AbstractFloat, AbstractVector{<:AbstractFloat}}; kwargs...) where {F<:AbstractFloat}

Generalized Autoregressive Conditional Heteroskedasticity

Generalized Autoregressive Conditional Heteroskedasticity (GARCH; Bollerslev, 1986) is a seminal model for modeling time-varying volatility in financial markets. PostForecasts.jl implements GARCH(1,1) specification, given by:

\[\hat{\sigma}^2_{t}=\alpha\varepsilon^2_{t-1}+\beta\hat{\sigma}^2_{t-1}+\omega,\]

where $\hat{\sigma}^2_{t}$ is the conditional prediction error variance at time $t$, $\varepsilon_{t-1}$ is the forecast error at the previous time step and $\alpha$, $\beta$ and $\omega$ are the model parameters.

Parameters are estimated via maximum likelihood, assuming the standardized forecast errors follow a normal distribution, $\varepsilon_t \sim \mathcal{N}(0, \hat{\sigma}^2_t)$. To reduce computational cost and improve estimation robustness against model misspecification, the asymptotic variance $\omega$ is determined using volatility targeting:

\[\omega=\bar{\sigma}^2(1-\alpha-\beta),\]

where $\bar{\sigma}$ is the sample standard deviation of the forecast errors in the calibration window.

Optimization is handled with the NLopt library, the default algorithm is set by the hyperparameter :garch_solver, see NLopt docs for details.

Multi-Step-Ahead Forecasting

For multi-step-ahead volatility forecasting, when he future error $\varepsilon_{t}$ is unknown, the model replaces its squared value with its conditional expectation, $\mathbb{E}[\varepsilon^2_{t}]=\hat{\sigma}^2_{t}$, leading to the multi-step formula:

\[\hat{\sigma}^2_{i,t+1}=(\alpha+\beta)\hat{\sigma}^2_{t} + \omega\]

Quantile Forecasts

The quantile forecasts can be generated in three different approaches, depending on two keyword argumnents of the GARCH constructor, filter and abs.

• filter=false (default): the quantile forecasts are generated from the $\mathcal{N}(0, \hat{\sigma}^2_t)$

• filter=false and abs=false: the quantile forecasts are derived from the empirical distribution of forecast errors standardized by the estimated time-varying volatility and then scaled by $\hat{\sigma}_t$ (you can think of this as a hybrid of GARCH and HS)

• filter=false and abs=true: the quantile forecasts are derived from the empriical distribution of absolute forecast errors standardized by the estimated time-varying volatility and then scaled by $\hat{\sigma}_t$ (you can think of this as a hybrid of GARCH and CP)

GARCH([type::Type{F}=Float64,] n::Integer; kwargs...) where {F<:AbstractFloat}

getparams(m::GARCH)

Training and prediction

train(m, X::AbstractVecOrMat{<:Number}, Y::AbstractVector{<:Number})

predict(m, input, quantiles)

predict!(m, output, input, quantiles)