Online Distributional Regression

Distributional Modeling and Forecasting of Natural Gas Prices

Reduces estimation time by 2-3 orders of magnitude

Maintains competitive forecasting accuracy

Real-World Validation in Energy Markets

Python Package ondil on Github and PyPi

Forecasting of Day-Ahead and Month-Ahead prices

To capture the full distribution of future prices

Extensive data analysis from 2011-2020

State-Space models, skewed Student’s t distribution

Python Package sstudentt on Github and PyPi

High-Resolution Peak Demand Estimation Using Generalized Additive Models and Deep Neural Networks

rcpptimer: Rcpp Tic-Toc Timer with OpenMP Support

Predict high-resolution electricity peaks using only low-resolution data

Combine GAMs and DNN’s for superior accuracy

Won Western Power Distribution Competition

Won Best-Student-Presentation Award

R Package rcpptimer on Github and CRAN

Provides Rcpp bindings for cpptimer

tic-toc class for timing C++ code

Supports nested, overlapping and scoped timers

Supports OpenMP parallelism

//[[Rcpp::depends(rcpptimer)]]
#include <rcpptimer.h>

//[[Rcpp::export]]
void ex()
{
  Rcpp::Timer timer;
  timer.tic();
  // Code to be timed
  timer.toc();
}

The Idea:

• Combine multiple forecasts instead of choosing one

• Combination weights may vary over time, over the distribution or both

2 Popular options for combining distributions:

• Combining across quantiles (this paper)
  • Horizontal aggregation, vincentization
• Combining across probabilities
  • Vertical aggregation
• Combining at an angle
  • Taylor & Meng (2023)

\(\ell_2\) loss:

\[\begin{equation} \ell_2(x, y) = | x -y|^2 \label{eq:elltwo} \end{equation}\]

Strictly proper for mean prediction

\(\ell_1\) loss:

\[\begin{equation} \ell_1(x, y) = | x -y| \label{eq:ellone} \end{equation}\]

Strictly proper for median predictions

Popular Aggregation Algorithms

The naive combination

\[\begin{equation} w_{t,k}^{\text{Naive}} = \frac{1}{K}\label{eq:naive_combination} \end{equation}\]

The exponentially weighted average forecaster (EWA)

\[\begin{equation} \begin{aligned} w_{t,k}^{\text{EWA}} & = \frac{e^{\eta R_{t,k}} }{\sum_{k = 1}^K e^{\eta R_{t,k}}}\\ & = \frac{e^{-\eta \ell(\widehat{X}_{t,k},Y_t)} w^{\text{EWA}}_{t-1,k} }{\sum_{k = 1}^K e^{-\eta \ell(\widehat{X}_{t,k},Y_t)} w^{\text{EWA}}_{t-1,k}} \end{aligned}\label{eq:exp_combination} \end{equation}\]

Risk

In stochastic settings, the cumulative Risk should be analyzed Wintenberger (2017):

\[\begin{align} &\underbrace{\widetilde{\mathcal{R}}_t = \sum_{i=1}^t \mathbb{E}[\ell(\widetilde{X}_{i},Y_i)|\mathcal{F}_{i-1}]}_{\text{Cumulative Risk of Forecaster}} \\ &\underbrace{\widehat{\mathcal{R}}_{t,k} = \sum_{i=1}^t \mathbb{E}[\ell(\widehat{X}_{i,k},Y_i)|\mathcal{F}_{i-1}]}_{\text{Cumulative Risk of Experts}} \label{eq_def_cumrisk} \end{align}\]

(7) expected loss of the algorithm (lower = better)

The selection problem

\[\begin{equation} \frac{1}{t}\left(\widetilde{\mathcal{R}}_t - \widehat{\mathcal{R}}_{t,\min} \right) \stackrel{t\to \infty}{\rightarrow} a \quad \text{with} \quad a \leq 0. \label{eq_opt_select} \end{equation}\] The forecaster is asymptotically not worse than the best expert.

The convex aggregation problem

\[\begin{equation} \frac{1}{t}\left(\widetilde{\mathcal{R}}_t - \widehat{\mathcal{R}}_{t,\pi} \right) \stackrel{t\to \infty}{\rightarrow} b \quad \text{with} \quad b \leq 0 . \label{eq_opt_conv} \end{equation}\] The forecaster is asymptotically not worse than the best convex combination in hindsight (oracle).

Optimal rates with respect to selection \(\eqref{eq_opt_select}\) and convex aggregation \(\eqref{eq_opt_conv}\) Wintenberger (2017):

\[\begin{align} \frac{1}{t}\left(\widetilde{\mathcal{R}}_t - \widehat{\mathcal{R}}_{t,\min} \right) & = \mathcal{O}\left(\frac{\log(K)}{t}\right)\label{eq_optp_select} \end{align}\]

\[\begin{align} \frac{1}{t}\left(\widetilde{\mathcal{R}}_t - \widehat{\mathcal{R}}_{t,\pi} \right) & = \mathcal{O}\left(\sqrt{\frac{\log(K)}{t}}\right) \label{eq_optp_conv} \end{align}\]

Algorithms can satisfy both \(\eqref{eq_optp_select}\) and \(\eqref{eq_optp_conv}\) depending on:

• The loss function
• Regularity conditions on \(Y_t\) and \(\widehat{X}_{t,k}\)
• The weighting scheme

Optimal Convergence

Requirements:

EWA satisfies optimal selection convergence \(\eqref{eq_optp_select}\) in a deterministic setting if:

Loss \(\ell\) is exp-concave

Learning-rate \(\eta\) is chosen correctly

Optimal convex aggregation convergence \(\eqref{eq_optp_conv}\) can be satisfied by applying the kernel-trick:

\[\begin{align} \ell^{\nabla}(x,y) = \ell'(\widetilde{X},y) x \end{align}\]

\(\ell'\) is the subgradient of \(\ell\) at forecast combination \(\widetilde{X}\)

Those results can be converted to any stochastic setting Wintenberger (2017)

Probabilistic Setting

An appropriate choice:

\[\begin{equation*} \text{CRPS}(F, y) = \int_{\mathbb{R}} {(F(x) - \mathbb{1}\{ x > y \})}^2 dx \label{eq:crps} \end{equation*}\]

It’s strictly proper (Gneiting & Raftery, 2007).

Using the CRPS, we can calculate time-adaptive weights \(w_{t,k}\). However, what if the experts’ performance varies in parts of the distribution?

Utilize this relation:

\[\begin{equation*} \text{CRPS}(F, y) = 2 \int_0^{1} \text{QL}_p(F^{-1}(p), y) dp.\label{eq_crps_qs} \end{equation*}\]

… to combine quantiles of the probabilistic forecasts individually using the quantile-loss QL.

Simple Example:

\[\begin{align} Y_t & \sim \mathcal{N}(0,\,1) \\ \widehat{X}_{t,1} & \sim \widehat{F}_{1} = \mathcal{N}(-1,\,1) \\ \widehat{X}_{t,2} & \sim \widehat{F}_{2} = \mathcal{N}(3,\,4) \label{eq:dgp_sim1} \end{align}\]

• True weights vary over \(p\)
• Figures show the ECDF and calculated weights using \(T=25\) realizations
• Pointwise solution creates rough estimates
• Pointwise is better than constant
• Smooth solution is better than pointwise

We extend the B-Smooth and P-Smooth procedures to the multivariate setting:

Data

• Day-Ahead electricity price forecasts from Marcjasz et al. (2023)
• Produced using probabilistic neural networks
• 24-dimensional distributional forecasts
• Distribution assumptions: JSU and Normal
• 8 experts (4 JSU, 4 Normal)
• 27th Dec. 2018 to 31st Dec. 2020 (736 days)
• We extract 99 quantiles (percentiles)

Setup

Evaluation: Exclude first 182 observations

Extensions: Penalized smoothing | Forgetting

Tuning strategies:

• Bayesian Fix
  • Sophisticated Baesian Search algorithm
• Online
  • Dynamic based on past performance
• Bayesian Online
  • First Bayesian Fix then Online

Computation Time: ~30 Minutes

Potential Downsides:

• Pointwise optimization can induce quantile crossing
  • Can be solved by sorting the predictions

Important:

• The choice of the learning rate is crucial
• The loss function has to meet certain criteria

Upsides:

• Pointwise learning outperforms the Naive solution significantly
• Online learning is much faster than batch methods
• Smoothing further improves the predictive performance
• Asymptotically not worse than the best convex combination

The profoc R Package:

• Implements all algorithms discussed above
• Is written using RcppArmadillo its fast
• Accepts vectors for most parameters
  • The best parameter combination is chosen online
• Implements
  • Forgetting, Fixed Share
  • Different loss functions + gradients

Pubications:

Berrisch, J., & Ziel, F. (2023). CRPS learning. Journal of Econometrics, 237(2), 105221.

Berrisch, J., & Ziel, F. (2024). Multivariate probabilistic CRPS learning with an application to day-ahead electricity prices. International Journal of Forecasting, 40(4), 1568-1586.

Motivation

Understanding European Emission Allowances (EUA)

Portfolio & Risk Management,

Sustainability Planing

Political decisions

EUA prices are connected to energy markets

How can the dynamics be characterized?

Several Questions arise:

• Data (Pre)processing
• Modeling Approach
• Evaluation

Data

Daily Observations: 03/15/2010 - 10/14/2022

EUA, Natural Gas, Brent Crude Oil, Coal

• normalized w.r.t. \(\text{CO}_2\) emissions
• Adjusted for inflation by Eurostat’s HICP, excluding energy

Emission-adjusted prices reflect one tonne of \(\text{CO}_2\)

Log transformation of the data to stabilize the variance

ADF Test: All series are stationary in first differences

Johansen’s likelihood ratio trace test suggests two cointegrating relationships (only in levels)

• Let \(\boldsymbol{X}_t\) be a \(K\)-dimensional vector at time \(t\)
• The forecasting target:
  • Conditional joint distribution
  • \(F_{\boldsymbol{X}_t|\mathcal{F}_{t-1}}\)
  • \(\mathcal{F}_{t}\) is the sigma field generated by all information available up to and including time \(t\)

Sklars theorem: decompose target into

• marginal distributions: \(F_{X_{k,t}|\mathcal{F}_{t-1}}\) for \(k=1,\ldots, K\), and
• copula function: \(C_{\boldsymbol{U}_{t}|\mathcal{F}_{t - 1}}\)

Let \(\boldsymbol{x}_t= (x_{1,t},\ldots, x_{K,t})^\intercal\) be the realized values

It holds that:

\[\begin{align} F_{\boldsymbol{X}_t|\mathcal{F}_{t-1}}(\boldsymbol{x}_t) = C_{\boldsymbol{U}_{t}|\mathcal{F}_{t - 1}}(\boldsymbol{u}_t) \nonumber \end{align}\]

with: \(\boldsymbol{u}_t =(u_{1,t},\ldots, u_{K,t})^\intercal\), \(u_{k,t} = F_{X_{k,t}|\mathcal{F}_{t-1}}(x_{k,t})\)

For brevity we drop the conditioning on \(\mathcal{F}_{t-1}\).

The model can be specified as follows

\[\begin{align} F(\boldsymbol{x}_t) = C \left[\mathbf{F}(\boldsymbol{x}_t; \boldsymbol{\mu}_t, \boldsymbol{ \sigma }_{t}^2, \boldsymbol{\nu}, \boldsymbol{\lambda}); \Xi_t, \Theta\right] \nonumber \end{align}\]

\(\Xi_{t}\) denotes time-varying dependence parameters \(\Theta\) denotes time-invariant dependence parameters

We take \(C\) as the \(t\)-copula

Individual marginal distributions:

\[\mathbf{F} = (F_1, \ldots, F_K)^{\intercal}\]

Generalized non-central t-distributions

• Time varying:
  • expectation \(\boldsymbol{\mu}_t = (\mu_{1,t}, \ldots, \mu_{K,t})^{\intercal}\)
  • variance: \(\boldsymbol{\sigma}_{t}^2 = (\sigma_{1,t}^2, \ldots, \sigma_{K,t}^2)^{\intercal}\)
• Time invariant
  • degrees of freedom: \(\boldsymbol{\nu} = (\nu_1, \ldots, \nu_K)^{\intercal}\)
  • noncentrality: \(\boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_K)^{\intercal}\)

VECM Model

\[\begin{align} \Delta \boldsymbol{\mu}_t = \Pi \boldsymbol{x}_{t-1} + \Gamma \Delta \boldsymbol{x}_{t-1} \nonumber \end{align}\]

where \(\Pi = \alpha \beta^{\intercal}\) is the cointegrating matrix of rank \(r\), \(0 \leq r\leq K\).

GARCH model

\[\begin{align} \sigma_{i,t}^2 = & \omega_i + \alpha^+_{i} (\epsilon_{i,t-1}^+)^2 + \alpha^-_{i} (\epsilon_{i,t-1}^-)^2 + \beta_i \sigma_{i,t-1}^2 \nonumber \end{align}\]

where \(\epsilon_{i,t-1}^+ = \max\{\epsilon_{i,t-1}, 0\}\) …

Positive vs. negative innovations (capture leverage effects).

Time-varying dependence parameters

\[\begin{align*} \Xi_{t} = & \Lambda\left(\boldsymbol{\xi}_{t}\right) \\ \xi_{ij,t} = & \eta_{0,ij} + \eta_{1,ij} \xi_{ij,t-1} + \eta_{2,ij} z_{i,t-1} z_{j,t-1}, \end{align*}\]

\(z_{i,t}\) is the \(i\)-th standardized residual from time series \(i\)

\(\Lambda(\cdot)\) is a link function:

• ensures that \(\Xi_{t}\) is a valid variance covariance matrix
• ensures that \(\Xi_{t}\) does not exceed its support space and remains semi-positive definite

Rolling-window forecasting study

• 3257 observations total
• Window size: 1000 days (~ four years)
• We sample 250 of 2227 starting points
• We draw \(2^{12}= 2048\) trajectories 30 steps ahead

Estimation

Joint (quasi) maximum likelihood estimation:

\[\begin{align*} f_{\mathbf{X}_t}(\mathbf{x}_t | \mathcal{F}_{t-1}) = c\left[\mathbf{F}(\mathbf{x}_t;\boldsymbol{\mu}_t, \boldsymbol{\sigma}_{t}^2, \boldsymbol{\nu}, \boldsymbol{\lambda});\Xi_t, \Theta\right] \cdot \\ \prod_{i=1}^K f_{X_{i,t}}(\mathbf{x}_t;\boldsymbol{\mu}_t, \boldsymbol{\sigma}_{t}^2, \boldsymbol{\nu}, \boldsymbol{\lambda}) \end{align*}\]

The copula density \(c\) can be derived analytically.

Evaluation

Our main objective is the Energy Score (ES)

\[\begin{align*} \text{ES}_t(F, \mathbf{x}_t) = \mathbb{E}_{F} \left(||\tilde{\mathbf{X}}_t - \mathbf{x}_t||_2\right) - \\ \frac{1}{2} \mathbb{E}_F \left(||\tilde{\mathbf{X}}_t - \tilde{\mathbf{X}}_t'||_2 \right) \end{align*}\]

where \(\mathbf{x}_t\) is the observed \(K\)-dimensional realization and \(\tilde{\mathbf{X}}_t\), respectively \(\tilde{\mathbf{X}}_t'\) are independent random vectors distributed according to \(F\)

For univariate cases the Energy Score becomes the Continuous Ranked Probability Score (CRPS)

Theoretical

Probabilistic Online Learning:

Aggregation
Regression

Practical

Applications

Energy Commodities
Electricity Prices
Electricity Load

Well received by the academic community:

of papers already published

108 citations since 2020 (Google Scholar)

Software

R Packages:

profoc, rcpptimer, dccpp

Python Packages:

ondil, sstudentt

Contributions to other projects:

RcppArmadillo
gamlss
NixOS/nixpkgs
OpenPrinting/foomatic-db

Awards:

Berrisch, J., Narajewski, M., & Ziel, F. (2023):

Won Western Power Distribution Competition
Won Best-Student-Presentation Award