Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
168,815 characters · 19 sections · 109 citation commands
Online Multivariate Regularized Distributional Regression for High-dimensional Probabilistic Electricity Price Forecasting
Short-term electricity markets play a key role in the integration of renewable energy sources and flexible generation in the electricity system. In Germany, the day-ahead auction is the major venue for physically delivered electricity. Trading volumes have grown with the increase of renewable generation capacity. To optimize decision-making and bidding strategies, market participants need accurate price forecasts. {\color{black} Since} electricity prices are characterized by high volatility, positive and negative spikes and skewness, research and industry have moved toward probabilistic electricity price forecasting (PEPF) to account for their stochastic nature nowotarski2018recent, dexter2024probabilistic. {\color{black} With 24 hourly prices per day, electricity prices are multivariate time series with a potentially complex dependency structure.} However, the multivariate dimension has received little attention for PEPF so far, while being of high importance for market participants in the context of the optimization of flexible assets and portfolio management lohndorf2023value, pena2024hedging, beykirch2022bidding, beykirch2024value. At the same time, the increasing availability of high-frequency data and the need for real-time decision-making in energy markets require online estimation methods for efficient model updating. This work presents an online, multivariate distributional regression model, which we apply to probabilistic day-ahead electricity price forecasting in Germany. Our work is among the first to treat the 24-dimensional hourly electricity prices as a multivariate distribution and the first to treat the problem in a strict online estimation setting, which makes the complex, high-dimensional distributional learning problem feasible on standard laptops. Our results show that modeling the dependence structure improves forecasting performance significantly compared to univariate approaches.
The literature on PEPF has evaluated a wide range of different statistical and machine learning methods, such as quantile regression, {\color{black} generalized autoregressive conditional heteroskedasticity models nowotarski2018recent, bille2023forecasting, marcjasz2023distributional,} conformal prediction methods, kath2021conformal, zaffran2022adaptive, lipiecki2024postprocessing, brusaferri2024line, distributional regression and neural network approaches ziel2021gamlss, marcjasz2023distributional, brusaferri2024nbmlss, hirsch2024online. However, these works treat each delivery hour as independent, univariate time series as in ziel2018day. {\color{black} We motivate the need for multivariate probabilistic forecasting approaches for the day-ahead electricity price by two simple plots. Figure (ref) shows a time series plot for the day-ahead electricity prices in Germany. The left inset panel shows each delivery hour as individual, daily series, emphasizing the daily co-movement. The right panel shows the cross-section, i.e. the daily shape for the first 90 days of 2017 -- here, the temporal correlation along the dimension of the delivery hours $h = 0 ,..., 23$ is visible. } Additionally, Figure (ref) shows the correlation matrix of the raw electricity prices {\color{black} (lower triangular) and} the residual correlation for a LASSO-estimated AutoRegressive model for the electricity price lago2021forecasting. We see a strong, statistically significant remaining residual cross-correlation, indicating that the resulting marginal error distributions, which are conditional on the mean, are not independent. On top of the statistical motivation, beykirch2022bidding, beykirch2024value describe the need for predicting joint distributions for the optimization of schedules and bidding curves in energy markets, further examples are provided by pena2024hedging, lohndorf2023value.
Work on multivariate probabilistic forecasting for day-ahead electricity prices is sparse in the literature and the majority of the existing work, e.g. maciejowska2024multiple, berrisch2024multivariate, han2023probabilistic, mashlakov2021assessing and agakishiev2025multivariate, does not evaluate multivariate scoring rules such as the {\color{black} Variogram score (VS), Dawid-Sebastiani-Score (DSS) or Energy Score (ES), but focus on the evaluation of the marginals of the multivariate distribution through the continuous ranked probability score (CRPS).} This reduces the problem to modeling 24 marginal distributions, taking only lagged cross-information into account. To the best knowledge of the author, only two studies truly model and evaluate the multivariate dependence structure. janke2020probabilistic approach the issue through implicit generative Copula models. grothe2023point employ the Schaake shuffle, a post-processing method for point forecasts. On the contrary, in the fields of probabilistic weather, renewable production bjerregaard2021introduction, soerensen2022recent, kolkmann2024modeling and probabilistic load forecasting gioia2022additive, browell2022covariance truly multivariate forecasting approaches have gained more attention.
The goal of distributional regression or “regression beyond the mean" kneib2023rage,klein2024distributional is modeling not only the conditional expectation, but all distribution parameters of the assumed parametric response distribution conditional on explanatory variables. The most prominent model in this regard is the original Generalized Additive Model for Location, Scale and Shape rigby2005generalized, of which numerous extensions have been developed over the last years kock2023truly,kneib2023rage,muschinski2022cholesky and distributional deep neural networks klein2021marginally, klein2023deep, rugamer2024semi. {\color{black} Due to the direct modeling of the variable's distribution}, this method is well suited for the generation of probabilistic forecasts and has been successfully applied in energy markets muniain2020probabilistic,gioia2022additive, serinaldi2011distributional,brusaferri2024nbmlss,marcjasz2023distributional. {\color{black} A drawback of fully distributional models is the computational effort. The need for efficient estimation approaches for (multivariate) distributional regression models has been recognized by umlauf2025scalable and gioia2025scalable, who propose efficient batch estimation approaches.}
For environments with large amounts of continuously incoming data, such as energy markets, online learning describes the task of updating the model given new data, without falling back on previous samples. Formally, in the strict online setting, after having seen $N$ samples of our data set, we fit a model, predict for step $N + 1$. Subsequently, we receive the realized values for $N+1$ and update our model, taking into account only the new row $N+1$. This approach allows an efficient processing of high-velocity data and results in greatly decreased computational effort. The principle is outlined in Figure (ref). Online learning for {\color{black} the absolute shrinkage and selection operator (LASSO) regularized regression} for the mean has been introduced in angelosante2009online, angelosante2010online and messner2019online. Univariate approaches suitable for probabilistic forecasting based stochastic gradient descent have been developed for specific distributions, pierrot2021adaptive, conformal prediction zaffran2022adaptive,gibbs2021adaptive,gibbs2024conformal and the generic online distributional regression in hirsch2024online. However, in the multivariate case, the literature remains sparse and focused on unconditional distributions and Copulae dasgupta2007line, zhao2022online, landgrebe2020online.
We add to the literature by presenting a generic, online, regularized, multivariate distributional regression model, allowing to model all distribution parameters conditional on explanatory variables and validate the approach in a forecasting study for the day-ahead electricity market in Germany. Our paper is one the first to tackle the issue of truly multivariate probabilistic electricity price forecasting. {\color{black} The contribution of this paper is therefore threefold:
Thereby, our contribution is valuable both from a methodological and applied perspective in research and industry alike. Reproduction code can be found on GitHub.\footnote{See: \url{https://github.com/simon-hirsch/online-mv-distreg}.} The remainder of the paper is structured as usual: The following Section (ref) introduces the multivariate, online, regularized distributional regression model. Section (ref) introduces the forecasting study, the used data and discusses multivariate forecast evaluation in detail. Section (ref) presents our results and analyzes the computational costs. {Section (ref)} concludes the paper.}
{\color{black} The following subsections introduce the building blocks for the multivariate, online distributional regression algorithm. The general structure of the algorithm and also the flow of the following section is outlined in Figure (ref). First, we like to introduce some intuition to distributional regression. Distributional regression or “regression beyond the mean" kneib2023rage aims at modeling not only the conditional expectation, but all distribution parameters of the assumed parametric response distribution conditioned on explanatory variables. Generally, we aim to model $$ \ensuremath{\bm{\mathrm{y}}} \sim \mathcal{F}(\ensuremath{\bm{\mathrm{\theta}}}_1, ..., \ensuremath{\bm{\mathrm{\theta}}}_K) \quad \text{where} \quad g_k(\ensuremath{\bm{\mathrm{\theta}}}_k) = \ensuremath{\bm{\mathrm{X}}}_k\ensuremath{\bm{\mathrm{\beta}}}_k $$ that is, we model the parameters $\ensuremath{\bm{\mathrm{\theta}}}_k$ of the distribution $\mathcal{F}$ of the response variable $\ensuremath{\bm{\mathrm{y}}}$ as regression based on the covariates in $\ensuremath{\bm{\mathrm{X}}}$. {\color{black} The link function $g(\cdot)$ ensures that the distribution parameters are in their domain.} Naturally, this distributional model is aligned with our goal of probabilistic forecasting.
{\color{black} Section (ref) gives a formal introduction to the univariate and multivariate case. To achieve an unconstrained estimation, we employ different parameterizations of the scale matrix in the multivariate case, which are briefly introduced in Section (ref). Section (ref) introduces our main contribution, the multivariate extension of online distributional regression. We further develop a path-based estimation along a sequence of the increasingly complex scale matrices, which allows for a regularized estimation and early stopping in Section (ref).} }on}.} }
We denote scalar float and integer values as lowercase letters (e.g. $a$), constants as large letters (e.g. $T$) vectors as bold, upright lower case letters (e.g. $\ensuremath{\bm{\mathrm{v}}}$) and matrices as bold upper case letters (e.g. $\ensuremath{\bm{\mathrm{A}}}$). The calligraphic $\mathcal{F}$ and $\mathcal{D}$ are reserved for (arbitrary) distributions, $\mathcal{N}$ denotes the normal distribution and $\mathcal{L}$ denotes the likelihood; other calligraphic letters (usually) denote index sets. Subscript values are usually indices in matrices, which we start with 0. Superscript indices (in square brackets) denote iterations and/or the number of samples received in the online setting.
{\color{black}
}
{\color{black} Starting from the univariate case,} distributional regression aims to model the conditional distribution parameters of the response vector $\ensuremath{\bm{\mathrm{y}}} = (y_1, ..., y_N) \in \mathbb{R}^{N \times 1}$, conditional on the covariate or explanatory data in $\ensuremath{\bm{\mathrm{X}}} \in \mathbb{R}^{N \times J}$ by adopting a parametric distribution $\ensuremath{\bm{\mathrm{y}}} \sim \mathcal{F}(\Theta)$, where $\mathbf{\Theta} = (\ensuremath{\bm{\mathrm{\theta}}}_1, ..., \ensuremath{\bm{\mathrm{\theta}}}_K)$ is a tuple of $K$ distribution parameters $\ensuremath{\bm{\mathrm{\theta}}}_k = (\theta_{k1}, ..., \theta_{kN})$. Each of the distribution parameters are linked to the covariate data through a known, twice differentiable link function $g_k(\cdot)$, leading to:
where $\ensuremath{\bm{\mathrm{\beta}}}_k$ is the coefficient vector to be estimated, relating the $J_k$ covariates in the design matrix $\ensuremath{\bm{\mathrm{X}}}_k = \ensuremath{(\ensuremath{\bm{\mathrm{x}}}_{k1}, ..., \ensuremath{\bm{\mathrm{x}}}_{kJ})^\top}$ to the distribution parameter $\ensuremath{\bm{\mathrm{\theta}}}_k$ through the link function $g_k(\cdot)$. Hence, we have:
and the probability density function $f(y_i \mid \theta_{1i}, ..., \theta_{Ki})$. The distributional regression framework therefore allows the modeling of all distribution parameters as linear regression equations of the design matrices $\ensuremath{\bm{\mathrm{X}}}_k$, which can be a subset or all of the available the covariate data $\ensuremath{\bm{\mathrm{X}}}$. Commonly additive models are employed, where $\ensuremath{\bm{\mathrm{\eta}}}_k = f_{k1}(\ensuremath{\bm{\mathrm{x}}}_{k1}) + ... + f_{kJ}(\ensuremath{\bm{\mathrm{x}}}_{kJ})$ {\color{black} and} the functions $f_{kj}(\cdot)$ can be linear terms, but also non-linear effects such as B-splines klein2024distributional,stasinopoulos2024generalized. Note that while the functions $f_{kj}(\cdot)$ might be non-linear, they can be represented by a combination of linear regression coefficients and $B$ basis functions $b(\cdot)$, i.e. $f_{kj}(\cdot) = \sum_{i=1}^B \beta_{kji}b_i(x_{kj})$. rigby2005generalized introduce iteratively reweighted least squares (IRLS), maximizing the penalized likelihood, to estimate $\ensuremath{\bm{\mathrm{\beta}}}_k$. It is important to note here that in the frequentist {\color{black} setting}, the IRLS algorithm is agnostic to the actual estimation technique stasinopoulos2024generalized. Different flavors of LASSO-type regularized estimation approaches have been introduced by groll2019lasso, muniain2020probabilistic, ziel2021gamlss, o2023variable. A regularized, incremental estimation approach using online coordinate descent has been proposed by hirsch2024online, which will form the basis for the multivariate approach proposed in this paper.
Moving to the multivariate setting, we are interested in learning the conditional distribution parameters of the $D$-dimensional response {variable $\ensuremath{\bm{\mathrm{Y}}} = (\ensuremath{\bm{\mathrm{y}}}_{1}, ..., \ensuremath{\bm{\mathrm{y}}}_{D})$}, conditional on the covariate {\color{black} matrix} $\ensuremath{\bm{\mathrm{X}}}$, by adopting a multivariate parametric distribution $\ensuremath{\bm{\mathrm{Y}}}_i \sim \mathcal{F}_D(\mathbf{\Theta}_i)$, where $\mathbf{\Theta}_i = (\ensuremath{\bm{\mathrm{\theta}}}_{i1}, ..., \ensuremath{\bm{\mathrm{\theta}}}_{iK})$ is a tuple of $K$ scalar, vector or matrix-valued distribution parameters. Each of the coordinates $m$ of the distribution parameter $\ensuremath{\bm{\mathrm{\theta}}}_k$ can again be related to its linear predictors by
This formulation is rather general. {\color{black} We note two important points here:
} The general setting introduced here includes the Gaussian multivariate distributional regression introduced by muschinski2022cholesky, the Copula-based multivariate distributional regression by kock2023truly and the MCD-based additive covariance models by gioia2022additive. The general estimation principle of repeatedly iterating through the distribution parameters until convergence translates to the multivariate case. The exact estimation algorithm will be introduced in Section (ref) and the following Section (ref) briefly discusses different options to parameterize the covariance respectively precision matrix.
{\color{black} To ensure the positive definiteness of the scale matrix, we propose three unconstrained parameterizations. To save computational costs, we parameterize the distributions in terms of the inverse covariance matrix $\ensuremath{\bm{\mathrm{\Sigma}}}^{-1} = \ensuremath{\bm{\mathrm{\Omega}}}$. This allows us to avoid matrix inversion in the evaluation of the (log-) likelihood function. We briefly review first the (modified) Cholesky decompositions (CD resp. MCD), which have been proposed in the context of multivariate distributional regression by e.g. pourahmadi2011covariance, muschinski2022cholesky, kock2023truly. Additionally, we employ a low-rank approximation of the precision matrix, which has been used by salinas2019high in the context of high-dimensional Gaussian processes and by marz2022multi and omalley2023probabilistic in the context of distributional gradient boosted random forests.}
{\color{black} The Cholesky-decomposition (CD) of the covariance matrix $\ensuremath{\bm{\mathrm{\Sigma}}}$ and the precision matrix $\ensuremath{\bm{\mathrm{\Omega}}}$ is defined as:}
muschinski2022cholesky parametrize the normal distribution in terms of $\ensuremath{\bm{\mathrm{A}}}^{-1}$ and kock2023truly choose $\ensuremath{\bm{\mathrm{A}}}$. {\color{black} Additionally the modified Cholesky-decomposition (MCD) can be used:
}align}} For the CD to yield a positive definite matrix, we require the diagonal of $\ensuremath{\bm{\mathrm{A}}}$ to be positive, which can be enforced by employing a log-link function. The lower diagonal of $\ensuremath{\bm{\mathrm{A}}}$ is unconstrained. {\color{black} The same holds for the MCD, where $\ensuremath{\bm{\mathrm{D}}}$ is a diagonal matrix with positive entries and $\ensuremath{\bm{\mathrm{L}}}$ is a unit lower triangular matrix with ones on the diagonal. Note that $\ensuremath{\bm{\mathrm{A}}}^{-1} = \ensuremath{\bm{\mathrm{D}}}^{-1/2}\ensuremath{\bm{\mathrm{L}}}$.} The low-rank approximation (LRA) is defined as
where $\ensuremath{\bm{\mathrm{A}}} = \operatorname{diag}(a_1, ..., a_D)$ and $\ensuremath{\bm{\mathrm{V}}}$ is a $D \times R$ matrix of rank $R$. The advantage of the LRA is that the dimensions of the parameters $\ensuremath{\bm{\mathrm{A}}}$ and $\ensuremath{\bm{\mathrm{V}}}$ scale linearly with the dimension $D$. However, the partial derivatives of the multivariate Gaussian and $t$-distribution with respect to the coordinates of $\ensuremath{\bm{\mathrm{A}}}$ and $\ensuremath{\bm{\mathrm{V}}}$ require inversion of the precision matrix. To ensure positive-definiteness for the LRA, we require the non-zero elements of $\ensuremath{\bm{\mathrm{A}}}$ to be positive, while $\ensuremath{\bm{\mathrm{V}}}$ is unconstrained. These requirements can easily be satisfied by choosing the log-link {\color{black} or the square root link} function for $\ensuremath{\bm{\mathrm{A}}}$.
This section introduces the online estimation algorithm for multivariate distributional regression. The algorithm is based on the combination of the iteratively reweighted least squares (IRLS) algorithm for the estimation of distributional regression models and online coordinate descent (OCD) for regularized estimation. rigby2005generalized introduce iteratively reweighted least squares for {\color{black} the estimation of GAMLSS models.} The RS algorithm consists of two nested loops, in which we cycle repeatedly through the distribution parameters and run a weighted fit of the {\color{black} working vector $\ensuremath{\bm{\mathrm{z}}}$} on the design matrix $\ensuremath{\bm{\mathrm{X}}}$ using the diagonal weight matrix $\ensuremath{\bm{\mathrm{W}}}$. The following paragraphs introduce the scoring vector and weights, the algorithm and the necessary modifications to move from a univariate case to the multivariate case. The score vector is defined as
where $\ell$ is the log-likelihood $\ell = \log(\mathcal{L})$ and $\eta = g(\ensuremath{\bm{\mathrm{\theta}}})$ is the linked predictor. The working vector for the Newton-Raphson or Fisher-Scoring algorithm is defined as
where the weights are defined as:
for Newton-Raphson and Fisher's scoring respectively. In the IRLS, estimation is done by repeatedly fitting the working vector $\ensuremath{\bm{\mathrm{z}}}$ on the design matrix $\ensuremath{\bm{\mathrm{X}}}$ using the weights $\ensuremath{\bm{\mathrm{W}}}$ until convergence. {\color{black} hirsch2024online use online LASSO estimation for linear models as building block to allow for online estimation of distributional regression models. The online LASSO is based on the Gramian matrices $$ \ensuremath{\bm{\mathrm{G}}} = \ensuremath{\ensuremath{\bm{\mathrm{X}}}^\top}\ensuremath{\bm{\mathrm{W}}}\ensuremath{\bm{\mathrm{\Gamma}}}\ensuremath{\bm{\mathrm{X}}} \quad \text{and} \quad \ensuremath{\bm{\mathrm{h}}} = \ensuremath{\ensuremath{\bm{\mathrm{X}}}^\top}\ensuremath{\bm{\mathrm{W}}}\ensuremath{\bm{\mathrm{\Gamma}}}\ensuremath{\bm{\mathrm{z}}}, $$ where $\ensuremath{\bm{\mathrm{W}}}$ and $\ensuremath{\bm{\mathrm{\Gamma}}}$ are estimation or sample and discounting weight matrices and $\ensuremath{\bm{\mathrm{G}}}$ and $\ensuremath{\bm{\mathrm{h}}}$ can be updated for the $n+1$-th observation in an online fashion as follows: $$ \ensuremath{\bm{\mathrm{G}}}^{[n+1]} = \gamma\ensuremath{\bm{\mathrm{G}}}^{[n]} + \ensuremath{\ensuremath{\bm{\mathrm{x}}}_{n+1}^\top}\ensuremath{\bm{\mathrm{x}}}_{n+1} \quad \text{and} \quad \ensuremath{\bm{\mathrm{h}}}^{[n+1]} = \gamma\ensuremath{\bm{\mathrm{h}}}^{[n]} + \ensuremath{\ensuremath{\bm{\mathrm{x}}}_{n+1}^\top}z_{n+1}, $$ where $\gamma \in (0, 1]$ is a forgetting factor, which allows to put more weight on recent data. The online LASSO estimation can then be performed by running coordinate descent on the updated $\ensuremath{\bm{\mathrm{G}}}$ and $\ensuremath{\bm{\mathrm{h}}}$, as outlined in Algorithm (ref). We denote this as $$ \ensuremath{\bm{\mathrm{\beta}}} = \textsc{OnlineLASSO}(\ensuremath{\bm{\mathrm{G}}}, \ensuremath{\bm{\mathrm{h}}}, \ensuremath{\bm{\mathrm{\lambda}}}). $$ and the details can be found in (ref) and angelosante2009online, messner2019online, hirsch2024online. } {\color{black} In the multivariate case, we extend the two nested loops of the IRLS algorithm with a third loop along all elements of the distribution parameter $k$. In the outer loop, we iterate through all distribution parameters. In the inner loop, we repeatedly run a weighted fit of the working vector $\ensuremath{\bm{\mathrm{z}}}$ on the design matrix $\ensuremath{\bm{\mathrm{X}}}$ using the weights $\ensuremath{\bm{\mathrm{W}}}$ until convergence. Algorithm (ref) gives an overview on the online estimation of multivariate distributional regression models. Note that in the inner loop, we run the weighted fit sequentially for all elements of the distribution parameter (see Lines 4 and 6). As discussed in Section (ref), we define the index sets $\mathcal{K} = \{1, ..., p\}$ for the number of parameters and $\mathcal{M}_k = \{1, .., M_k\}$ for the number of elements of each parameter.
} {\color{black} In the following paragraphs, we will discuss the details of the algorithm and the necessary modifications to move from the univariate case to the multivariate case.}
In the original GAMLSS, rigby2005generalized use Fisher's scoring {\color{black} in Equation (ref)}. Our approach generally uses Newton-Raphson scoring for the multivariate case, since the derivation of the expected value of second derivatives can be intractable, especially for more complex parameterizations of the precision matrix. Newton-Raphson scoring requires the partial derivatives of the log-likelihood function with respect to the predictors. While many previous works on distributional regression employ Newton-Raphson scoring, each derive the partial derivatives for specific combinations of distribution function and link function only o2023variable, muschinski2022cholesky. {\color{black} To facilitate the computational implementation in a mix-and-match fashion, we propose to use the first and second derivative of the log-likelihood with respect to the parameter and the first and second derivative of the link function and relate both to the necessary derivatives for Newton-Raphson scoring using the equalities
The proof is straight-forward and utilizes the chain and quotient rules and can be found in the online supplement}. We provide the necessary first and second partial derivatives of the log-likelihood with respect to the distribution parameter's coordinates, ${\partial \ell}/{\partial \theta}$ and ${\partial^2 \ell}/{\partial \theta^2}$, for all parameters for the multivariate normal and multivariate $t$-distribution given in Table (ref). The derivations can be found in {\color{black} the online supplement}.
{\color{black} For each inner iteration $i$, the update of the Gramian matrices starts at the Gramian matrices of $\ensuremath{\bm{\mathrm{G}}}_{km}^{[n]}$ and $\ensuremath{\bm{\mathrm{h}}}_{km}^{[n]}$ and the new information enters through the update {\color{black} of the weights $\ensuremath{\bm{\mathrm{W}}}$ and the working vector $\ensuremath{\bm{\mathrm{z}}}$}. However, the weights are also updated iteratively along each inner and outer iteration $i$ and $r$ due to the Newton-Raphson step toward the optimal coefficients. The weights can only be updated for the current update step $n+1$, while previous weights remain fixed. In a pure batch case, all weights are updated within each Newton-Raphson step. This introduces an approximation error for the online case, which can be controlled by the forget parameter $\gamma$ as shown in hirsch2024online. }
{\color{black} For each element of the distribution parameter, we estimate a regularization path $\ensuremath{\bm{\mathrm{\beta}}}_{km}$ for a decreasing $\ensuremath{\bm{\mathrm{\lambda}}}_{km}$.} This raises the issue of model selection, i.e. the selection of the optimal regularization $\lambda_{mk}^\text{opt}$. We propose to use information criteria (IC), as it is well-aligned to the likelihood-based framework of distributional regression. Define a generalized IC
where $\ell$ is the log-likelihood under the model, $K$ is the number of parameters in the model and $N$ the number of seen observations. We can recover Akaikes Information Criterion (AIC), the Bayesian Information Criterion (BIC), and the Hannan-Quinn Information Criterion (HQC) by setting $\nu_0, \nu_1, \nu_2$ accordingly. The optimal regularization parameter is then selected as $\lambda_{mk}^\text{opt} = \operatorname{argmin}_{\lambda} \operatorname{IC}$. Since the evaluation of the likelihood can be costly for high-dimensional data, we propose to employ the first derivative of the log-likelihood, i.e. calculate
where the superscript ${[\lambda_i]}$ denotes the model with the regularization parameter $\lambda_i$. The approximation is valid for small changes in the regularization parameter and avoids the costly re-evaluation of the likelihood.
The algorithm runs iteratively along all coordinates of the distribution parameter {\color{black} (line 6). Since the partial derivatives are not information orthogonal, }the coordinates of the distribution parameters might impact each other, e.g. in the matrix multiplication of the CD-based scale matrix {\color{black} (see also the definition of the derivatives in Table (ref) and the online supplement)}. To stabilize the estimation, we propose to update the values in the very first iteration $i$ by a “dampened" version, i.e. taking
Hence, the predictions from the first iteration will be the average of the first fitted values and the initialization. This feature is mainly important for the scale matrix, whose coordinates are usually not orthogonal and less so for the location and (scalar) tail parameters. {\color{black} For the same reason, the options for parallelization remain limited unfortunately.} For the multivariate normal and $t$-distribution used in this paper, only the estimation of the location parameter can be parallelized, as well as the estimation of the coordinates of the LRA matrix {$\ensuremath{\bm{\mathrm{A}}} = \operatorname{diag}(a_1, ..., a_D)$} for the normal distribution. {\color{black} This would allow to parallelize lines 6-10.} For the $t$-distribution, the estimation of $\ensuremath{\bm{\mathrm{A}}}$ can only be parallelized for sufficiently high degrees of freedom. {\color{black} A potential strategy for parallelization with non-orthogonal parameters could be using a step size smaller than one, i.e. using a convex combination of the previous and the newly estimated coefficients (see Equation (ref)), which, in combination with a Cole-Green-algorithm green1984iteratively, cole1992smoothing, could allow to parallelize the estimation of all elements of all distribution parameters (lines 5-10). However, this strategy introduces an additional hyperparameter and we incur further questions with respect to individual or joint regularization and model selection, as we have to select the optimal regularization parameter for each coordinate of each distribution parameter. Hence, we leave the question of parallelization for future research. }
{\color{black} Using LASSO-type regularization, the algorithm can handle high-dimensional covariate spaces for each coordinate of the distribution parameter. However, for high-dimensional response variables, i.e. large $D$, the number of parameters in the scale matrix grows quadratically in $D$ {\color{black} for most unconstrained parameterizations}. To alleviate this issue and allow for parsimonious modeling for large $D$, we propose a path-based estimation approach that starts with a simple (highly regularized) structure of the scale matrix and gradually increases its complexity.
{\color{black} This approach is inspired by path-based estimation methods in high-dimensional statistics friedman2008sparse}.}rse}}.} We exploit that in many cases, some structure can be imposed on the scale matrix, i.e. in spatial or temporal data, which has a clear dependence pattern along the diagonal. In these cases, the scale matrix can be regularized by systematically setting off-diagonal elements to zero gabriel1962ante,zimmerman1997structured,zimmerman1998computational. While both, the CD-based and the LRA-based scale matrix parameterization lend themselves to this type of regularization, the approach is mainly popular with the Cholesky-based parameterizations due to the relationship between the elements of the CD and the temporal correlation for longitudinal data under the name {\color{black} antedependence-$r$ (AD-$r$)} regularization. However, such regularization is commonly applied a-priori and not in a data-driven fashion, see e.g. muschinski2022cholesky,zimmerman1997structured. On the other side, in coordinate descent estimation of regularized problems such as (graphical) LASSO, path-based estimation starting from a strongly regularized solution toward an (almost) not regularized solution has proven itself as an efficient solution approach. In this section, we aim to combine these two principles by introducing path-based estimation for the regularized scale matrix.
{\color{black} On a high level, our algorithm starts with an “independence-parameterization" of the scale matrix and subsequently adds more non-zero elements and thus complexity to the parameterization of the scale matrix. Figure (ref) illustrates how the path-based estimation uses increasingly complex specifications for the scale matrix $\ensuremath{\bm{\mathrm{\Sigma}}}$ respectively $\ensuremath{\bm{\mathrm{\Omega}}}$. Formally, for the regularization parameter $\alpha$,
We can use warm-starting for all previously fitted elements of the scale matrix, however, due to the non-orthogonality of the elements, we need to re-estimate all elements of the scale matrix in each iteration. We discuss some details in the following:}
Currently, the Algorithm will add only full off-diagonals {\color{black} (for Cholesky-based approaches)} respectively columns (LRA). The implementation however could also work for block-wise schemes cai2012adaptive or user-defined regularization patterns. {\color{black} The development of smart selection schemes for choosing the next elements to include in the scale matrix} would be beneficial for the speed of the algorithm. {\color{black} We provide an analysis of the in-sample selection of $\alpha$ and the out-of-sample performance for various $\alpha$ in the forecasting study in (ref) and leave the development of more advanced selection schemes for future research.}
For electricity produced on day $t$ and in hour $h$, the short-term electricity market in Germany is split in three major parts: The daily day-ahead auction on $t-1$ at 12:00 hours for 24 hourly delivery periods $h \in \{0, ..., 23\}$, the afternoon auction with quarter-hourly delivery periods on $t-1$, at 15:00 hours and the continuous intraday market. The daily procedure for the day-ahead auction, which is the focus of this paper, is shown in Figure (ref). The market is organized by EPEX SPOT and Nordpool in the joint single day-ahead coupling (SDAC) as a pay-as-cleared auction, resembling the merit-order model for the electricity market bille2023forecasting,hirsch2024online,viehmann2017state.
At each day $t$, before the day-ahead auction at 12:00, we aim to generate forecasts for day $t+1$. Prior to forecasting, we update our models by taking into account the realized prices and forecasts for delivery day $t$. {\color{black} Figure (ref) contrasts online learning with the expanding and rolling window batch learning, two popular schemes for forecasting studies. Note that in the online learning scheme, we only use the new observation for updating the model, while in the batch learning schemes, we re-use all or a fixed window of previous observations. The rolling window scheme is popular in the EPF community lago2021forecasting,nowotarski2018recent while online learning schemes are a rather recent development in EPF berrisch2024multivariate, zaffran2022adaptive,hirsch2024online.}
{\color{black} We use the same data set as in lipiecki2024postprocessing, which consists of electricity prices for the German day-ahead market from 2015-01-01 to 2024-01-01 exclusive. In line with previous works, we use the data until 2018-12-26 as initial training set, leaving 1831 observations and therefore more than 4 years, as it is best practice lago2021forecasting, for out-of-sample testing. This dataset has been employed in various studies and to enable comparability, we further split the dataset in two sub-samples, that is 2018-12-27 to 2020-12-31 (736 days) and 2021-01-01 to 2023-12-31 (1095 days).} Additionally, the data set contains day-ahead renewable production forecasts, load forecasts and prices for fundamental commodities. All features are briefly described in Table (ref). {\color{black} We use an incremental mean-variance scaling to {\color{black} pre}-process the electricity prices. We denote the electricity price for day $t$ and hour $h \in 0, ..., 23$ as $P_{th}$. We have $y_{th} = (P_{th} - \tilde{\mu}_{th}) / \tilde{\sigma}_{th}$, where $\tilde{\mu}_{th} = 1/t \sum_i^t P_{ih}$ and $\tilde{\sigma}_{th} = \sqrt{1/t \sum_i^t (P_{ih} - \tilde{\mu}_{ih})^2}$ are the mean and standard deviation up to observation $t,h$ and therefore have $\ensuremath{\bm{\mathrm{Y}}} = (\ensuremath{\bm{\mathrm{y}}}_0, ..., \ensuremath{\bm{\mathrm{y}}}_{23})$ as the 24-dimensional ($D=H=24$) response matrix. We find that this stabilizes the estimation of covariance matrices and the normalization can be re-applied after the estimation, i.e. the covariance of $\ensuremath{\bm{\mathrm{P}}}_t$ corresponds to $\operatorname{diag}(\tilde{\ensuremath{\bm{\mathrm{\sigma}}}}_t) \ensuremath{\bm{\mathrm{\Sigma}}} \operatorname{diag}(\tilde{\ensuremath{\bm{\mathrm{\sigma}}}}_t)$, where $\ensuremath{\bm{\mathrm{\Sigma}}}$ is estimated based on $\ensuremath{\bm{\mathrm{y}}}_t$. Incremental updates of the mean and variance are straight-forward in an online learning setting using Welford's method welford1962note.} Note that, in this application study, we have $T$ corresponding to $N$ and $H$ corresponding to $D$ in the general notation.
We propose modeling the multivariate distribution of the day-ahead electricity prices in increasing complexity. {\color{black} If possible}, models are updated online, i.e. using only the new data for each day. We differentiate between an adaptive estimation, which is updating a single, unconditional distributional parameter and the full conditional estimation linking the distribution parameter to explanatory variables. {\color{black} Table (ref)} shows the increasing model complexity. {\color{black} For all three complexity levels, we include a reference model that assumes independence to showcase the value-added of modeling the dependence structure.} We start with the established LEAR models and naively estimate the unconditional residual distribution (denoted as LEAR-N(0, $\sigma$) and LEAR-N(0, $\Sigma$)). {\color{black} Additionally, we employ two univariate benchmarks. First, we use a conformal prediction (CP) approach, which is based on the LEAR model and uses a combination of the split-conformal prediction (SCP, {\color{black} during initial training}) and adaptive conformal prediction gibbs2021adaptive,gibbs2024conformal. We use absolute residual scores and a calibration set of {\color{black} 200 days per delivery hour}. {\color{black} This model is denoted as LEAR-CP.} Conformal prediction is generally a univariate post-processing method, since it is centered on issuing prediction intervals and therefore not suitable for the generation of ensemble forecasts. Multivariate conformal prediction approaches exist, but are in their infancy and focussed on a notion of multivariate quantiles.\footnote{As a simple workaround, we have tried to combine CP with a copula-based approach using a {\color{black} probability integral transform (PIT).} However, the conformal predictive density approximated using 199 quantiles was not sufficient to generate useful ensembles in the inverse transformation, i.e. the step from the simulated copula on $\mathcal{U}(0, 1)$ to the original domain failed.} Second, we use a GARCH(1,1) model with a Gaussian distributional assumption (denoted as LEAR-GARCH). We estimate the GARCH model on the residuals of the LEAR model using the arch package in Python kevin_sheppard_2024_14035889. {\color{black} This is the only model that we cannot run online in the study.}}
We increase complexity by moving to a full distributional model for the marginals and adding an adaptive estimation of the dependence structure using the Gaussian copula {\color{black} (denoted as ODR-Copula), where the Copula is either the independence Copula (ODR-IND), or a (sparse) Gaussian copula (ODR-GC and ODR-SPGC).} Lastly, we estimate the full multivariate distribution in a conditional way using the proposed multivariate online distributional regression approach {\color{black} (denoted as MODR-Method-Parameterization-Dependence)}. We describe the full model in the following and note that we additionally describe the hyperparameters in (ref).
We model the mean/location for all regression models, {\color{black} also for the LEAR models,} by
We model the scale parameters for univariate distributional models, as well as the elements of the Cholesky-factor $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{(\ensuremath{\bm{\mathrm{A}}}^{-1})^\top}(\ensuremath{\bm{\mathrm{A}}}^{-1})$, and the elements of the diagonal matrix $\ensuremath{\bm{\mathrm{A}}}$ in the LRA-based scale matrices by
where $\operatorname{SignedSquare}(a) = \operatorname{sign}(a)\sqrt{|a|}$ is the signed square root and $\ensuremath{\bm{\mathrm{\Sigma}}}^{[t-1:t-7]}$ is the rolling empirical covariance matrix of $\ensuremath{\bm{\mathrm{y}}}_t$ for the last 7 days. {\color{black} Note that we use the inverse of the signed square, since we are working on the precision matrix. For the univariate models, we replace this accordingly with the empirical rolling standard deviation.} For the LRA-based parameterization, we choose $r=2$ and model the elements of $\ensuremath{\bm{\mathrm{V}}}$ as
that is, the first rank takes most of the fundamental variables, while the second rank contains the weekday binary variables. The degrees of freedom are modeled as
The univariate models are therefore a slight simplification compared to the models used in hirsch2024online, thereby the multivariate distributional regression models and the Copula-based approaches are better comparable. Lastly, let us remark on the online tracking of the Gaussian copula. The probability density function (PDF) for the Gaussian copula is given by:
where $\ensuremath{\bm{\mathrm{u}}}$ are the pseudo-observations on the $\mathcal{U}(0, 1)$ space, $\ensuremath{\bm{\mathrm{n}}} = \Phi^{-1}(\ensuremath{\bm{\mathrm{u}}})$, $\Phi$ is the CDF of the standard normal distribution and $\tilde{\ensuremath{\bm{\mathrm{\Sigma}}}}$ is the covariance matrix $\ensuremath{\bm{\mathrm{\Sigma}}}$ scaled to the correlation matrix, $\ensuremath{\bm{\mathrm{I}}}$ is the identity matrix and $p(x_d \mid \ensuremath{\bm{\mathrm{\theta}}}_d)$ is the likelihood of the observation $y_d$ under the (conditional) marginal distribution kock2023truly,arbenz2013bayesian. We fit the Copula model by transforming the in-sample data to the uniform space $\ensuremath{\bm{\mathrm{u}}}$ by the probability integral transform (PIT) and subsequently transforming to the $\mathcal{N}(0, 1)$ space $\ensuremath{\bm{\mathrm{n}}}$, on which we can fit the dependence structure. We update the scale matrix of the Gaussian copula by taking
where $\ensuremath{\bm{\mathrm{n}}}$ are the PIT-transformed in-sample values and the superscript $[t]$ denotes the observations available in the online learning dasgupta2007line. Samples are drawn from the Gaussian copula in the usual manner. We use the same principle to track the residual covariance structure for the LEAR models under the normality assumption. We employ a second model, where we sparsify the estimated dependence matrix of the Gaussian copula {\color{black} using} the graphical LASSO friedman2008sparse.
{\color{black} Forecast evaluation should follow the well-known principle of sharpness subject to calibration gneiting2007probabilistic}. {\color{black} We check the calibration of the forecasts by calculating joint prediction bands (JPB) from the simulations. {\color{black} JPB aim to cover the true price vector with probability $1-\alpha$ and therefore generalize marginal prediction intervals to the multivariate case staszewskabystrova2011bootstrap, luetkepohl2015comparison. They have been used in energy market forecasting by serafin2022trading, chen2025probabilistic before and proven effective in supporting trading decision-making.} As it is standard in probabilistic forecasting, we evaluate both marginal and multivariate quality of the forecasts. Therefore, we} employ the root mean square error (RMSE), mean absolute error (MAE) and the continuous ranked probability score (CRPS), {\color{black} which focus on point prediction and the marginal distribution.} We employ four well-established multivariate probabilistic scoring rules: The Energy Score (ES), the Dawid-Sebastiani Score (DSS), the Variogram Score (VS) and the Log-Score (LS). {\color{black} Let us briefly review some recent results with respect to the four multivariate scores:
} Additionally, we test for statistically significant score differences using the well-established Diebold-Mariano test. The following paragraphs introduce {\color{black} JPBs and the scoring rules} and are largely based on gneiting2007probabilistic, gneiting2007strictly, nowotarski2018recent, marcotte2023regions, ziel2019multivariate as well as the references mentioned for the individual scores. Denote the true price vector as $\ensuremath{\bm{\mathrm{Y}}} = (\ensuremath{\bm{\mathrm{y}}}_0, ..., \ensuremath{\bm{\mathrm{y}}}_H)$ of shape $T \times H$ and the ensemble forecast as $\ensuremath{\bm{\mathrm{F}}}$ of shape $T \times H \times M$ of $M = 2500$ samples.
{\color{black} A joint prediction band for the $1-\alpha$ coverage is defined by the lower bound $\ensuremath{\bm{\mathrm{l}}}_t = (l_{t,0}, ..., l_{t,23})$ and upper bound $\ensuremath{\bm{\mathrm{u}}}_t = (u_{t,0}, ..., u_{t,23})$, such that
that is, the true price trajectory $\ensuremath{\bm{\mathrm{y}}}_t$ is fully covered by the JPB with probability $1-\alpha$. This is in difference to marginal prediction intervals based on predicted quantiles, which consider element-wise coverage. There are multiple algorithms to construct such bands and following serafin2022trading, chen2025probabilistic, we use the neighboring paths method described by staszewskabystrova2011bootstrap, luetkepohl2015comparison. The method is based on iteratively removing paths from the ensemble forecast $\ensuremath{\bm{\mathrm{F}}}$, such that the envelope of the remaining paths covers the true price trajectory $\ensuremath{\bm{\mathrm{y}}}_t$ with probability $1-\alpha$. It should be noted that JPB are in general not unique, i.e. there are multiple sets of lower and upper bounds that cover the true trajectory with the desired probability. We compare the mis-coverage of the {\color{black} joint prediction band}, which is defined as
and gives a measure of the multivariate calibration of the joint prediction. Additionally, we evaluate the mean width of the {\color{black} JPB}, which is defined as
which gives a measure of the efficiency of the prediction bands and can be interpreted as the area of the prediction band divided by the dimensions.} The RMSE is defined as
where $\hat{\ensuremath{\bm{\mathrm{\mu}}}}_t = \frac{1}{M}\sum_{m=0}^M \ensuremath{\bm{\mathrm{F}}}_{t,h,m}$ is the mean prediction vector. The MAE is defined as
where $\operatorname{median}\left(\ensuremath{\bm{\mathrm{F_t}}}\right)$ denotes the median trajectory for each day $t$. The CRPS is estimated from the forecast ensemble by using the probability-weighted moment estimator of zamo2018estimation:
The CRPS is a strictly proper scoring rule for the marginal distribution. Note that many works on energy price forecasting report the average pinball loss (APS) as CRPS, which needs to be rescaled $\operatorname{CRPS} = 2 \cdot\operatorname{APS}$ to be comparable. The energy score gneiting2007strictly is defined as
The energy score is a strictly proper scoring rule. We aggregate the ES by taking the average: $\operatorname{ES} = \frac{1}{T} \sum_{t=0}^T \operatorname{ES}_t$. The Log-Score (LS) is defined as
where $\mathcal{L}$ is the underlying likelihood or probability density function of the distribution $\mathcal{D}$ and $\hat{\ensuremath{\bm{\mathrm{\theta}}}}^\mathcal{D}_t$ is the estimated parameter vector. Again, we aggregate the LS by simple averaging over all points in the test set $\operatorname{LS} = \frac{1}{T} \sum_{t=0}^T \operatorname{LS}_t$. It is a strictly proper scoring rule. The Dawid-Sebastiani-Score dawid1999coherent is defined as
where $\hat{\Sigma}_F$ denotes the empirical covariance of the forecast ensemble $\ensuremath{\bm{\mathrm{F}}}$ and $\hat{\ensuremath{\bm{\mathrm{\mu}}}}$ denotes the mean ensemble as above. We aggregate the $\operatorname{DSS} = \frac{1}{T} \sum_{t=0}^T \operatorname{DSS}_t$ by simple averaging. The DSS is a proper scoring rule for the first and second moment and strictly proper for the Gaussian predictive distribution, since it is a linear transformation of Gaussian log-likelihood. The Variogram Score scheuerer2015variogram is defined as
and is a proper scoring rule. We aggregate the VS by taking the average and normalize the score by dividing by $H^2$, i.e. $\operatorname{VS} = \frac{1}{TH^2}\sum_{t=0}^T{\operatorname{VS}_t}$ to make the scales of the score comparable. The scoring rules used are implemented in the Python package scoringrules zanetta2024scoringrules.
Conclusions on the performance of forecasting models cannot be {\color{black} derived by} looking at aggregate scores alone, but need to be drawn by evaluating whether the differential between the loss series of two models is statistically significantly from zero diebold2002comparing, diebold2015comparing. For the DM-test, we evaluate the differential of two score series $\Delta\ensuremath{\bm{\mathrm{s}}}^{\mathcal{A},\mathcal{B}} = \ensuremath{\bm{\mathrm{s}}}^\mathcal{A} - \ensuremath{\bm{\mathrm{s}}}^\mathcal{B}$, where $\ensuremath{\bm{\mathrm{s}}}^\mathcal{A} = (s_0^\mathcal{A}, ..., s_T^\mathcal{A})$ are the scores for each scoring rule at $t$ for model $\mathcal{A}$ respectively $\mathcal{B}$. We provide two one-sided and hence complementary tests. {\color{black} To ensure validity of the DM-test, we check stationarity of the differential series $\Delta\ensuremath{\bm{\mathrm{s}}}^{\mathcal{A},\mathcal{B}}$ by the augmented Dickey-Fuller test dickey1979distribution, cheung1995lag.}
{\color{black} We present the results of our forecasting study in three parts. First, we analyze the predicted time-varying dependence structure and discuss the economic interpretation. Second, we present the forecasting accuracy of all models using the scoring rules presented above. Lastly, we discuss the computational costs of the different models and compare the online and batch estimation. }
{\color{black}
} {\color{black} Figure (ref) plots the evolution of the predicted volatility (standard deviation) and correlation matrix for the {MODR-OLS-MCD-$\Sigma$} model over time. The lower three panels show the 1st, 2nd and 3rd off-diagonal of the correlation matrix over time. We see that the correlation is the lowest around the hours 5-7, which correspond to start of the morning ramp. The correlation between hours $h$ and $h+i$ decays stronger for days with lower levels of volatility and vice versa. There are some weekly patterns visible: The dependence between the hours 5-7 and 16-19 tapers off more strongly on working days. This is likely driven by stronger shapes in prices and residual demand. Overall, we see that the dependence structure is not constant over time, which underlines the importance of including a time-varying dependence structure in multivariate electricity price forecasting. }
\FloatBarrier {\color{black}
}
{\color{black} {\color{black} We first examine the calibration of the joint prediction bands in Figure (ref). The figure shows miscoverage and mean JPB width for the $50\%, 55\%, ..., 95\%$ nominal coverage joint prediction bands. Table (ref) gives the results for the scoring rules for each model. Figure (ref) gives daily skills scores of the multivariate distributional regression models over the benchmarks {LEAR-N(0,$\Sigma$)}. As a mental guidance for the increasing complexity of the online regression models, remember Section (ref) and Table (ref): The LEAR-based models have an adaptive, but unconditional estimation for the probabilistic forecast. The univariate online distributional models (ODR) employ conditional estimation for all marginal distributional parameters and an adaptive, but unconditional estimation for the Copula, while the multivariate online distributional regression (MODR) models yield an online estimated conditional multivariate distribution. We summarize the main results before we discuss the results in more detail. }
}e}}
The $p$-values of the DM-test {\color{black} (see Figure (ref) in (ref))} largely confirm the statistical significance of the aforementioned results. We note the strong performance of the Copula-based models for the ES and the statistically significant superior performance of the multivariate distributional regression for the DSS, LS, and the VS. {\color{black} We note the disconnect between the results of the ES on the one hand and the VS, DSS and LS on the other hand, which is in line with the findings of the asymmetry of scoring rules by buchweitz2025asymmetric and the different behaviors observed by pinson2013discrimination, marcotte2023regions and alexander2024evaluating.}
{\color{black} Looking at Figure (ref), we see that skill scores are also well correlated with the current market regime. Skill scores are defined as $$ \operatorname{SS}_\text{model} = 1 - \frac{\operatorname{LS}_\text{model}}{\operatorname{LS}_\text{baseline}}, $$ where we use the LEAR-N(0,$\Sigma$) model as baseline. {\color{black} The MODR-LASSO models achieve up to 10% improvement to the baseline LEAR-N(0,$\Sigma$) and even provide small gains during the high-volatility period of 2022. Especially during this period, the MODR models outperform the univariate ODR-models.} Additionally, we see that the models assuming independence have an increasingly better performance in recent years. This can be interpreted as a sign that the dependence structure has both weakened and changed with the {\color{black} turbulent} market conditions. }
{\color{black} A common theme for distributional regression models is the trade-off between accuracy in the point predictions and the distributional error metrics. marcjasz2023distributional and hirsch2024online note this observation in the univariate case, trading RMSE for CRPS. This behavior is rooted in the fact that observations with predicted high variance are weighted down in the estimation of the mean parameters, which leads to slightly worse point predictions, but better distributional forecasts. At the same time, these observations with high variance are especially costly in terms of the RMSE. We see that the behavior extends in the multivariate case as well, trading off marginal accuracy for somewhat better modeling of the dependence structure in otherwise equal models. So far, this issue has not been discussed for multivariate distributional regression and yet deserves further attention in future research gioia2022additive,muschinski2022cholesky. }
{\color{black} For the multivariate online distributional regression models, we see that the models using the Cholesky-based parameterizations {\color{black} (CD and MCD)} provide better performance. This is likely due to the fact that these parameterizations are closer to the natural, time-based structure of the (conditional) dependence than the LRA and therefore the path-based regularization allows to select a parsimonious model, which still reflects the shape of the dependence structure. {\color{black} We further discuss the scale matrix regularization in (ref). For regularization with respect to the feature space,} using LASSO yields slight improvements in the scoring rules . {\color{black} The rather marginal gain for aggregate scores might be explained by the fact that all fundamental variables are known to be relevant for electricity price formation, and hence the regularization does not remove many variables, but helps to stabilize the estimates in the online learning setting.} }
{\color{black} Additionally, we analyze the performance of the models for extreme price spikes and for days with large spreads during the day. These days are interesting from a risk-management respectively from a battery optimization perspective. We define a price spike as a day with a minimum or maximum price exceeding the 5% resp. 95%-quantile of all minimum resp. maximum prices. A large spread event is defined as the min-to-max spread exceeding the 90%-quantile of all min-to-max spreads. Formally:
The results for the CRPS, LS and the miscoverage are given in Table (ref). For the CRPS, we see that the ordering of the models is similar to the overall results, and remains largely unchanged between spike and no-spike events. For the miscoverage, we see that the multivariate distributional regression models yield the best results for both regimes for large spread events, and the best results for the price spike events, while the Copula-based models yield the best results for the no-spike events. Overall, these results highlight the robustness of the distributional regression models also for extreme events. {\color{black} The results for the ES, DSS and VS can be found in Table (ref) and mirror roughly the results for the CRPS and LS. Figures (ref) and (ref) give the miscoverage of the $50\%, 55\%, ..., 95\%$-joint prediction bands for the spike and large spread events. The disaggregated analysis reveals undercoverage during spikes for almost all models for smaller nominal coverage levels, where the LEAR and the ODR-SPGC deliver robust performance, while for larger nominal coverage levels, the undercoverage is less severe. }}
{\color{black} Before concluding the results, we want to emphasize some limitations of this study.
Overall, our results highlight that neglecting the dependence structure by relying solely on marginal, univariate models yields subpar probabilistic forecasting performance. We note that for the truly multivariate approaches, using both, copula-based combinations of univariate models and the fully multivariate distributional regression yield statistically significant performance improvements. In this regard, this paper is the first to carry out a comprehensive, multivariate probabilistic forecasting study on the day-ahead market, including also the challenging years of the COVID-19 pandemic and the energy crisis following the Russian invasion of Ukraine in the test set. }
\FloatBarrier {\color{black}
}
Table (ref) gives computation times for all experiments. The initial fit for the multivariate distributional regression model takes a few minutes, the update algorithm can be executed in seconds. {\color{black} Using online estimation methods, the experiments can be run in a few hours on a standard laptop.} An estimate for the benefit of online vs. repeated batch fitting can be achieved by multiplying the initial fit duration with 1831 days of out-of-sample and comparing this to the total time of the online study: $$ \text{Speedup} = \frac{\text{Initial Fit} \times T}{\text{Total Time}}. $$ By this (albeit simple) measure, the online learning improves computation by a factor {\color{black} of 80 to 600}. These estimates are in line with benefits reported in hirsch2024online for the univariate online distributional regression case in {\color{black} direct} comparison {\color{black} between online estimation and repeated batch estimation}.
{\color{black} To further analyze the trade-off between computation time and accuracy, we take the OLS-estimated distributional regression model and estimate the model every 1, 7, 14, 30, 60, 180, 365 days online (using a mini-batch online update) and, vice versa, re-estimate the full model in a repeated batch fitting every 7, 14, ..., 365 days, using the first subsample of the test data.\footnote{Note that already estimating the batch model every 7 days takes more than 10 hours and hence we did not run the experiment for daily estimation in the batch case. {\color{black} We expect daily refitting scheme to deliver further improvements in the predictive accuracy, but at about seven times the computational time; the linear relationship between the computation time and update frequency is clearly visible in Figure (ref).}} We compare the results for online learning, rolling and expanding window batch estimation in Figure (ref). {\color{black} For online and batch estimation, we see that higher update frequency improves forecasting accuracy, but also increases computation times.} Crucially, the “efficient frontiers“ of both approaches never intersect. We see that for low computation time budgets, online learning approaches give strictly better results than repeated batch fitting approaches and only for large computation time budgets, expanding window batch estimation takes the lead. A second important observation is that the online learning algorithm is almost constant in the time per update, since we always process a fixed amount of data points. On the other hand the number of data points in the expanding window batch estimation increases, which increases the computation time per fit. Interestingly, the rolling window estimation scheme, which is often favored in the EPF literature, delivers worse predictive accuracy than an expanding window scheme. These results further confirm the estimated speed-up in Table (ref). } {\color{black} Summarizing the results on the computational effort, our online update algorithm makes the approach practically viable for researchers and data scientists without access to specialized high-performance computation centers. {\color{black} Since one good model is usually preceeded by fitting a large amount of bad models, taking} into account the trade-off between accuracy and computational effort is of practical importance in many industrial settings, where analysts work under time pressure and data arrives at high velocity. Furthermore, time saved in raw computation can be used for better data exploration and feature engineering. }
\FloatBarrier
Distributional learning algorithms such as GAMLSS and deep distributional networks have been used successfully for probabilistic electricity price forecasting muniain2020probabilistic,hirsch2024online,marcjasz2023distributional. However, even for univariate distributions, these models are computationally expensive. At the same time, the literature on probabilistic electricity price forecasting has largely focused on modeling the hourly marginal distributions only, leaving the dependence structure neglected. Against this background, we develop an online estimation algorithm for multivariate distributional regression models, {\color{black} making the use of these algorithms feasible on a standard laptop, even for high-dimensional problems such as the 24-dimensional distribution of electricity prices}. We benchmark our implementation in a forecasting study for the German day-ahead electricity market and thereby provide the first study exclusively focused on online learning for multivariate PEPF.
{\color{black} Our results show that modeling the dependence structure in the day-ahead market improves probabilistic forecasting performance significantly. First, we see that calibration for prediction bands of the 24-dimensional price path is improved significantly by modeling the dependence structure and moving toward proper joint prediction bands staszewskabystrova2011bootstrap instead of marginal, quantile-based prediction bands. The online, multivariate distributional regression models deliver strong predictive accuracy across a range of multivariate scoring rules. Additionally, we like to highlight the importance of regularization to avoid overfitting in a high-dimensional setting and conduct two experiments to validate the regularization of the scale matrix. Distributional regression models are interpretable and therefore allow us to discuss the economic interpretation of the time-varying dependence structure.
We analyze the trade-off between computation time and forecasting accuracy. Building an efficient frontier between accuracy and computation time we show that online learning, compared to repeated batch fitting, yields better results for a given computation budget -- with speed-ups of 2-3 orders of magnitude.} {\color{black} Our algorithm is implemented in a fairly generic manner, allowing e.g. for different distributional assumptions and keeping a familiar, sklearn-like API to facilitate the usage by other researchers and data scientists pedregosa2011scikit and contributed the implementation to the ondil package hirsch2024online.\footnote{See: \url{https://github.com/simon-hirsch/ondil}}}
Our research opens multiple avenues for future work. First, further research on the driving forces of the dependence structure in the German electricity market is necessary to improve the forecasting performance and guide decision-making processes in electricity trading. Modeling the dependence structure in electricity markets is a rather open field and has implications beyond forecasting, concerning also risk and portfolio management and asset optimization pena2024hedging,lohndorf2023value, beykirch2022bidding, beykirch2024value. From an algorithmic perspective, we note that while our algorithm is already quite fast, {\color{black} further improvements could be achieved by distributed and parallel computing. We outline a potential strategy for future research in Section (ref). A related issue }is model selection -- while the regularized online estimation is fast, the models are still quite complex and can be prone to overfitting. Lastly, due to the generic nature of our implementation, the usage for other high-dimensional forecasting problems such as probabilistic wind, solar, and load forecasting can be explored.
\FloatBarrier
Simon Hirsch is employed as an industrial PhD student by Statkraft Trading GmbH and gratefully acknowledges the support and funding received. This work contains the author’s opinions and does not necessarily reflect Statkraft’s position. Simon Hirsch is grateful for helpful discussions with Florian Ziel, Daniel Gruhlke, Christoph Hanck and the participants of the IWSM 2025 in Limerick.
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The data used in this paper has been provided by lipiecki2024postprocessing and cannot be shared further by the author. Therefore, for any requests with respect to the data, please contact lipiecki2024postprocessing directly. A shorter data set of the same explanatory variables is available at the GitHub repository of marcjasz2023distributional.\footnote{See: \url{https://github.com/gmarcjasz/distributionalnn}.} For the differences between both data sets please consult lipiecki2024postprocessing.
During the preparation of this work the author used generative AI tools such as ChatGPT and GitHub Copilot in order to improve the quality of the language and of the code {\color{black} (esp. automatic PR reviews, documentation, creation of figures)}. After using this tool/service, the author reviewed and edited the content as needed and takes full responsibility for the content of the publication