The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
168,791 characters
Online Multivariate Regularized Distributional Regression for High-dimensional Probabilistic Electricity Price Forecasting
\begin{frontmatter}
\title{Online Multivariate Regularized Distributional Regression for High-dimensional Probabilistic Electricity Price Forecasting}
\author[sk,ude]{Simon Hirsch\textsuperscript{\orcidlink{0009-0008-1409-9677}}}
\ead{[email removed]}
\ead{[email removed]}
\affiliation[sk]{
organization={Statkraft Trading GmbH},
country={Germany},
}
\affiliation[ude]{
organization={University of Duisburg-Essen, House of Energy Markets and Finance},
country={Germany},
}
\date{This Version: March 24, 2026}
\begin{abstract}
\noindent {\color{black}
Probabilistic electricity price forecasting (PEPF) is vital for short-term electricity markets, yet the multivariate nature of day-ahead prices — spanning 24 consecutive hours — remains underexplored. At the same time, real-time decision-making requires methods that are both accurate and fast.
We introduce an online algorithm for multivariate distributional regression models, allowing efficient modeling of the conditional means, variances, and dependence structures of electricity prices. The approach combines multivariate distributional regression with online coordinate descent and LASSO-type regularization (absolute shrinkage and selection operator), enabling scalable estimation in high-dimensional covariate spaces. Additionally, we propose a regularized estimation path over increasingly complex dependence structures, allowing for early stopping and avoiding overfitting.
{\color{black}
In a case study using historical data from the German day-ahead market, the proposed method yields interpretable and well-calibrated joint prediction intervals for the 24-dimensional price distribution and provides robust performance across a range of proper scoring rules. The results underscore the importance of modeling the dependence structure of electricity prices.
}
Furthermore, we analyze the trade-off between predictive accuracy and computational costs for batch and online estimation and provide a high-performing open-source Python implementation in the \texttt{ondil} package.
}
\end{abstract}
\begin{keyword}
online learning \sep
multivariate distributional regression \sep
probabilistic electricity price forecasting \sep
LASSO regularization \sep
day-ahead electricity market
\end{keyword}
\end{frontmatter}
\section{Introduction}
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 \cite[see e.g.][]{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 \citep{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.
\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{figures/price_time_series_fig_in_fig.png}
\caption{Time series plot for day-ahead electricity prices $P_{d, h}$ in Germany. {\color{black} The main panel shows the time series of electricity prices. The inset panels plot two different slicing approaches: Plotting 24 time series for each delivery hour $h$ (left) and plotting the cross-section of 24 hourly prices for each day $d$ (right). The joint movement over days and the temporal correlation along the dimension of the delivery hours $h = 0 ,..., 23$ is visible.}}
\label{fig:electricity_prices_spot_timeseries}
\end{figure}
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 \citep[GARCH,][]{nowotarski2018recent, bille2023forecasting, marcjasz2023distributional},} conformal prediction methods, \citep[{\color{black} CP}, see e.g.][]{kath2021conformal, zaffran2022adaptive, lipiecki2024postprocessing, brusaferri2024line}, distributional regression and neural network approaches \citep[e.g.][]{ziel2021gamlss, marcjasz2023distributional, brusaferri2024nbmlss, hirsch2024online}. However, these works treat each delivery hour as \emph{independent}, univariate time series as in \cite{ziel2018day}. {\color{black}
We motivate the need for multivariate probabilistic forecasting approaches for the day-ahead electricity price by two simple plots. Figure \ref{fig:electricity_prices_spot_timeseries} 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{fig:electricity_prices_spot_correlation} 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 \citep[{\color{black} LEAR, upper triangular, }][]{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, \cite{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 \cite{pena2024hedging, lohndorf2023value}.
\begin{wrapfigure}{r}{0.5\textwidth}
\centering
\includegraphics[width=\linewidth]{figures/correlation.png}
\caption{Correlation Matrix for day-ahead electricity prices $P_{d, h}$ in Germany. The lower triangle gives the hourly Pearson correlation $\rho$ for electricity prices. The upper triangle gives the hourly correlation of residuals $\varepsilon_{d, h} = P_{d, h} - \hat{\mu}_{d,h}$ for a standard LASSO-estimated autoregressive model \cite[LEAR, see e.g.][and Eq. \ref{eq:mean_model}]{lago2021forecasting}. The high degree of residual correlation, especially around the noon hours is visible. All correlation coefficients are statistically significant to the $\alpha = 0.01$ confidence level.}
\label{fig:electricity_prices_spot_correlation}
\end{wrapfigure}
Work on \emph{multivariate} probabilistic forecasting for day-ahead electricity prices is sparse in the literature and the majority of the existing work, e.g. \cite{maciejowska2024multiple, berrisch2024multivariate, han2023probabilistic, mashlakov2021assessing} and \cite{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. \cite{janke2020probabilistic} approach the issue through implicit generative Copula models. \cite{grothe2023point} employ the Schaake shuffle, a post-processing method for point forecasts. On the contrary, in the fields of probabilistic weather, renewable production \citep{bjerregaard2021introduction, soerensen2022recent, kolkmann2024modeling} and probabilistic load forecasting \citep{gioia2022additive, browell2022covariance} truly multivariate forecasting approaches have gained more attention.
The goal of distributional regression or ``regression beyond the mean" \citep{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 \citep[GAMLSS, ][]{rigby2005generalized}, of which numerous extensions have been developed over the last years \citep{kock2023truly,kneib2023rage,muschinski2022cholesky} and distributional deep neural networks \citep[DDNN, e.g.][]{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 \citep{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 \cite{umlauf2025scalable} and \cite{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{fig:online_forecasting_study}. Online learning for {\color{black} the absolute shrinkage and selection operator (LASSO) regularized regression} for the mean has been introduced in \cite{angelosante2009online, angelosante2010online} and \cite{messner2019online}. Univariate approaches suitable for probabilistic forecasting based stochastic gradient descent have been developed for specific distributions, \cite[see e.g.][]{pierrot2021adaptive}, conformal prediction \citep[see e.g.][]{zaffran2022adaptive,gibbs2021adaptive,gibbs2024conformal} and the generic online distributional regression in \cite{hirsch2024online}. However, in the multivariate case, the literature remains sparse and focused on unconditional distributions and Copulae \citep[see e.g.][]{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:
\begin{itemize}
\item \emph{Methodological Contribution:} We develop a regularized online estimation method for multivariate distributional regression based on the univariate work by \cite{hirsch2024online}. Our algorithm allows for two layers of regularization: \begin{itemize}
\item By leveraging online coordinate descent and LASSO-type penalties for each individual distribution parameter, we allow for high-dimensional covariate spaces.
\item By exploiting the structure in scale matrix of the multivariate distribution, we develop a path-based estimation along increasingly complex dependency structures, allowing for parsimonious estimation and early stopping. We validate the trade-off between model complexity and predictive accuracy in our application study.
\end{itemize}
Our algorithm is generic and can be applied to any parametric multivariate distribution, as long as the (log-)likelihood function and its derivatives are available. We implement the multivariate normal and $t$-distributions with three different parameterizations of the scale matrix. Further distributions can be easily added.
\item \emph{Applied Contribution:} We apply the method to probabilistic forecasting of the joint distribution of spot electricity prices. We benchmark multivariate distributional regression models to LASSO models, conformal prediction, GARCH, univariate distributional regressions approaches, partly combined with Copula constructions. \begin{itemize}
\item {\color{black} We provide the first comprehensive study on multivariate probabilistic forecasting of day-ahead electricity prices. We evaluate the full multivariate distribution through (strictly) proper scoring rules and the calibration of the 24-dimensional joint prediction bands (JPB) for a wide range of benchmark models and including the volatile years 2021 to 2023 in our test set.}
\item {\color{black} We show that the multivariate distributional regression models, which allow modeling all distributional parameters, i.e. the mean, but also the dependence structure, conditional on explanatory variables such as renewable in-feed or past prices provide well-calibrated prediction bands and superior forecasting performance in terms of the Log Score (LS), Variogram Score (VS) and Dawid-Sebastiani Score (DSS) compared to advanced benchmark models.}
\item We showcase the computational advantage of online estimation by benchmarking with repeated batch estimation on various re-estimation frequencies. These results provide an ``efficient frontier" of computational costs against predictive accuracy.
\end{itemize}
\item \emph{Software:} We provide a high-performing Python implementation using just-in-time compilation and providing a familiar, \texttt{scikit-learn}-like API to facilitate the usage of our package for other researchers. We contributed our code to the \texttt{ondil} package by \cite{hirsch2024online}.
\end{itemize}
Thereby, our contribution is valuable both from a methodological and applied perspective in research and industry alike. Reproduction code can be found on \texttt{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{sec:online_multivariate_distreg} introduces the multivariate, online, regularized distributional regression model. Section \ref{sec:forecasting_study} introduces the forecasting study, the used data and discusses multivariate forecast evaluation in detail. Section~\ref{sec:results} presents our results and analyzes the computational costs. {Section \ref{sec:conclusion}} concludes the paper.}
\section{Online Multivariate Distribution Regression}\label{sec:online_multivariate_distreg}
{\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{fig:paper_overview}. First, we like to introduce some intuition to distributional regression. Distributional regression or ``regression beyond the mean" \citep{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.
\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{sketch_overview.png}
\caption{Online Distributional Regression. For univariate distributions (top panel), online distributional regression has been introduced by \citep{hirsch2024online} by a combination of the iteratively reweighted least squares algorithm with online coordinate descent. We extend the approach to multivariate distributions by utilizing unconstrained scale matrix parameterizations and introducing a path-based regularization along increasingly complex dependency structures.}
\label{fig:paper_overview}
\end{figure}
{\color{black} Section~\ref{sec:distributional_regression_setting} 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{sec:covariance_parameterization}. Section \ref{sec:online_estimation} 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{sec:path_based_estimation}.} }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} \subsection{Distributional Regression Setting} \label{sec:distributional_regression_setting}}
{\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: \begin{equation}
g_k(\ensuremath{\bm{\mathrm{\theta}}}_k) = \ensuremath{\bm{\mathrm{\eta}}}_k = {\ensuremath{\bm{\mathrm{X}}}_k}\ensuremath{\bm{\mathrm{\beta}}}_k
\end{equation}
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: \begin{equation}
y_i \sim \mathcal{F}(\theta_{1i}, ..., \theta_{Ki}) \quad \text{and} \quad \theta_{ki} = g^{-1}\left(\ensuremath{{\ensuremath{\bm{\mathrm{x}}}_{ki}}^\top}\beta_k\right)
\end{equation}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 \citep{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})$. \cite{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 \citep[see e.g. p. 113 in][]{stasinopoulos2024generalized}. Different flavors of LASSO-type regularized estimation approaches have been introduced by \cite{groll2019lasso, muniain2020probabilistic, ziel2021gamlss, o2023variable}. A regularized, incremental estimation approach using online coordinate descent has been proposed by \cite{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 \begin{equation}
g_{km}(\ensuremath{\bm{\mathrm{\theta}}}_{km}) = \ensuremath{\bm{\mathrm{\eta}}}_{km} = \ensuremath{\bm{\mathrm{X}}}_{km}\ensuremath{\bm{\mathrm{\beta}}}_{km}.
\end{equation}
This formulation is rather general. {\color{black} We note two important points here:
\begin{itemize}
\item The distribution parameters are subject to constraints. This especially holds for the scale (respectively precision) matrix, which is often required to be positive semi-definite. To ensure this holds, unconstrained parameterizations such as the (modified) Cholesky-decomposition ({\color{black} CD resp. MCD}) or the low-rank approximation are often used \citep[{\color{black} LRA}][]{pourahmadi2011covariance, muschinski2022cholesky, salinas2019high}. The parameterization of the scale matrix will be discussed in Section \ref{sec:covariance_parameterization}.
\item In practice, the different distribution parameters $\ensuremath{\bm{\mathrm{\theta}}}_{1}, ..., \ensuremath{\bm{\mathrm{\theta}}}_{K}$ can have many different shapes. Take, e.g. the multivariate $t$-distribution, parameterized using the Cholesky factor of the precision matrix $\ensuremath{\bm{\mathrm{\Sigma}}} = \ensuremath{\ensuremath{\bm{\mathrm{L}}}^\top}\ensuremath{\bm{\mathrm{L}}}$, denoted as $t_D(\ensuremath{\bm{\mathrm{\theta}}}_1, \ensuremath{\bm{\mathrm{\theta}}}_2, \ensuremath{\bm{\mathrm{\theta}}}_3) \Leftrightarrow t_D(\ensuremath{\bm{\mathrm{\mu}}}, \ensuremath{\bm{\mathrm{L}}}, \ensuremath{\bm{\mathrm{\nu}}})$. Then $\ensuremath{\bm{\mathrm{\mu}}}$ is a ${N \times D}$ matrix, $\ensuremath{\bm{\mathrm{L}}}$ is a ${N \times D \times D}$ cube (of which each vertical slice is a triangular matrix) and $\ensuremath{\bm{\mathrm{\nu}}}$ is a $N \times 1$ vector. Accordingly, the index set $\mathcal{M}_k$ of coordinates spans $\mathcal{M}_1 = \{1, .., D\}$, $\mathcal{M}_2 = \{(1, 1), ..., (D, D)\}$ and $\mathcal{M}_3 = \{1\}$ and its cardinality is given by the product of the parameter's dimension beyond $N$.
\end{itemize}} The general setting introduced here includes the Gaussian multivariate distributional regression introduced by \cite{muschinski2022cholesky}, the Copula-based multivariate distributional regression by \cite{kock2023truly} and the MCD-based additive covariance models by \cite{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{sec:online_estimation} and the following Section \ref{sec:covariance_parameterization} briefly discusses different options to parameterize the covariance respectively precision matrix.
\subsection{Parameterization of the Precision Matrix}\label{sec:covariance_parameterization}
{\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. \cite{pourahmadi2011covariance, muschinski2022cholesky, kock2023truly}. Additionally, we employ a low-rank approximation of the precision matrix, which has been used by \cite{salinas2019high} in the context of high-dimensional Gaussian processes and by \cite{marz2022multi} and \cite{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:}
\begin{align}
\ensuremath{\bm{\mathrm{\Sigma}}} = \ensuremath{\bm{\mathrm{A}}}\ensuremath{\ensuremath{\bm{\mathrm{A}}}^\top} && {\color{black} \text{and}} && \ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{(\ensuremath{\bm{\mathrm{A}}}^{-1})^\top}(\ensuremath{\bm{\mathrm{A}}}^{-1})
\end{align}
\cite{muschinski2022cholesky} parametrize the normal distribution in terms of~$\ensuremath{\bm{\mathrm{A}}}^{-1}$ and \cite{kock2023truly} choose~$\ensuremath{\bm{\mathrm{A}}}$. {\color{black} Additionally the modified Cholesky-decomposition (MCD) can be used:
\begin{align}
\ensuremath{\bm{\mathrm{\Sigma}}} = \ensuremath{(\ensuremath{\bm{\mathrm{L}}}^{-1})^\top}\ensuremath{\bm{\mathrm{D}}}\ensuremath{\bm{\mathrm{L}}}^{-1} && {\color{black} \text{and}} && \ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{\ensuremath{\bm{\mathrm{L}}}^\top}(\ensuremath{\bm{\mathrm{D}}}^{-1})\ensuremath{\bm{\mathrm{L}}}
\end{align}}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 \begin{equation}
\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{\bm{\mathrm{A}}} + \ensuremath{\bm{\mathrm{V}}}\ensuremath{\ensuremath{\bm{\mathrm{V}}}^\top},
\end{equation}
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}}}$.
\subsection{Online Estimation Algorithm}\label{sec:online_estimation}
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. \citet[{\color{black} RS}]{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
\begin{equation}\label{eq:rs_score_vector}
\ensuremath{\bm{\mathrm{u}}} = \frac{\partial \ell}{\partial \eta},
\end{equation}
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
\begin{equation}\label{eq:rs_working_vector}
\ensuremath{\bm{\mathrm{z}}} = g\big(\hat{\ensuremath{\bm{\mathrm{\theta}}}}\big) + \frac{\partial \ell}{\partial \eta} \ensuremath{\bm{\mathrm{W}}}^{-1} \Leftrightarrow \ensuremath{\bm{\mathrm{z}}} = \ensuremath{\bm{\mathrm{\eta}}} + \frac{\partial \ell}{\partial \eta} \ensuremath{\bm{\mathrm{W}}}^{-1},
\end{equation}
where the weights are defined as:
\begin{align}\label{eq:rs_weights}
\ensuremath{\bm{\mathrm{W}}} = - \frac{\partial^2 \ell}{\partial \eta^2} & & \text{or} & & \ensuremath{\bm{\mathrm{W}}} = - \mathbb{E}\left[\frac{\partial^2 \ell}{\partial \eta^2} \right]
\end{align}
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}
\cite{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{alg:online_lasso}. 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{sec:online_coordinate_descent} and \cite{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{alg:online_mv_gamlss} 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{sec:distributional_regression_setting}, 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.
\begin{algorithm}[htb]
\raggedright
\caption{Online regularized multivariate distributional regression.}
\label{alg:online_mv_gamlss}
\DontPrintSemicolon
\KwIn{$\ensuremath{\bm{\mathrm{y}}}^{[n+1]}, \ensuremath{\bm{\mathrm{X}}}_{k,m}^{[n+1]}$ and the stored Gramian matrices $\ensuremath{\bm{\mathrm{G}}}_{km}^{[n]}, \ensuremath{\bm{\mathrm{h}}}_{km}^{[n]}$.}
Initialize the fitted values $\widehat{\ensuremath{\bm{\mathrm{\theta}}}}_{km}^{[n+1]} = \hat{\ensuremath{\bm{\mathrm{\beta}}}}_{km}^{[n]}\ensuremath{(\ensuremath{\bm{\mathrm{X}}}_{km}^{[n+1]})^\top}$ for $k,m \in \mathcal{K} \times \mathcal{M}$. \;
Evaluate the linear predictors ${\widehat{\ensuremath{\bm{\mathrm{\eta}}}}^{[n+1]}_{km} = g_{km}(\hat{\theta}_{km}^{[n+1]})}$ for $k,m \in \mathcal{K} \times \mathcal{M}$. \;
\For{\normalfont outer iteration $i = 0, ...$ \normalfont\textbf{until} convergence}{
\ForAll{\normalfont distribution parameter $k \in \mathcal{K}$}{
\For{\normalfont inner iteration $r = 0, 1, ...$ \normalfont\textbf{until} convergence}{
\ForAll{\normalfont distribution element $m \in \mathcal{M}_k$}{
Evaluate $u_{km}^{[n+1]}$, $w_{km}^{[n+1]}$ and $z_{km}^{[n+1]}$ using Equations \eqref{eq:rs_score_vector}, \eqref{eq:rs_working_vector} and \eqref{eq:rs_weights}. \;
Update $\ensuremath{\bm{\mathrm{G}}}_{km}^{[n+1]}$ and $\ensuremath{\bm{\mathrm{h}}}_{km}^{[n+1]}$ based on $u_{km}^{[n+1]}$, $w_{km}^{[n+1]}$ and $z_{km}^{[n+1]}$. \;
Update $\hat{\ensuremath{\bm{\mathrm{\beta}}}}_{km}^{[n+1, \ensuremath{\bm{\mathrm{\lambda}}}]} \gets \textsc{OnlineLASSO}(\ensuremath{\bm{\mathrm{G}}}_{km}^{[n+1]}, \ensuremath{\bm{\mathrm{h}}}_{km}^{[n+1]}, \ensuremath{\bm{\mathrm{\lambda}}})$. \;
Model selection using IC for $\hat{\ensuremath{\bm{\mathrm{\beta}}}}_{km}^{[n+1]}$ and calculate the updated $\hat{\ensuremath{\bm{\mathrm{\eta}}}}_{km}^{[n+1]}$ \;
}
}
}
}
\KwOut{$\widehat{\ensuremath{\bm{\mathrm{\beta}}}}_{k, n+1}$ and $\widehat{\ensuremath{\bm{\mathrm{\Theta}}}}^{[n+1]} = (\widehat{\ensuremath{\bm{\mathrm{\theta}}}}_0^{[n+1]}, ..., \widehat{\ensuremath{\bm{\mathrm{\theta}}}}_p^{[n+1]})$ and the updated $\ensuremath{\bm{\mathrm{G}}}_{km}^{[n+1]}$ and $\ensuremath{\bm{\mathrm{h}}}_{km}^{[n+1]}$.}
\end{algorithm}
}
{\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, \cite{rigby2005generalized} use Fisher's scoring {\color{black} in Equation \ref{eq:rs_weights}}. 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 \emph{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 \cite[see e.g.][]{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 \emph{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
\begin{align}
\frac{\partial \ell}{\partial \eta} &= \frac{\partial \ell}{\partial \theta} \left( \frac{\partial g(\theta)}{\partial\theta}\right)^{-1} \quad\label{eq:newton_raphson_generic_first} \quad \text{and}\\
\frac{\partial^2 \ell}{\partial \eta^2} &= \left(
\frac{\partial^2 \ell}{\partial \theta^2} \frac{\partial g(\theta)}{\partial\theta} - \frac{\partial \ell}{\partial \theta} \frac{\partial^2 g(\theta)}{\partial\theta^2}
\right) \left( \frac{\partial g(\theta)}{\partial\theta}\right)^{-3}. \label{eq:newton_raphson_generic_second}
\end{align}
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{tab:implemented_dists}. The derivations can be found in {\color{black} the online supplement}.
\begin{table}[htb]
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{lllllll}
\toprule
Distribution & \multicolumn{2}{c}{Location} & \multicolumn{2}{c}{Scale $\boldsymbol{\Sigma}$ resp. Precision $\boldsymbol{\Omega}$} & \multicolumn{2}{c}{Shape} \\
\cmidrule(lr){2-3}
\cmidrule(lr){4-5}
\cmidrule(lr){6-7}
& Param. & Dim. & Param. & Dim. & Param. & Dim. \\
\midrule
Multivariate Gaussian & $\ensuremath{\bm{\mathrm{\mu}}}$ & $N \times D$ & $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{(\ensuremath{\bm{\mathrm{A}}}^{-1})^\top}(\ensuremath{\bm{\mathrm{A}}}^{-1})$ & $N \times \text{triangular}(D \times D)$ & - & - \\
Multivariate Gaussian & $\ensuremath{\bm{\mathrm{\mu}}}$ & $N \times D$ & $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{\ensuremath{\bm{\mathrm{L}}}^\top}(\ensuremath{\bm{\mathrm{D}}}^{-1})\ensuremath{\bm{\mathrm{L}}}$ & $N \times \text{triangular}(D \times D)$, $N \times \text{diag}(D)$ & - & -\\
Multivariate Gaussian & $\ensuremath{\bm{\mathrm{\mu}}}$ & $N \times D$ & $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{\bm{\mathrm{A}}} + \ensuremath{\bm{\mathrm{V}}}\ensuremath{\ensuremath{\bm{\mathrm{V}}}^\top}$ & $N \times \text{diag}(D)$, $D \times r $ & - & -\\
Multivariate-$t$ & $\ensuremath{\bm{\mathrm{\mu}}}$ & $N \times D$ & $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{(\ensuremath{\bm{\mathrm{A}}}^{-1})^\top}(\ensuremath{\bm{\mathrm{A}}}^{-1})$ & $N \times \text{triangular}(D \times D)$ & $\nu$ & $N \times 1$\\
Multivariate-$t$ & $\ensuremath{\bm{\mathrm{\mu}}}$ & $N \times D$ & $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{\ensuremath{\bm{\mathrm{L}}}^\top}(\ensuremath{\bm{\mathrm{D}}}^{-1})\ensuremath{\bm{\mathrm{L}}}$ & $N \times \text{triangular}(D \times D)$, $N \times \text{diag}(D)$ & $\nu$ & $N\times 1$ \\
Multivariate-$t$ & $\ensuremath{\bm{\mathrm{\mu}}}$ & $N \times D$ & $\ensuremath{\bm{\mathrm{\Omega}}} = \ensuremath{\bm{\mathrm{A}}} + \ensuremath{\bm{\mathrm{V}}}\ensuremath{\ensuremath{\bm{\mathrm{V}}}^\top}$ & $N \times \text{diag}(D)$, $D \times r $ & $\nu$ & $N\times 1$ \\
\bottomrule
\end{tabular}}
\caption{Overview of multivariate distributions and scale matrix parameterization (Param.) implemented in the paper and the respective dimensions (Dim.) for input data $\ensuremath{\bm{\mathrm{Y}}}$ of shape $N \times D$. Note that the number of parameters for the CD-based parameterization grows quadratically in $D$, but the LRA-based parameterizations grow linear in $D$ for fixed $r$. }
\label{tab:implemented_dists}
\end{table}
{\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 \cite{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 \begin{equation}
\operatorname{IC} =
-2\ell\left(\ensuremath{\bm{\mathrm{Y}}} \mid \widehat{\ensuremath{\bm{\mathrm{\Theta}}}}\right)
+ \nu_0 K
+ \nu_1 K \log(N)
+ \nu_2 K \log\left(\log(N)\right)
\end{equation}
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
\begin{equation}
\ell\left(\ensuremath{\bm{\mathrm{Y}}} \mid \widehat{\ensuremath{\bm{\mathrm{\Theta}}}}^{[\lambda_i]} \right) \approx
\ell\left(\ensuremath{\bm{\mathrm{Y}}} \mid \widehat{\ensuremath{\bm{\mathrm{\Theta}}}}^{[\lambda_0]}\right)
+ \frac{\partial \ell}{\partial \theta} \left( \widehat{\theta}_{km}^{[\lambda_0]} - \widehat{\theta}_{km}^{[\lambda_i]}\right)
\end{equation}
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{tab:implemented_dists} 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 \begin{equation}\label{eq:dampening}
\hat{\eta}_m^{[0, i]} \leftarrow g_m^{-1}\left((i+1)\hat{\theta}_m^{[0, i]} + \hat{\theta}_m^{[0, i-1]}) / (i+1)\right)
\end{equation} 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{eq:dampening}), which, in combination with a Cole-Green-algorithm \citep[CG, see][]{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.
}
\subsection{Path-based Regularized Estimation for the Scale Matrix}\label{sec:path_based_estimation}
{\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.
\begin{figure}[htb]
\centering
\resizebox{\linewidth}{!}{\begin{tikzpicture}[
align=center,
squarednodegreen/.style={rectangle, draw=Gray, fill=Gray!01, very thick, minimum size=5mm},
squarednodeblue/.style={rectangle, draw=Gray, fill=Gray!01, very thick, minimum size=5mm},
]
\node[squarednodegreen] (cholfac_00) {
1st iteration $\ensuremath{\bm{\mathrm{L}}}^{[\textcolor{Red}{0}]}$ \\[1em]
$\begin{pmatrix}
l_{1,1} & & & & 0\\
0 & l_{2,2} & & \\
... & 0 & l_{3, 3} & & \\
& ...& ... & ... & \\
0 & & ... &0 & l_{d, d}
\end{pmatrix}$};
\node[squarednodegreen] (cholfac_01) [right=of cholfac_00] {
2nd iteration $\ensuremath{\bm{\mathrm{L}}}^{[\textcolor{Red}{1}]}$ \\[1em]
$\begin{pmatrix}
l_{1,1} & & & & 0\\
\color{Red}{l_{2,1}} & l_{2,2} & & \\
& \color{Red}{l_{3,2}} & l_{3, 3} & & \\
&& \color{Red}{...} & ... & \\
0&&& \color{Red}{l_{d,d-1}}& l_{d, d}
\end{pmatrix}$};
\node[squarednodegreen] (cholfac_02) [right=of cholfac_01] {
3rd iteration $\ensuremath{\bm{\mathrm{L}}}^{[\textcolor{Red}{2}]}$ \\[1em]
$\begin{pmatrix}
l_{1,1} & & & & 0\\
\color{Green}{l_{2,1}} & l_{2,2} & & \\
\color{Red}{l_{3,1}} & \color{Green}{l_{3,2}} & l_{3, 3} & & \\
& \color{Red}{...} & \color{Green}{...} & ... & \\
0&& \color{Red}{l_{d, d-2}} & \color{Green}{l_{d,d-1}}& l_{d, d}
\end{pmatrix}$};
\node[squarednodegreen] (cholfac_03) [right=of cholfac_02] {
Last iteration $\ensuremath{\bm{\mathrm{L}}}^{[\textcolor{Red}{d-1}]}$\\[1em]
$\begin{pmatrix}
l_{1,1} & & & & 0\\
\color{Green}{l_{2,1}} & l_{2,2} & & \\
\color{Green}{l_{3,1}} & \color{Green}{l_{3,2}} & l_{3, 3} & & \\
\color{Green}{...} & \color{Green}{...} & \color{Green}{...} & ... & \\
\color{Red}{l_{d,d}} & \color{Green}{...} & \color{Green}{l_{d, d-2}} & \color{Green}{l_{d,d-1}}& l_{d, d}
\end{pmatrix}$};
\node[squarednodeblue] (lowrank0) [below=2cm of cholfac_00] {
1st iteration $\ensuremath{\bm{\mathrm{V}}}^{[\textcolor{Red}{0}]}$ \\[1em]
$\begin{pmatrix}
0 & ... & 0 \\
0 & ... & 0 \\
... & ... & ... \\
0 & ... & 0 \\
\end{pmatrix}$};
\node[squarednodeblue] (lowrank1) [below=2cm of cholfac_01] {
2nd iteration $\ensuremath{\bm{\mathrm{V}}}^{[\textcolor{Red}{1}]}$ \\[1em]
$\begin{pmatrix}
\color{Red}{v_{1,1}} & 0 & ... & 0 \\
\color{Red}{v_{2,1}} & 0 & ... & 0 \\
\color{Red}{...} & ... & ... & ... \\
\color{Red}{v_{d,1}} & 0 & ... & 0 \\
\end{pmatrix}$};
\node[squarednodeblue] (lowrank2) [below=2cm of cholfac_02] {
3rd iteration $\ensuremath{\bm{\mathrm{V}}}^{[\textcolor{Red}{2}]}$ \\[1em]
$\begin{pmatrix}
\color{Green}{v_{1,1}} & \color{Red}{v_{2,1}} & 0 &... & 0 \\
\color{Green}{v_{2,1}} & \color{Red}{v_{2,2}} & 0 &... & 0 \\
\color{Green}{...} & \color{green}{...} & .. &... & ... \\
\color{Green}{v_{d,1}} & \color{Red}{v_{d,2}} & 0 &... & 0 \\
\end{pmatrix}$};
\node[squarednodeblue] (lowrank3) [below=2cm of cholfac_03] {
Last iteration $\ensuremath{\bm{\mathrm{V}}}^{[\textcolor{Red}{d}]}$\\[1em]
$\begin{pmatrix}
\color{Green}{v_{1,1}} & \color{Green}{v_{2,1}} & ... & \color{Green}{v_{1,r-1}} & \color{Red}{v_{2,r}} \\
\color{Green}{v_{2,1}} & \color{Green}{v_{2,2}} & ... & \color{Green}{v_{2,r-1}} & \color{Red}{v_{2,r}} \\
\color{Green}{...} & \color{Green}{...} & ... & ... & \color{Red}{...} \\
\color{Green}{v_{d,1}} & \color{Green}{v_{d,2}} & ... & \color{Green}{v_{d,r-1}} & \color{Red}{v_{2,r}}\\
\end{pmatrix}$};
\draw[very thick][->] (cholfac_00.east) -- (cholfac_01.west);
\draw[very thick][->] (cholfac_01.east) -- (cholfac_02.west);
\draw[very thick, dotted][->] (cholfac_02.east) -- (cholfac_03.west);
\node[draw, very thick, draw=gray!25, fit=(cholfac_00) (cholfac_01) (cholfac_02) (cholfac_03), inner sep=10pt, minimum width=3cm, minimum height=2cm, label=above:{\large Cholesky-based scale matrix parameterization $\ensuremath{\bm{\mathrm{\Omega}}}^{[\alpha]} = \ensuremath{(L^{[\alpha]})^\top} (L^{[\alpha]})$}] (box1) {};
\draw[very thick][->] (lowrank0.east) -- (lowrank1.west);
\draw[very thick][->] (lowrank1.east) -- (lowrank2.west);
\draw[very thick, dotted][->] (lowrank2.east) -- (lowrank3.west);
\node[draw, very thick, draw=gray!25, fit=(lowrank0) (lowrank1) (lowrank2) (lowrank3), inner sep=10pt, minimum width=3cm, minimum height=2cm, label=above:{\large Low-rank approximation-based scale matrix parameterization, keeping $\ensuremath{\bm{\mathrm{A}}} = \operatorname{diag}(a_0, ..., a_d)$}] (box2) {};
\end{tikzpicture}}
\caption{
Path-based estimation along increasingly complex scale matrix parameterizations.
{\color{black}
The top panel shows the AD-$r$ regression for a Cholesky-based parameterization $\ensuremath{((\ensuremath{\bm{\mathrm{L}}}^{[\alpha]})^{-1})^\top}(\ensuremath{\bm{\mathrm{L}}}^{[\alpha]})^{-1}$.
The superscript $[\alpha]$ denotes the iteration regularization.
Black elements are the state after the initial fit assuming independence.
{\color{black} Red elements of the vectors are added in iteration $\alpha$.
Every arrow denotes that we add one off-diagonal respectively column to the parameterization of the scale matrix.
Green elements have been added in previous iterations ($\alpha-1$ and before) and are warm-started for the current iteration.}
The lower panel shows the estimation along the LRA-based parameterization $\ensuremath{\bm{\mathrm{A}}} + (\ensuremath{\bm{\mathrm{V}}}^{[\alpha]})\ensuremath{(\ensuremath{\bm{\mathrm{V}}}^{[\alpha]})^\top}$,
}
where $\ensuremath{\bm{\mathrm{A}}} = \operatorname{diag}(a_1, ..., a_d)$ is not regularized and the $D \times r$ matrix $\ensuremath{\bm{\mathrm{V}}}$ is filled column-wise with non-zero elements.
}
\label{fig:path_based_estimation}
\end{figure}
{\color{black} This approach is inspired by path-based estimation methods in high-dimensional statistics \citep[such as the graphical LASSO,][]{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 \citep{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. \cite{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{fig:path_based_estimation} 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$, \begin{enumerate}
\item Fit the distributional regression model using Algorithm \ref{alg:online_mv_gamlss} for the current regularization level $\alpha$, where we set the elements of the scale to zero for \begin{itemize}
\item the {\color{black} Cholesky}-based parameterizations if the indices {\color{black} of the diagonal matrices~$\ensuremath{\bm{\mathrm{A}}}$ or~$\ensuremath{\bm{\mathrm{L}}}$,} $i,j$ are such that $|i-j| > \alpha$.
\item the LRA-based parameterization if the indices $d,r$ of $\ensuremath{\bm{\mathrm{V}}}$ are such that $r \geq \alpha$.
\end{itemize}
\item Evaluate the log-likelihood for the current regularization level $\alpha$.
\item Stop early if the log-likelihood (or information criteria) does not increase sufficiently.
\end{enumerate}
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:} \begin{itemize}
\item For a (small) fixed maximum regularization size, the number of parameters in the CD {\color{black} and MCD}-based parameterization grows (almost) linearly in $D$, alleviating the disadvantage of quadratic complexity.
\item For the multivariate $t$-distribution, independence is only achieved as $\nu \rightarrow \infty$. We therefore set a high initial guess ($\nu = 10^6$) for the first outer iteration of $\ensuremath{\bm{\mathrm{\mu}}}$ and $\ensuremath{\bm{\mathrm{\Omega}}}$ to ensure numerical stability for the first iteration and subsequently choose a lower initial guess for the first iteration of $\nu$, since the Newton-Raphson algorithm relies on appropriate start values and tends to alternate between extrema otherwise \cite[see e.g.][on the impact of initial values for Newton-Raphson algorithms.]{casella2021choice,kornerup2006choosing}\footnote{We have found the algorithm to iterate between $\nu = 2$ and $\nu > 10^{10}$ for too large start values for the degrees of freedom. The proposed approach however has proved stable through the full simulation study with highly volatile electricity prices.}
\item We can employ the path-based estimation to early stop the estimation, if the log-likelihood {\color{black} or} an information criterion does not increase sufficiently by adding more non-zero elements. This allows for both, implicit regularization and decreased estimation time. However, {\color{black} once the algorithm stops early}, we cannot increase the complexity of the parameterization in the online estimation, but need to treat it as fixed.
\end{itemize}
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 \cite[see e.g. the adaptive block structure in][]{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{sec:overfitting} and leave the development of more advanced selection schemes for future research.}
\section{Forecasting Study}\label{sec:forecasting_study}
\subsection{Day-ahead Electricity Market and Data}
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{fig:day_ahead_market}. 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 \citep{bille2023forecasting,hirsch2024online,viehmann2017state}.
\begin{figure}[htb]
\centering
\resizebox{\linewidth}{!}{\begin{tikzpicture}[
node_delivery_f1/.style={regular polygon,regular polygon sides=4, draw=Red, fill=Red!10, very thick, minimum size=12mm},
node_delivery_f2/.style={regular polygon,regular polygon sides=4, draw=Blue, fill=Blue!10, very thick, minimum size=12mm},
node_delivery_f3/.style={regular polygon,regular polygon sides=4, draw=Green, fill=Green!10, very thick, minimum size=12mm},
node_delivery_e/.style={regular polygon,regular polygon sides=4, draw=White, fill=White, very thick, minimum size=12mm},
node_auction_f1/.style={circle, draw=Red, fill=Red!10, very thick, minimum size=10mm},
node_auction_f2/.style={circle, draw=Blue, fill=Blue!10, very thick, minimum size=10mm},
node_auction_f3/.style={circle, draw=Green, fill=Green!10, very thick, minimum size=10mm},
]
\node[node_delivery_e] (delivery_f0) [] {};
\node[node_delivery_e] (delivery_f1) [right=0.25cm of delivery_f0] {};
\node[node_delivery_e] (delivery_f2) [right=0.25cm of delivery_f1] {};
\node[node_delivery_e] (delivery_f3) [right=0.25cm of delivery_f2] {...};
\node[node_auction_f1, label=above:{Auction $t-1$, 12:00}] (auction_0) [above=1cm of delivery_f0] {};
\node[node_delivery_f1] (delivery_00) [right=0.25cm of delivery_f3] {0};
\node[node_delivery_f1] (delivery_01) [right=0.25cm of delivery_00] {1};
\node[node_delivery_f1] (delivery_02) [right=0.25cm of delivery_01] {2};
\node[node_delivery_e, label=below:{Delivery periods $h=0, ..., 23$ on day $t$}] (delivery_0m) [right=0.25cm of delivery_02] {...};
\node[node_delivery_f1] (delivery_03) [right=0.25cm of delivery_0m] {22};
\node[node_delivery_f1] (delivery_04) [right=0.25cm of delivery_03] {22};
\node[node_delivery_f1] (delivery_05) [right=0.25cm of delivery_04] {23};
\node[node_auction_f2, label=above:{Auction $t$, 12:00}] (auction_1) [above=1cm of delivery_0m] {};
\node[node_delivery_f2] (delivery_10) [right=0.25cm of delivery_05] {0};
\node[node_delivery_f2] (delivery_11) [right=0.25cm of delivery_10] {1};
\node[node_delivery_f2] (delivery_12) [right=0.25cm of delivery_11] {2};
\node[node_delivery_e, label=below:{Delivery periods $h=0, ..., 23$ on day $t+1$}] (delivery_1m) [right=0.25cm of delivery_12] {...};
\node[node_delivery_f2] (delivery_13) [right=0.25cm of delivery_1m] {22};
\node[node_delivery_f2] (delivery_14) [right=0.25cm of delivery_13] {22};
\node[node_delivery_f2] (delivery_15) [right=0.25cm of delivery_14] {23};
\node[node_auction_f3, label=above:{Auction $t+1$, 12:00}] (auction_2) [above=1cm of delivery_1m] {};
\node[node_delivery_f3] (delivery_20) [right=0.25cm of delivery_15] {0};
\node[node_delivery_f3] (delivery_21) [right=0.25cm of delivery_20] {1};
\node[node_delivery_f3] (delivery_22) [right=0.25cm of delivery_21] {2};
\node[node_delivery_e, label=below:{Delivery periods $h=0, ..., 23$ on day $t+2$}] (delivery_2m) [right=0.25cm of delivery_22] {...};
\node[node_delivery_f3] (delivery_23) [right=0.25cm of delivery_2m] {22};
\node[node_delivery_f3] (delivery_24) [right=0.25cm of delivery_23] {22};
\node[node_delivery_f3] (delivery_25) [right=0.25cm of delivery_24] {23};
\draw[->, dotted, thick] (auction_0) to (delivery_00.north west);
\draw[->, dotted, thick] (auction_0) to (delivery_01.north west);
\draw[->, dotted, thick] (auction_0) to (delivery_02.north west);
\draw[->, dotted, thick] (auction_0) to (delivery_03.north west);
\draw[->, dotted, thick] (auction_0) to (delivery_04.north west);
\draw[->, dotted, thick] (auction_0) to (delivery_05.north west);
\draw[->, dotted, thick] (auction_1) to (delivery_10.north west);
\draw[->, dotted, thick] (auction_1) to (delivery_11.north west);
\draw[->, dotted, thick] (auction_1) to (delivery_12.north west);
\draw[->, dotted, thick] (auction_1) to (delivery_13.north west);
\draw[->, dotted, thick] (auction_1) to (delivery_14.north west);
\draw[->, dotted, thick] (auction_1) to (delivery_15.north west);
\draw[->, dotted, thick] (auction_2) to (delivery_20.north west);
\draw[->, dotted, thick] (auction_2) to (delivery_21.north west);
\draw[->, dotted, thick] (auction_2) to (delivery_22.north west);
\draw[->, dotted, thick] (auction_2) to (delivery_23.north west);
\draw[->, dotted, thick] (auction_2) to (delivery_24.north west);
\draw[->, dotted, thick] (auction_2) to (delivery_25.north west);
\draw[very thick, ->] (auction_0.east) to node[midway, above] {Results at 12:42, $t-1$} (auction_1.west);
\draw[very thick, ->] (auction_1.east) to node[midway, above] {Results at 12:42, $t$} (auction_2.west);
\end{tikzpicture}}
\caption{Structure of the day-ahead electricity market in Germany. {\color{black} Based on information from the EPEX SPOT website and \cite{viehmann2017state}}.}
\label{fig:day_ahead_market}
\end{figure}
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{fig:online_forecasting_study} 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 \citep[see e.g.][]{lago2021forecasting,nowotarski2018recent} while online learning schemes are a rather recent development in EPF \citep[e.g.][]{berrisch2024multivariate, zaffran2022adaptive,hirsch2024online}.}
\begin{figure}[htb]
\centering
\resizebox{\linewidth}{!}{\begin{tikzpicture}[
node_initial_batch/.style={rectangle, draw=Blue, fill=Blue!10, very thick, minimum size=10mm},
node_initial_batch_last/.style={rectangle, draw=White, fill=Blue!5, very thick, minimum size=10mm},
node_step/.style={rectangle, draw=Green, fill=Green!10, very thick, minimum size=10mm},
node_step_last/.style={rectangle, draw=White, fill=Green!5, very thick, minimum size=10mm},
]
\node[node_initial_batch, minimum width=5cm] (batch_0) [
label={[label distance=0.25cm, align=left] above:\raggedright\text{\Large Batch Expanding Window}}
] {Initial set $i = 0, ..., N$};
\node[node_initial_batch, minimum width=6cm] (batch_1) [below=0.25cm of batch_0.south west, anchor=north west] {Dataset $i = 0, ..., N+1$};
\node[node_initial_batch, minimum width=7cm] (batch_2) [below=0.25cm of batch_1.south west, anchor=north west] {Dataset $i = 0, ..., N+2$};
\node[node_initial_batch, minimum width=8cm] (batch_3) [below=0.25cm of batch_2.south west, anchor=north west] {Dataset $i = 0, ..., N+3$};
\node[node_initial_batch, minimum width=9cm] (batch_4) [below=0.25cm of batch_3.south west, anchor=north west] {Dataset $i = 0, ..., N+4$};
\node[node_initial_batch_last, minimum width=10cm] (batch_5) [below=0.25cm of batch_4.south west, anchor=north west] {...};
\draw[-] (9, 1) to (9, -7.5);
\node[node_initial_batch, minimum width=5cm] (slide_0) [
label={[label distance=0.25cm, align=left] above:\raggedright\text{\Large Batch Rolling Window}},
right=8cm of batch_0
] {Initial set $i = 0, ..., N$};
\node[node_initial_batch, minimum width=5cm] (slide_1) [below right=0.25cm and 1cm of slide_0.south west, anchor=north west] {Dataset $i = 1, ..., N+1$};
\node[node_initial_batch, minimum width=5cm] (slide_2) [below right=0.25cm and 1cm of slide_1.south west, anchor=north west] {Dataset $i = 2, ..., N+2$};
\node[node_initial_batch, minimum width=5cm] (slide_3) [below right=0.25cm and 1cm of slide_2.south west, anchor=north west] {Dataset $i = 3, ..., N+3$};
\node[node_initial_batch, minimum width=5cm] (slide_4) [below right=0.25cm and 1cm of slide_3.south west, anchor=north west] {Dataset $i = 4, ..., N+4$};
\node[node_initial_batch_last, minimum width=5cm] (slide_5) [below right=0.25cm and 1cm of slide_4.south west, anchor=north west] {...};
\draw[-] (22, 1) to (22, -7.5);
\node[node_initial_batch, minimum width=5cm] (online_0) [
label={[label distance=0.25cm, align=left]above:\text{\Large Online Learning}},
right=8cm of slide_0
] {Initial set $i = 0, ..., N$};
\node[node_step] (online_1) [below=0.25cm of online_0.south east, anchor=north west] {$N$+1};
\node[node_step] (online_2) [
label={[label distance=0.25cm, align=left]left:{New observations \\ added row-by-row.}},
below=0.25cm of online_1.south east, anchor=north west] {$N$+2};
\node[node_step] (online_3) [below=0.25cm of online_2.south east, anchor=north west] {$N$+3};
\node[node_step] (online_4) [below=0.25cm of online_3.south east, anchor=north west] {$N$+4};
\node[node_step_last] (online_5) [below=0.25cm of online_4.south east, anchor=north west] {...};
\end{tikzpicture}
}
\caption{Repeated Batch Learning vs. Online Learning for the forecasting study.}
\label{fig:online_forecasting_study}
\end{figure}
{\color{black} We use the same data set as in \cite{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 \citep{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{tab:fundamental_variables}. {\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 \citep{welford1962note}.} Note that, in this application study, we have $T$ corresponding to $N$ and $H$ corresponding to $D$ in the general notation.
\begin{table}[thb]
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{llll}
\toprule
Variable & Description & Resolution & Source \\
\midrule
$\operatorname{ResLoad}_{t,h}$ & Day-ahead residual load forecast & Hourly & ENTSO-E \\
$\overline{\operatorname{ResLoad}}_{t}$ & Day-ahead baseload residual load forecast $\frac{1}{H}\sum_{h=1}^H\operatorname{ResLoad}_{t,h}$ & Daily & ENTSO-E \\
$\operatorname{EUA}_{t}$ & Log European Union Emission Allowance (EUA) prices & Daily & \url{investing.com} \\
$\operatorname{Gas}_{t}$ & Log natural gas prices & Daily & \url{investing.com} \\
$\operatorname{Coal}_{t}$ & Log coal prices & Daily & \url{investing.com} \\
$\operatorname{Oil}_{t}$ & Log oil prices & Daily & \url{investing.com} \\
$\operatorname{WD}_{t}$ & Weekday dummies & Daily & Calendar \\
\bottomrule
\end{tabular}
}
\caption{Variables from the data set of \cite{marcjasz2023distributional} and \cite{lipiecki2024postprocessing}. {\color{black} Wind power production, solar power production and system load forecasts are summarized in the residual load and are retrieved from the European Network of Transmission System Operators for Electricity (ENTSO-E).}}
\label{tab:fundamental_variables}
\end{table}
\subsection{Model Definition and Benchmarks}
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 \emph{adaptive} estimation, which is updating a single, unconditional distributional parameter and the full \emph{conditional} estimation linking the distribution parameter to explanatory variables. {\color{black} Table~\ref{tab:models}} 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 \citep[ACI, see][{\color{black} during the online training}]{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 \texttt{arch} package in Python \citep{kevin_sheppard_2024_14035889}. {\color{black} This is the only model that we cannot run online in the study.}}
\begin{table}
\resizebox{\textwidth}{!}{
\begin{tabular}{p{0.4\textwidth}p{0.95\textwidth}}
\toprule
Abbreviation & Short Description \\
\midrule
LEAR-\{structure\} & LEAR Model (Equation \ref{eq:mean_model}) for the mean and probabilistic model according to {structure}, either $N(0, \sigma$), $N(0,\Sigma)$, conformal prediction (CP) or GARCH.
Estimated online for all but the LEAR-GARCH. \\
ODR-\{dep\} & (Univariate) distributional regression model (Equation \ref{eq:mean_model}, \ref{eq:model-scale}, \ref{eq:model-dof}) and independence or (sparse) Gaussian copula (dep=GC and dep=SPGC, Equation \ref{eq:gauss_copula}). Estimated online. \\
MODR-\{method\}-\{param\}-\{dep\} & Multivariate distributional regression models (Equation \ref{eq:mean_model}, \ref{eq:model-scale}, \ref{eq:model-dof}) estimated with method (OLS, LASSO) and the parameterization (see Section \ref{sec:covariance_parameterization}).
Assuming independence (dep=IND) or modelling the dependence structure (dep=$\Sigma$).
Estimated online. \\
\bottomrule
\end{tabular}
}
\caption{Overview and brief descriptions of the models included in the forecasting study.}\label{tab:models}
\end{table}
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{app:hyperparameters}.
We model the mean/location for all regression models, {\color{black} also for the LEAR models,} by \begin{equation}\label{eq:mean_model}
{\color{black}
\begin{aligned}
g_\mu(\mu_{t, h}) &= \beta_{\mu, 0, h}
+ \sum_{l=1}^{L=7} \beta_{\mu, l, h} y_{t-l, h}
+ \sum_{i \in \{0, ..., 23\} \setminus h} \beta_{\mu, 8+i, h} y_{t-1, i}
+ \sum_{w=1}^{W=6} \beta_{\mu, 30+w, h} \operatorname{WD}_{t, h} \\
&+ \beta_{\mu, 37, h} \operatorname{min}(\ensuremath{\bm{\mathrm{y}}}_{t-1})
+ \beta_{\mu, 38, h} \operatorname{max}(\ensuremath{\bm{\mathrm{y}}}_{t-1})
+ \beta_{\mu, 39, h} \operatorname{Q}_{10}(\ensuremath{\bm{\mathrm{y}}}_{t-1})
+ \beta_{\mu, 40, h} \operatorname{Q}_{90}(\ensuremath{\bm{\mathrm{y}}}_{t-1}) \\
&+ \beta_{\mu, 41, h} \operatorname{ResLoad}_{t, h}
+ \beta_{\mu, 42, h} \operatorname{EUA}_{t}
+ \beta_{\mu, 43, h} \operatorname{Gas}_{t}
+ \beta_{\mu, 44, h} \operatorname{Coal}_{t}
+ \beta_{\mu, 45, h} \operatorname{Oil}_{t}.
\end{aligned}}
\end{equation}
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 \begin{equation}
{\color{black} \begin{aligned}\label{eq:model-scale}
g_\theta(\theta_{t, h, h})
&= \beta_{\theta, 0, h}
+ \beta_{\theta, 1, h} \operatorname{SignedSquare}\left(\ensuremath{\bm{\mathrm{\Sigma}}}_{h,h}^{[t-1:t-7]}\right)^{-1}
+ \beta_{\theta, 2, h} \operatorname{ResLoad}_{t, h} \\
&+ \beta_{\theta, 4, h} \operatorname{EUA}_{t}
+ \beta_{\theta, 5, h} \operatorname{Gas}_{t}
+ \beta_{\theta, 6, h} \operatorname{Coal}_{t}
+ \beta_{\theta, 7, h} \operatorname{Oil}_{t},
\end{aligned}}
\end{equation}
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
\begin{equation}
{\color{black} \begin{aligned}\label{eq:model-scale-v}
g_v(v_{t, h, 0})
&= \beta_{v, 0, h}
+ \beta_{v, 1, h} \operatorname{SignedSquare}\left(\ensuremath{\bm{\mathrm{\Sigma}}}_{h,h}^{[t-1:t-7]}\right)
+ \beta_{v, 2, h} \overline{\operatorname{ResLoad}}_{t} \\
&+ \beta_{v, 3, h} \operatorname{EUA}_{t}
+ \beta_{v, 4, h} \operatorname{Gas}_{t}
+ \beta_{v, 5, h} \operatorname{Coal}_{t}
+ \beta_{v, 6, h} \operatorname{Oil}_{t},
\end{aligned}}
\end{equation}
\begin{equation}
{\color{black} \begin{aligned}\label{eq:model-v-2}
g_v(v_{t, h, 1}) = \sum_w^{W=6} \beta_{v, 14+w, h} \operatorname{WD}_{t, h}
\end{aligned}}
\end{equation}
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\begin{equation}
\begin{aligned}\label{eq:model-dof}
g_\nu(\nu_{t})
&= \beta_{\nu, 0}
+ \beta_{\nu, 1} \operatorname{mean}(\ensuremath{\bm{\mathrm{y}}}_{t-1})
+ \sum_w^{W=6} \beta_{\nu, 1+w, h} \operatorname{WD}_{t, h}
+ \beta_{\nu, 8} \overline{\operatorname{ResLoad}}_{t} \\
&+ \beta_{\nu, 9} \operatorname{EUA}_{t}
+ \beta_{\nu, 10} \operatorname{Gas}_{t}
+ \beta_{\nu, 11} \operatorname{Coal}_{t}
+ \beta_{\nu, 12} \operatorname{Oil}_{t}.
\end{aligned}
\end{equation} The univariate models are therefore a slight simplification compared to the models used in \cite{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: \begin{equation}
\ell(\ensuremath{\bm{\mathrm{u}}} \mid \ensuremath{\bm{\mathrm{\Sigma}}}) = \frac{1}{| \tilde{\ensuremath{\bm{\mathrm{\Sigma}}}}|^{1/2}} \exp\left(- \frac{1}{2} \ensuremath{\ensuremath{\bm{\mathrm{n}}}^\top} (\tilde{\ensuremath{\bm{\mathrm{\Sigma}}}}^{-1} - \ensuremath{\bm{\mathrm{I}}})\ensuremath{\bm{\mathrm{n}}}\right)\prod_{d=0}^Dp(y_d \mid \ensuremath{\bm{\mathrm{\theta_d}}})
\end{equation}
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 \citep[see][]{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 \begin{equation}\label{eq:gauss_copula}
\widehat{\ensuremath{\bm{\mathrm{\Sigma}}}}^{[t+1]} = \frac{t-1}{t} \widehat{\ensuremath{\bm{\mathrm{\Sigma}}}}^{[t]} + \frac{1}{t} \left(\ensuremath{\ensuremath{\bm{\mathrm{n}}}^{[t+1]}(\ensuremath{\bm{\mathrm{n}}}^{[t+1]})^\top} \right)
\end{equation}
where $\ensuremath{\bm{\mathrm{n}}}$ are the PIT-transformed in-sample values and the superscript $[t]$ denotes the observations available in the online learning \citep[see e.g.][]{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 \citep{friedman2008sparse}.
\subsection{Forecast Evaluation and Scoring Rules}
{\color{black} Forecast evaluation should follow the well-known principle of sharpness subject to calibration \citep{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 \citep{staszewskabystrova2011bootstrap, luetkepohl2015comparison}. They have been used in energy market forecasting by \cite{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: \begin{itemize}
\item The ES, DSS and LS are all able to reliably detect mis-specifications in the mean structure of the multivariate distribution \citep{marcotte2023regions}, while the VS (by design) is not sensitive. However, \cite{pinson2013discrimination,alexander2024evaluating} discuss the discrimination ability of the ES. While widely used for multivariate forecast evaluation, the ES has been shown to have low discrimination ability with respect to misspecified covariance structures, {\color{black} especially in dimensions $D > 10$ and \cite{alexander2024evaluating} recommend to use the VS in addition to the ES.}
\item In a similar vein, \cite{marcotte2023regions} show that the ES has low reliability, i.e. statistical power to discriminate between a correctly and incorrectly specified dependence model for multivariate forecasts compared to the DSS and the LS, with the VS being somewhat in the middle. Their analysis also reflects the role of the number of test samples (which, with more than 1800 out-of-sample days, is not a concern in our study) and the number of sample paths $M$. Interestingly, they find that the ES and the VS are complimentary, i.e. the VS is more reliable in cases where the ES is not and vice versa, and they find the reliability of the VS to increase non-monotonically with respect to the dimension $D$ and the number of samples $M$.
\item In contrast, \cite{ziel2019multivariate} find that the ES, in combination with the Diebold-Mariano test, is able to detect the correct model specification in a simulation study. Nevertheless, they still recommend the use of multiple scores.
\end{itemize}}
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 \cite{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 \begin{equation*}
\mathbb{P}_t\left(\ensuremath{\bm{\mathrm{l}}}_{t} \leq \ensuremath{\bm{\mathrm{y}}}_{t} \leq \ensuremath{\bm{\mathrm{u}}}_{t} \; \forall \; h \in \{0, ..., 23\} \right) = 1 - \alpha,
\end{equation*}
that is, the true price trajectory~$\ensuremath{\bm{\mathrm{y}}}_t$ is \emph{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 \cite{serafin2022trading, chen2025probabilistic}, we use the neighboring paths method described by \cite{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 \begin{equation}
\operatorname{MC}_{1-\alpha} = \frac{1}{T}\sum^T_{t=0}\boldsymbol{1}\left(\ensuremath{\bm{\mathrm{l}}}_{t} \leq \ensuremath{\bm{\mathrm{y}}}_{t} \leq \ensuremath{\bm{\mathrm{u}}}_{t} \; \forall \; h \in \{0, ..., 23\}\right) - (1-\alpha),
\end{equation}
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 \begin{equation}
\operatorname{JPBW}_{1-\alpha} = \frac{1}{TH}\sum^T_{t=0}\sum^H_{h=0}\left(\ensuremath{\bm{\mathrm{u}}}_{t,h} - \ensuremath{\bm{\mathrm{l}}}_{t,h}\right),
\end{equation} 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
\begin{equation}
\operatorname{RMSE} = \sqrt{\frac{1}{TH} \sum_{t=0}^T\sum_{h=0}^H\left( \ensuremath{\bm{\mathrm{y}}}_{t} - \hat{\ensuremath{\bm{\mathrm{\mu_t}}}}\right)^2}
\end{equation}
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 \begin{equation}
\operatorname{MAE} = \frac{1}{TH} \sum_{t=0}^T\sum_{h=0}^H\left| \ensuremath{\bm{\mathrm{y}}}_{t} - \operatorname{median}\left(\ensuremath{\bm{\mathrm{F_t}}}\right)\right|
\end{equation}
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 \cite{zamo2018estimation}:\begin{equation}
\operatorname{CRPS}_t = \frac{1}{M} \sum_{m=0}^M |\ensuremath{\bm{\mathrm{F}}}_{t,h,m} - y_{t.h}| + \frac{1}{M} \sum_{m=0}^M\ensuremath{\bm{\mathrm{F}}}_{t,h,m} + \frac{1}{M(M-1)}\sum_{m=0}^M m\ensuremath{\bm{\mathrm{F}}}_{t,h,m}.
\end{equation}
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 \citep[ES,][]{gneiting2007strictly} is defined as \begin{equation}
\operatorname{ES}_t = \frac{1}{M}\sum_{m=0}^M\lVert \ensuremath{\bm{\mathrm{y}}}_{t} - \ensuremath{\bm{\mathrm{F}}}_{t,m}\rVert_2^2 - \frac{1}{M^2}\sum_{i=0}^M\sum_{j=i+1}^M \lVert \ensuremath{\bm{\mathrm{F}}}_{t,i} - \ensuremath{\bm{\mathrm{F}}}_{t,j}\rVert_2^2.
\end{equation}
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 \begin{equation}
\operatorname{LS}_t = - \log\left(\mathcal{L}(\ensuremath{\bm{\mathrm{y}}}_t \mid \hat{\ensuremath{\bm{\mathrm{\theta}}}}^\mathcal{D}_t)\right),
\end{equation}
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 \citeyearpar[DSS,][]{dawid1999coherent} is defined as \begin{equation}
\operatorname{DSS}_t = \log\left(\det(\widehat{\ensuremath{\bm{\mathrm{\Sigma}}}}_F)\right) + (\ensuremath{\bm{\mathrm{y_t}}} - \ensuremath{\bm{\mathrm{\hat{\mu}_t}}})\widehat{\ensuremath{\bm{\mathrm{\Sigma}}}}_F^{-1}(\ensuremath{\bm{\mathrm{y_t}}} - \ensuremath{\bm{\mathrm{\hat{\mu}_t}}}),
\end{equation}
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 \citep[VS,][]{scheuerer2015variogram} is defined as \begin{equation}
\operatorname{VS}^p_t = \sum_{i=0}^H \sum_{j=0}^H \left(\frac{1}{M}\sum_{m=0}^M |\ensuremath{\bm{\mathrm{F}}}_{t,i,m} - \ensuremath{\bm{\mathrm{F}}}_{t,j,m}|^p - |y_{t, i} - y_{t,j}|^p\right)^2
\end{equation}
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 \texttt{scoringrules} \citep{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 \citep{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 \citep[ADF,][]{dickey1979distribution, cheung1995lag}.}
\section{Results}\label{sec:results}
{\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} \subsection{Forecasting the Time-varying Dependence Structure}}
{\color{black}
Figure~\ref{fig:insample_analysis} 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.
}
\begin{figure}[htb]
\centering
\includegraphics[width=\textwidth]{figures/evoluation_cov_and_corr.pdf}
\caption{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-6, which corresponds to the 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.}
\label{fig:insample_analysis}
\end{figure}
\FloatBarrier
{\color{black}
\subsection{Forecasting Accuracy}
}
{\color{black}
{\color{black}
We first examine the calibration of the joint prediction bands in Figure~\ref{fig:calibration_bands}. The figure shows miscoverage and mean JPB width for the $50\%, 55\%, ..., 95\%$ nominal coverage joint prediction bands. Table \ref{tab:results_scores} gives the results for the scoring rules for each model. Figure~\ref{fig:skill_scores} 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{sec:forecasting_study} and Table \ref{tab:models}: 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.
}
\begin{itemize}
\item {\color{black} \emph{Calibration:} In terms of the calibration of the joint prediction bands, the multivariate distributional regression models yield the best coverage across all widths. The MODR models for the CD and MCD parameterization yield prediction intervals closest to the nominal coverage, while the LRA-based models generally yield too small JPB. For the lower nominal coverage levels, the copula-based models yield too conservative (too wide) bands, while for higher nominal coverage levels, the LEAR-based models exhibit severe under-coverage. The JPB width of the multivariate distributional regression models is also lower than the univariate distributional regression models, however, the under-covering LEAR-based models yield the narrowest bands. }
\item \emph{Marginal Scores:} Here, the LEAR models yield the best performance in terms of RMSE, MAE and CRPS. These models yield sharp predictions, but at the cost of poor coverage of the joint prediction bands. This is also visible in the very small JPBW.
\item \emph{Multivariate Scores:} The Copula-based models yield the best performance in terms of the ES. The MCD-based multivariate distributional regression model yields the best overall performance in terms of the LS and DSS as well as a {\color{black} second-best performance in terms of the VS}.
\end{itemize}}e}}
\begin{figure}
\centering
\includegraphics[width=\textwidth]{figures/calibration_prediction_bands.png}
\caption{Calibration of the $50\%, 55\%, ..., 95\%$-joint prediction bands. A positive miscoverage indicates over-coverage, while a negative value indicates under-coverage. The second panel shows the mean width of the $50\%, 55\%, ..., 95\%$-joint prediction bands.}
\label{fig:calibration_bands}
\end{figure}
\begin{table}[h!]
\begin{center}
\centering \small Panel A: Full sample (2018-12-27 to 2023-12-31)
\resizebox{\textwidth}{!}{\begin{tabular}{R{6cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}}
\toprule
& RMSE & MAE & CRPS & $\text{VS}_{p=0.5}$ & $\text{VS}_{p=1}$ & ES & DSS & LS & $\text{MC}_{0.95}$ & $\text{JPBW}_{0.95}$ \\
\midrule
LEAR-N(0, $\sigma$) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 31.177 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 18.214 & {\cellcolor[HTML]{2A71B2}} \color[HTML]{F1F1F1} 13.677 & {\cellcolor[HTML]{D7634F}} \color[HTML]{F1F1F1} 4.709 & {\cellcolor[HTML]{2A71B2}} \color[HTML]{F1F1F1} 887.202 & {\cellcolor[HTML]{2267AC}} \color[HTML]{F1F1F1} 81.201 & {\cellcolor[HTML]{E8896C}} \color[HTML]{F1F1F1} 176.367 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 110.038 & {\cellcolor[HTML]{FBCEB7}} \color[HTML]{000000} -0.134 & {\cellcolor[HTML]{F9EDE5}} \color[HTML]{000000} 131.129 \\
LEAR-N(0, $\Sigma$) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 31.182 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 18.219 & {\cellcolor[HTML]{2A71B2}} \color[HTML]{F1F1F1} 13.678 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 3.396 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 654.953 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 78.603 & {\cellcolor[HTML]{195696}} \color[HTML]{F1F1F1} 124.345 & {\cellcolor[HTML]{9BC9E0}} \color[HTML]{000000} 84.057 & {\cellcolor[HTML]{F7B799}} \color[HTML]{000000} -0.161 & {\cellcolor[HTML]{F0F4F6}} \color[HTML]{000000} 118.914 \\
LEAR-CP & {\cellcolor[HTML]{063264}} \color[HTML]{F1F1F1} 31.268 & {\cellcolor[HTML]{0C3D73}} \color[HTML]{F1F1F1} 18.374 & {\cellcolor[HTML]{6BACD1}} \color[HTML]{F1F1F1} 14.402 & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.655 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 53.751 \\
LEAR-GARCH & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 31.181 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 18.218 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 12.992 & {\cellcolor[HTML]{F09C7B}} \color[HTML]{000000} 4.536 & {\cellcolor[HTML]{246AAE}} \color[HTML]{F1F1F1} 850.636 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 78.312 & {\cellcolor[HTML]{FAE8DE}} \color[HTML]{000000} 159.545 & {\cellcolor[HTML]{D35A4A}} \color[HTML]{F1F1F1} 101.790 & {\cellcolor[HTML]{BBDAEA}} \color[HTML]{000000} -0.040 & {\cellcolor[HTML]{FBE3D4}} \color[HTML]{000000} 140.459 \\
ODR-IND & {\cellcolor[HTML]{8AC0DB}} \color[HTML]{000000} 36.040 & {\cellcolor[HTML]{2C75B4}} \color[HTML]{F1F1F1} 19.104 & {\cellcolor[HTML]{3B88BE}} \color[HTML]{F1F1F1} 13.980 & {\cellcolor[HTML]{D55D4C}} \color[HTML]{F1F1F1} 4.722 & {\cellcolor[HTML]{4291C2}} \color[HTML]{F1F1F1} 1046.549 & {\cellcolor[HTML]{4F9BC7}} \color[HTML]{F1F1F1} 84.445 & {\cellcolor[HTML]{FBE5D8}} \color[HTML]{000000} 160.422 & {\cellcolor[HTML]{E6866A}} \color[HTML]{F1F1F1} 99.266 & {\cellcolor[HTML]{408FC1}} \color[HTML]{F1F1F1} 0.015 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 332.441 \\
ODR-GC & {\cellcolor[HTML]{3885BC}} \color[HTML]{F1F1F1} 33.901 & {\cellcolor[HTML]{114781}} \color[HTML]{F1F1F1} 18.511 & {\cellcolor[HTML]{1F63A8}} \color[HTML]{F1F1F1} 13.501 & {\cellcolor[HTML]{327CB7}} \color[HTML]{F1F1F1} 3.591 & {\cellcolor[HTML]{1B5A9C}} \color[HTML]{F1F1F1} 793.156 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 78.372 & {\cellcolor[HTML]{1B5A9C}} \color[HTML]{F1F1F1} 124.869 & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 77.934 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries -0.001 & {\cellcolor[HTML]{900D26}} \color[HTML]{F1F1F1} 294.100 \\
ODR-SPGC & {\cellcolor[HTML]{3A87BD}} \color[HTML]{F1F1F1} 33.964 & {\cellcolor[HTML]{124984}} \color[HTML]{F1F1F1} 18.522 & {\cellcolor[HTML]{1F63A8}} \color[HTML]{F1F1F1} 13.508 & {\cellcolor[HTML]{337EB8}} \color[HTML]{F1F1F1} 3.597 & {\cellcolor[HTML]{1B5A9C}} \color[HTML]{F1F1F1} 793.049 & {\cellcolor[HTML]{063264}} \color[HTML]{F1F1F1} 78.417 & {\cellcolor[HTML]{1A5899}} \color[HTML]{F1F1F1} 124.646 & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 77.964 & {\cellcolor[HTML]{134C87}} \color[HTML]{F1F1F1} 0.004 & {\cellcolor[HTML]{900D26}} \color[HTML]{F1F1F1} 294.163 \\
MODR-OLS-CD-IND & {\cellcolor[HTML]{E0ECF3}} \color[HTML]{000000} 38.898 & {\cellcolor[HTML]{8AC0DB}} \color[HTML]{000000} 20.207 & {\cellcolor[HTML]{D5E7F1}} \color[HTML]{000000} 15.410 & {\cellcolor[HTML]{87BEDA}} \color[HTML]{000000} 3.787 & {\cellcolor[HTML]{185493}} \color[HTML]{F1F1F1} 771.899 & {\cellcolor[HTML]{BDDBEA}} \color[HTML]{000000} 89.011 & {\cellcolor[HTML]{EF9979}} \color[HTML]{000000} 174.092 & {\cellcolor[HTML]{F9C4A9}} \color[HTML]{000000} 95.446 & {\cellcolor[HTML]{CAE1EE}} \color[HTML]{000000} -0.046 & {\cellcolor[HTML]{E98B6E}} \color[HTML]{F1F1F1} 191.574 \\
MODR-OLS-CD-$\Sigma$ & {\cellcolor[HTML]{DC6E57}} \color[HTML]{F1F1F1} 46.742 & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} 23.493 & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} 17.463 & {\cellcolor[HTML]{71B0D3}} \color[HTML]{F1F1F1} 3.747 & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 797.964 & {\cellcolor[HTML]{FAC8AF}} \color[HTML]{000000} 97.948 & {\cellcolor[HTML]{2267AC}} \color[HTML]{F1F1F1} 126.357 & {\cellcolor[HTML]{124984}} \color[HTML]{F1F1F1} 76.013 & {\cellcolor[HTML]{4C99C6}} \color[HTML]{F1F1F1} -0.017 & {\cellcolor[HTML]{AB162A}} \color[HTML]{F1F1F1} 273.946 \\
MODR-LASSO-CD-$\Sigma$ & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 52.791 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 26.483 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 20.094 & {\cellcolor[HTML]{E9F0F4}} \color[HTML]{000000} 4.070 & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 860.208 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 111.976 & {\cellcolor[HTML]{2870B1}} \color[HTML]{F1F1F1} 127.573 & {\cellcolor[HTML]{0E4179}} \color[HTML]{F1F1F1} 75.543 & {\cellcolor[HTML]{71B0D3}} \color[HTML]{F1F1F1} -0.023 & {\cellcolor[HTML]{C43B3C}} \color[HTML]{F1F1F1} 240.903 \\
MODR-OLS-MCD-IND & {\cellcolor[HTML]{E0ECF3}} \color[HTML]{000000} 38.937 & {\cellcolor[HTML]{87BEDA}} \color[HTML]{000000} 20.196 & {\cellcolor[HTML]{D2E6F0}} \color[HTML]{000000} 15.379 & {\cellcolor[HTML]{A7D0E4}} \color[HTML]{000000} 3.861 & {\cellcolor[HTML]{1B5A9C}} \color[HTML]{F1F1F1} 792.223 & {\cellcolor[HTML]{B8D8E9}} \color[HTML]{000000} 88.771 & {\cellcolor[HTML]{EA8E70}} \color[HTML]{F1F1F1} 175.635 & {\cellcolor[HTML]{F9C6AC}} \color[HTML]{000000} 95.274 & {\cellcolor[HTML]{A5CEE3}} \color[HTML]{000000} -0.034 & {\cellcolor[HTML]{E58368}} \color[HTML]{F1F1F1} 196.792 \\
MODR-OLS-MCD-$\Sigma$ & {\cellcolor[HTML]{DE735C}} \color[HTML]{F1F1F1} 46.450 & {\cellcolor[HTML]{F6AF8E}} \color[HTML]{000000} 23.376 & {\cellcolor[HTML]{F7B596}} \color[HTML]{000000} 17.255 & {\cellcolor[HTML]{337EB8}} \color[HTML]{F1F1F1} 3.596 & {\cellcolor[HTML]{144E8A}} \color[HTML]{F1F1F1} 747.954 & {\cellcolor[HTML]{FDDBC7}} \color[HTML]{000000} 96.739 & {\cellcolor[HTML]{0A3B70}} \color[HTML]{F1F1F1} 121.365 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 74.907 & {\cellcolor[HTML]{B1D5E7}} \color[HTML]{000000} -0.037 & {\cellcolor[HTML]{C94741}} \color[HTML]{F1F1F1} 234.422 \\
MODR-LASSO-MCD-$\Sigma$ & {\cellcolor[HTML]{AE172A}} \color[HTML]{F1F1F1} 50.034 & {\cellcolor[HTML]{BF3338}} \color[HTML]{F1F1F1} 25.078 & {\cellcolor[HTML]{BB2A34}} \color[HTML]{F1F1F1} 18.923 & {\cellcolor[HTML]{ACD2E5}} \color[HTML]{000000} 3.872 & {\cellcolor[HTML]{1D5FA2}} \color[HTML]{F1F1F1} 804.906 & {\cellcolor[HTML]{C53E3D}} \color[HTML]{F1F1F1} 105.786 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 119.963 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 74.646 & {\cellcolor[HTML]{F8F4F2}} \color[HTML]{000000} -0.080 & {\cellcolor[HTML]{E58368}} \color[HTML]{F1F1F1} 196.544 \\
MODR-OLS-LRA-IND & {\cellcolor[HTML]{E4EEF4}} \color[HTML]{000000} 39.197 & {\cellcolor[HTML]{93C6DE}} \color[HTML]{000000} 20.292 & {\cellcolor[HTML]{E6EFF4}} \color[HTML]{000000} 15.710 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 5.189 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} $>10,000$ & {\cellcolor[HTML]{F8F3F0}} \color[HTML]{000000} 93.714 & {\cellcolor[HTML]{E98B6E}} \color[HTML]{F1F1F1} 175.884 & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 95.602 & {\cellcolor[HTML]{D2E6F0}} \color[HTML]{000000} -0.049 & {\cellcolor[HTML]{E27B62}} \color[HTML]{F1F1F1} 201.897 \\
MODR-OLS-LRA-$\Sigma$ & {\cellcolor[HTML]{FBE4D6}} \color[HTML]{000000} 41.642 & {\cellcolor[HTML]{DDEBF2}} \color[HTML]{000000} 21.292 & {\cellcolor[HTML]{FBE4D6}} \color[HTML]{000000} 16.520 & {\cellcolor[HTML]{C7E0ED}} \color[HTML]{000000} 3.940 & {\cellcolor[HTML]{5FA5CD}} \color[HTML]{F1F1F1} 1161.734 & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} 94.398 & {\cellcolor[HTML]{EC9374}} \color[HTML]{F1F1F1} 174.901 & {\cellcolor[HTML]{FDD9C4}} \color[HTML]{000000} 93.993 & {\cellcolor[HTML]{F2F5F6}} \color[HTML]{000000} -0.072 & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} 176.183 \\
MODR-LASSO-LRA-$\Sigma$ & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} 39.664 & {\cellcolor[HTML]{A7D0E4}} \color[HTML]{000000} 20.550 & {\cellcolor[HTML]{FDDCC9}} \color[HTML]{000000} 16.699 & {\cellcolor[HTML]{6EAED2}} \color[HTML]{F1F1F1} 3.741 & {\cellcolor[HTML]{175290}} \color[HTML]{F1F1F1} 764.797 & {\cellcolor[HTML]{FAEAE1}} \color[HTML]{000000} 94.777 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 203.316 & {\cellcolor[HTML]{FCD5BF}} \color[HTML]{000000} 94.322 & {\cellcolor[HTML]{FBD0B9}} \color[HTML]{000000} -0.131 & {\cellcolor[HTML]{FCD7C2}} \color[HTML]{000000} 149.184 \\
\bottomrule
\end{tabular}
} \\
\vspace{0.25cm}
\centering \small Panel B: First sub-sample (2018-12-27 to 2020-12-31, $n = 736$)
\resizebox{\textwidth}{!}{\begin{tabular}{R{6cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}}
\toprule
& RMSE & MAE & CRPS & $\text{VS}_{p=0.5}$ & $\text{VS}_{p=1}$ & ES & DSS & LS & $\text{MC}_{0.95}$ & $\text{JPBW}_{0.95}$ \\
\midrule
LEAR-N(0, $\sigma$) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 7.568 & {\cellcolor[HTML]{F5F6F7}} \color[HTML]{000000} 4.790 & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} 3.588 & {\cellcolor[HTML]{BB2A34}} \color[HTML]{F1F1F1} 1.163 & {\cellcolor[HTML]{F9EBE3}} \color[HTML]{000000} 59.095 & {\cellcolor[HTML]{FCDFCF}} \color[HTML]{000000} 21.940 & {\cellcolor[HTML]{CC4C44}} \color[HTML]{F1F1F1} 120.417 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 82.117 & {\cellcolor[HTML]{E0ECF3}} \color[HTML]{000000} -0.063 & {\cellcolor[HTML]{FACAB1}} \color[HTML]{000000} 39.127 \\
LEAR-N(0, $\Sigma$) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 7.568 & {\cellcolor[HTML]{F5F6F7}} \color[HTML]{000000} 4.791 & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} 3.588 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 0.912 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 52.607 & {\cellcolor[HTML]{B1D5E7}} \color[HTML]{000000} 21.368 & {\cellcolor[HTML]{4F9BC7}} \color[HTML]{F1F1F1} 85.142 & {\cellcolor[HTML]{D4E6F1}} \color[HTML]{000000} 64.467 & {\cellcolor[HTML]{F7F6F6}} \color[HTML]{000000} -0.080 & {\cellcolor[HTML]{FBE4D6}} \color[HTML]{000000} 36.139 \\
LEAR-CP & {\cellcolor[HTML]{4291C2}} \color[HTML]{F1F1F1} 7.751 & {\cellcolor[HTML]{D05548}} \color[HTML]{F1F1F1} 4.959 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 3.831 & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.525 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 18.469 \\
LEAR-GARCH & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 7.567 & {\cellcolor[HTML]{F5F6F7}} \color[HTML]{000000} 4.791 & {\cellcolor[HTML]{F7F5F4}} \color[HTML]{000000} 3.601 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 1.258 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 65.391 & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 22.099 & {\cellcolor[HTML]{D6604D}} \color[HTML]{F1F1F1} 118.522 & {\cellcolor[HTML]{760521}} \color[HTML]{F1F1F1} 81.369 & {\cellcolor[HTML]{C5DFEC}} \color[HTML]{000000} -0.049 & {\cellcolor[HTML]{F6AF8E}} \color[HTML]{000000} 42.002 \\
ODR-IND & {\cellcolor[HTML]{90C4DD}} \color[HTML]{000000} 7.849 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 4.552 & {\cellcolor[HTML]{0A3B70}} \color[HTML]{F1F1F1} 3.399 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.994 & {\cellcolor[HTML]{98C8E0}} \color[HTML]{000000} 56.194 & {\cellcolor[HTML]{4393C3}} \color[HTML]{F1F1F1} 21.024 & {\cellcolor[HTML]{F09C7B}} \color[HTML]{000000} 112.649 & {\cellcolor[HTML]{E37E64}} \color[HTML]{F1F1F1} 74.195 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} -0.006 & {\cellcolor[HTML]{D05548}} \color[HTML]{F1F1F1} 50.476 \\
ODR-GC & {\cellcolor[HTML]{7EB8D7}} \color[HTML]{000000} 7.828 & {\cellcolor[HTML]{063264}} \color[HTML]{F1F1F1} 4.551 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 3.394 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.877 & {\cellcolor[HTML]{1F63A8}} \color[HTML]{F1F1F1} 53.687 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 20.579 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 76.162 & {\cellcolor[HTML]{1A5899}} \color[HTML]{F1F1F1} 56.599 & {\cellcolor[HTML]{5CA3CB}} \color[HTML]{F1F1F1} -0.025 & {\cellcolor[HTML]{F3A481}} \color[HTML]{000000} 43.055 \\
ODR-SPGC & {\cellcolor[HTML]{78B4D5}} \color[HTML]{000000} 7.819 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 4.548 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 3.390 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.877 & {\cellcolor[HTML]{1F63A8}} \color[HTML]{F1F1F1} 53.706 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 20.557 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 75.768 & {\cellcolor[HTML]{195696}} \color[HTML]{F1F1F1} 56.561 & {\cellcolor[HTML]{4291C2}} \color[HTML]{F1F1F1} -0.021 & {\cellcolor[HTML]{F3A481}} \color[HTML]{000000} 43.102 \\
MODR-OLS-CD-IND & {\cellcolor[HTML]{F0F4F6}} \color[HTML]{000000} 8.034 & {\cellcolor[HTML]{65A9CF}} \color[HTML]{F1F1F1} 4.664 & {\cellcolor[HTML]{D2E6F0}} \color[HTML]{000000} 3.557 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.952 & {\cellcolor[HTML]{A9D1E5}} \color[HTML]{000000} 56.556 & {\cellcolor[HTML]{F8F2EF}} \color[HTML]{000000} 21.777 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 133.695 & {\cellcolor[HTML]{F7B799}} \color[HTML]{000000} 71.549 & {\cellcolor[HTML]{DEEBF2}} \color[HTML]{000000} -0.061 & {\cellcolor[HTML]{EB9172}} \color[HTML]{F1F1F1} 44.810 \\
MODR-OLS-CD-$\Sigma$ & {\cellcolor[HTML]{F7F6F6}} \color[HTML]{000000} 8.050 & {\cellcolor[HTML]{ACD2E5}} \color[HTML]{000000} 4.713 & {\cellcolor[HTML]{B3D6E8}} \color[HTML]{000000} 3.536 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.916 & {\cellcolor[HTML]{68ABD0}} \color[HTML]{F1F1F1} 55.459 & {\cellcolor[HTML]{B1D5E7}} \color[HTML]{000000} 21.373 & {\cellcolor[HTML]{4393C3}} \color[HTML]{F1F1F1} 84.530 & {\cellcolor[HTML]{124984}} \color[HTML]{F1F1F1} 56.039 & {\cellcolor[HTML]{B8D8E9}} \color[HTML]{000000} -0.045 & {\cellcolor[HTML]{D7634F}} \color[HTML]{F1F1F1} 49.244 \\
MODR-LASSO-CD-$\Sigma$ & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 8.579 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 5.062 & {\cellcolor[HTML]{6D0220}} \color[HTML]{F1F1F1} 3.827 & {\cellcolor[HTML]{65A9CF}} \color[HTML]{F1F1F1} 1.012 & {\cellcolor[HTML]{F19E7D}} \color[HTML]{000000} 61.332 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 23.007 & {\cellcolor[HTML]{81BAD8}} \color[HTML]{000000} 88.095 & {\cellcolor[HTML]{114781}} \color[HTML]{F1F1F1} 55.962 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 0.005 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 66.428 \\
MODR-OLS-MCD-IND & {\cellcolor[HTML]{F7F5F4}} \color[HTML]{000000} 8.056 & {\cellcolor[HTML]{5CA3CB}} \color[HTML]{F1F1F1} 4.659 & {\cellcolor[HTML]{8AC0DB}} \color[HTML]{000000} 3.509 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.961 & {\cellcolor[HTML]{A2CDE3}} \color[HTML]{000000} 56.386 & {\cellcolor[HTML]{DBEAF2}} \color[HTML]{000000} 21.557 & {\cellcolor[HTML]{C43B3C}} \color[HTML]{F1F1F1} 122.110 & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 70.944 & {\cellcolor[HTML]{3480B9}} \color[HTML]{F1F1F1} -0.017 & {\cellcolor[HTML]{D55D4C}} \color[HTML]{F1F1F1} 49.627 \\
MODR-OLS-MCD-$\Sigma$ & {\cellcolor[HTML]{F8F4F2}} \color[HTML]{000000} 8.059 & {\cellcolor[HTML]{9BC9E0}} \color[HTML]{000000} 4.699 & {\cellcolor[HTML]{84BCD9}} \color[HTML]{000000} 3.505 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.900 & {\cellcolor[HTML]{3E8CBF}} \color[HTML]{F1F1F1} 54.724 & {\cellcolor[HTML]{98C8E0}} \color[HTML]{000000} 21.277 & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 80.581 & {\cellcolor[HTML]{08366A}} \color[HTML]{F1F1F1} 55.291 & {\cellcolor[HTML]{98C8E0}} \color[HTML]{000000} -0.036 & {\cellcolor[HTML]{E6866A}} \color[HTML]{F1F1F1} 45.742 \\
MODR-LASSO-MCD-$\Sigma$ & {\cellcolor[HTML]{CF5246}} \color[HTML]{F1F1F1} 8.380 & {\cellcolor[HTML]{F19E7D}} \color[HTML]{000000} 4.903 & {\cellcolor[HTML]{F6B191}} \color[HTML]{000000} 3.678 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.965 & {\cellcolor[HTML]{FBCCB4}} \color[HTML]{000000} 60.209 & {\cellcolor[HTML]{EC9374}} \color[HTML]{F1F1F1} 22.292 & {\cellcolor[HTML]{175290}} \color[HTML]{F1F1F1} 78.513 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 54.989 & {\cellcolor[HTML]{DEEBF2}} \color[HTML]{000000} -0.061 & {\cellcolor[HTML]{FDDBC7}} \color[HTML]{000000} 37.680 \\
MODR-OLS-LRA-IND & {\cellcolor[HTML]{F6F7F7}} \color[HTML]{000000} 8.047 & {\cellcolor[HTML]{62A7CE}} \color[HTML]{F1F1F1} 4.662 & {\cellcolor[HTML]{D8E9F1}} \color[HTML]{000000} 3.564 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.941 & {\cellcolor[HTML]{C0DCEB}} \color[HTML]{000000} 56.963 & {\cellcolor[HTML]{F9EDE5}} \color[HTML]{000000} 21.830 & {\cellcolor[HTML]{AB162A}} \color[HTML]{F1F1F1} 126.458 & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} 71.156 & {\cellcolor[HTML]{ACD2E5}} \color[HTML]{000000} -0.041 & {\cellcolor[HTML]{EA8E70}} \color[HTML]{F1F1F1} 44.971 \\
MODR-OLS-LRA-$\Sigma$ & {\cellcolor[HTML]{F7F5F4}} \color[HTML]{000000} 8.057 & {\cellcolor[HTML]{7EB8D7}} \color[HTML]{000000} 4.681 & {\cellcolor[HTML]{DAE9F2}} \color[HTML]{000000} 3.566 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.952 & {\cellcolor[HTML]{B1D5E7}} \color[HTML]{000000} 56.712 & {\cellcolor[HTML]{F9EBE3}} \color[HTML]{000000} 21.840 & {\cellcolor[HTML]{CC4C44}} \color[HTML]{F1F1F1} 120.469 & {\cellcolor[HTML]{FDDCC9}} \color[HTML]{000000} 69.491 & {\cellcolor[HTML]{90C4DD}} \color[HTML]{000000} -0.034 & {\cellcolor[HTML]{E58368}} \color[HTML]{F1F1F1} 45.885 \\
MODR-LASSO-LRA-$\Sigma$ & {\cellcolor[HTML]{F9EFE9}} \color[HTML]{000000} 8.078 & {\cellcolor[HTML]{4997C5}} \color[HTML]{F1F1F1} 4.649 & {\cellcolor[HTML]{9BC9E0}} \color[HTML]{000000} 3.519 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.993 & {\cellcolor[HTML]{E7F0F4}} \color[HTML]{000000} 58.064 & {\cellcolor[HTML]{E1EDF3}} \color[HTML]{000000} 21.594 & {\cellcolor[HTML]{CF5246}} \color[HTML]{F1F1F1} 119.889 & {\cellcolor[HTML]{FDD9C4}} \color[HTML]{000000} 69.719 & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} -0.006 & {\cellcolor[HTML]{C13639}} \color[HTML]{F1F1F1} 53.847 \\
\bottomrule
\end{tabular}
} \\
\vspace{0.25cm}
\centering \small Panel C: Second sub-sample (2021-01-01 to 2023-12-31, $n = 1095$)
\resizebox{\textwidth}{!}{\begin{tabular}{R{6cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}L{1.9cm}}
\toprule
& RMSE & MAE & CRPS & $\text{VS}_{p=0.5}$ & $\text{VS}_{p=1}$ & ES & DSS & LS & $\text{MC}_{0.95}$ & $\text{JPBW}_{0.95}$ \\
\midrule
LEAR-N(0, $\sigma$) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 39.835 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 27.237 & {\cellcolor[HTML]{2B73B3}} \color[HTML]{F1F1F1} 20.459 & {\cellcolor[HTML]{E48066}} \color[HTML]{F1F1F1} 7.093 & {\cellcolor[HTML]{2A71B2}} \color[HTML]{F1F1F1} 1443.811 & {\cellcolor[HTML]{246AAE}} \color[HTML]{F1F1F1} 121.033 & {\cellcolor[HTML]{F7B596}} \color[HTML]{000000} 213.973 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 128.804 & {\cellcolor[HTML]{F6B394}} \color[HTML]{000000} -0.182 & {\cellcolor[HTML]{F9F0EB}} \color[HTML]{000000} 192.967 \\
LEAR-N(0, $\Sigma$) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 39.841 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 27.245 & {\cellcolor[HTML]{2B73B3}} \color[HTML]{F1F1F1} 20.460 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 5.065 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 1059.817 & {\cellcolor[HTML]{0A3B70}} \color[HTML]{F1F1F1} 117.073 & {\cellcolor[HTML]{0F437B}} \color[HTML]{F1F1F1} 150.695 & {\cellcolor[HTML]{7BB6D6}} \color[HTML]{000000} 97.225 & {\cellcolor[HTML]{F19E7D}} \color[HTML]{000000} -0.215 & {\cellcolor[HTML]{ECF2F5}} \color[HTML]{000000} 174.550 \\
LEAR-CP & {\cellcolor[HTML]{063264}} \color[HTML]{F1F1F1} 39.931 & {\cellcolor[HTML]{08366A}} \color[HTML]{F1F1F1} 27.391 & {\cellcolor[HTML]{62A7CE}} \color[HTML]{F1F1F1} 21.507 & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} \color{white} & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.742 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 77.465 \\
LEAR-GARCH & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 39.840 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 27.243 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 19.304 & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 6.739 & {\cellcolor[HTML]{246AAE}} \color[HTML]{F1F1F1} 1378.435 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 116.095 & {\cellcolor[HTML]{DDEBF2}} \color[HTML]{000000} 187.118 & {\cellcolor[HTML]{EC9374}} \color[HTML]{F1F1F1} 115.516 & {\cellcolor[HTML]{9DCBE1}} \color[HTML]{000000} -0.034 & {\cellcolor[HTML]{FBE6DA}} \color[HTML]{000000} 206.636 \\
ODR-IND & {\cellcolor[HTML]{8AC0DB}} \color[HTML]{000000} 46.157 & {\cellcolor[HTML]{327CB7}} \color[HTML]{F1F1F1} 28.885 & {\cellcolor[HTML]{4291C2}} \color[HTML]{F1F1F1} 21.092 & {\cellcolor[HTML]{D86551}} \color[HTML]{F1F1F1} 7.228 & {\cellcolor[HTML]{4291C2}} \color[HTML]{F1F1F1} 1712.213 & {\cellcolor[HTML]{5FA5CD}} \color[HTML]{F1F1F1} 127.073 & {\cellcolor[HTML]{F0F4F6}} \color[HTML]{000000} 192.533 & {\cellcolor[HTML]{E98B6E}} \color[HTML]{F1F1F1} 116.117 & {\cellcolor[HTML]{8AC0DB}} \color[HTML]{000000} 0.029 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 521.962 \\
ODR-GC & {\cellcolor[HTML]{3885BC}} \color[HTML]{F1F1F1} 43.365 & {\cellcolor[HTML]{15508D}} \color[HTML]{F1F1F1} 27.894 & {\cellcolor[HTML]{246AAE}} \color[HTML]{F1F1F1} 20.295 & {\cellcolor[HTML]{3681BA}} \color[HTML]{F1F1F1} 5.415 & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 1290.187 & {\cellcolor[HTML]{0C3D73}} \color[HTML]{F1F1F1} 117.218 & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 157.607 & {\cellcolor[HTML]{2C75B4}} \color[HTML]{F1F1F1} 92.275 & {\cellcolor[HTML]{3B88BE}} \color[HTML]{F1F1F1} 0.014 & {\cellcolor[HTML]{8D0C25}} \color[HTML]{F1F1F1} 462.839 \\
ODR-SPGC & {\cellcolor[HTML]{3A87BD}} \color[HTML]{F1F1F1} 43.450 & {\cellcolor[HTML]{175290}} \color[HTML]{F1F1F1} 27.915 & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 20.309 & {\cellcolor[HTML]{3783BB}} \color[HTML]{F1F1F1} 5.425 & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 1289.995 & {\cellcolor[HTML]{0D3F76}} \color[HTML]{F1F1F1} 117.308 & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 157.499 & {\cellcolor[HTML]{2C75B4}} \color[HTML]{F1F1F1} 92.349 & {\cellcolor[HTML]{5CA3CB}} \color[HTML]{F1F1F1} 0.021 & {\cellcolor[HTML]{8D0C25}} \color[HTML]{F1F1F1} 462.912 \\
MODR-OLS-CD-IND & {\cellcolor[HTML]{E0ECF3}} \color[HTML]{000000} 49.867 & {\cellcolor[HTML]{93C6DE}} \color[HTML]{000000} 30.654 & {\cellcolor[HTML]{D8E9F1}} \color[HTML]{000000} 23.377 & {\cellcolor[HTML]{7EB8D7}} \color[HTML]{000000} 5.693 & {\cellcolor[HTML]{185493}} \color[HTML]{F1F1F1} 1252.714 & {\cellcolor[HTML]{C5DFEC}} \color[HTML]{000000} 134.201 & {\cellcolor[HTML]{FBE6DA}} \color[HTML]{000000} 201.244 & {\cellcolor[HTML]{FACAB1}} \color[HTML]{000000} 111.508 & {\cellcolor[HTML]{A2CDE3}} \color[HTML]{000000} -0.035 & {\cellcolor[HTML]{EB9172}} \color[HTML]{F1F1F1} 290.222 \\
MODR-OLS-CD-$\Sigma$ & {\cellcolor[HTML]{DB6B55}} \color[HTML]{F1F1F1} 60.081 & {\cellcolor[HTML]{F2A17F}} \color[HTML]{000000} 36.115 & {\cellcolor[HTML]{F19E7D}} \color[HTML]{000000} 26.825 & {\cellcolor[HTML]{71B0D3}} \color[HTML]{F1F1F1} 5.650 & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 1297.036 & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} 149.418 & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 154.471 & {\cellcolor[HTML]{124984}} \color[HTML]{F1F1F1} 89.439 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 0.002 & {\cellcolor[HTML]{AB162A}} \color[HTML]{F1F1F1} 424.979 \\
MODR-LASSO-CD-$\Sigma$ & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 67.902 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 40.882 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 31.029 & {\cellcolor[HTML]{E1EDF3}} \color[HTML]{000000} 6.126 & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 1397.169 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 171.775 & {\cellcolor[HTML]{1A5899}} \color[HTML]{F1F1F1} 154.108 & {\cellcolor[HTML]{0C3D73}} \color[HTML]{F1F1F1} 88.704 & {\cellcolor[HTML]{B6D7E8}} \color[HTML]{000000} -0.041 & {\cellcolor[HTML]{CB4942}} \color[HTML]{F1F1F1} 358.176 \\
MODR-OLS-MCD-IND & {\cellcolor[HTML]{E0ECF3}} \color[HTML]{000000} 49.915 & {\cellcolor[HTML]{93C6DE}} \color[HTML]{000000} 30.639 & {\cellcolor[HTML]{D8E9F1}} \color[HTML]{000000} 23.357 & {\cellcolor[HTML]{A0CCE2}} \color[HTML]{000000} 5.810 & {\cellcolor[HTML]{1B5A9C}} \color[HTML]{F1F1F1} 1286.812 & {\cellcolor[HTML]{C0DCEB}} \color[HTML]{000000} 133.948 & {\cellcolor[HTML]{F8BFA4}} \color[HTML]{000000} 211.612 & {\cellcolor[HTML]{FAC8AF}} \color[HTML]{000000} 111.628 & {\cellcolor[HTML]{C0DCEB}} \color[HTML]{000000} -0.045 & {\cellcolor[HTML]{E98B6E}} \color[HTML]{F1F1F1} 295.708 \\
MODR-OLS-MCD-$\Sigma$ & {\cellcolor[HTML]{DE735C}} \color[HTML]{F1F1F1} 59.700 & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} 35.929 & {\cellcolor[HTML]{F6AF8E}} \color[HTML]{000000} 26.497 & {\cellcolor[HTML]{3480B9}} \color[HTML]{F1F1F1} 5.407 & {\cellcolor[HTML]{144E8A}} \color[HTML]{F1F1F1} 1213.906 & {\cellcolor[HTML]{FBD0B9}} \color[HTML]{000000} 147.461 & {\cellcolor[HTML]{08366A}} \color[HTML]{F1F1F1} 148.778 & {\cellcolor[HTML]{063264}} \color[HTML]{F1F1F1} 88.092 & {\cellcolor[HTML]{A9D1E5}} \color[HTML]{000000} -0.038 & {\cellcolor[HTML]{C94741}} \color[HTML]{F1F1F1} 361.243 \\
MODR-LASSO-MCD-$\Sigma$ & {\cellcolor[HTML]{AE172A}} \color[HTML]{F1F1F1} 64.334 & {\cellcolor[HTML]{BD2D35}} \color[HTML]{F1F1F1} 38.639 & {\cellcolor[HTML]{BA2832}} \color[HTML]{F1F1F1} 29.170 & {\cellcolor[HTML]{A5CEE3}} \color[HTML]{000000} 5.827 & {\cellcolor[HTML]{1D5FA2}} \color[HTML]{F1F1F1} 1305.452 & {\cellcolor[HTML]{C2383A}} \color[HTML]{F1F1F1} 161.906 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 147.824 & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} \bfseries 87.859 & {\cellcolor[HTML]{F8F1ED}} \color[HTML]{000000} -0.092 & {\cellcolor[HTML]{E48066}} \color[HTML]{F1F1F1} 303.324 \\
MODR-OLS-LRA-IND & {\cellcolor[HTML]{E4EEF4}} \color[HTML]{000000} 50.255 & {\cellcolor[HTML]{9BC9E0}} \color[HTML]{000000} 30.798 & {\cellcolor[HTML]{EAF1F5}} \color[HTML]{000000} 23.874 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 8.044 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} $>10,000$ & {\cellcolor[HTML]{F9F0EB}} \color[HTML]{000000} 142.030 & {\cellcolor[HTML]{FACAB1}} \color[HTML]{000000} 209.106 & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 112.033 & {\cellcolor[HTML]{D4E6F1}} \color[HTML]{000000} -0.054 & {\cellcolor[HTML]{E37E64}} \color[HTML]{F1F1F1} 307.374 \\
MODR-OLS-LRA-$\Sigma$ & {\cellcolor[HTML]{FBE4D6}} \color[HTML]{000000} 53.441 & {\cellcolor[HTML]{E3EDF3}} \color[HTML]{000000} 32.456 & {\cellcolor[HTML]{FCE0D0}} \color[HTML]{000000} 25.226 & {\cellcolor[HTML]{C2DDEC}} \color[HTML]{000000} 5.949 & {\cellcolor[HTML]{62A7CE}} \color[HTML]{F1F1F1} 1904.470 & {\cellcolor[HTML]{FAEAE1}} \color[HTML]{000000} 143.167 & {\cellcolor[HTML]{F8BFA4}} \color[HTML]{000000} 211.488 & {\cellcolor[HTML]{FCD7C2}} \color[HTML]{000000} 110.463 & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} -0.098 & {\cellcolor[HTML]{F6AF8E}} \color[HTML]{000000} 263.763 \\
MODR-LASSO-LRA-$\Sigma$ & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} 50.860 & {\cellcolor[HTML]{B1D5E7}} \color[HTML]{000000} 31.239 & {\cellcolor[HTML]{FCD7C2}} \color[HTML]{000000} 25.558 & {\cellcolor[HTML]{5CA3CB}} \color[HTML]{F1F1F1} 5.588 & {\cellcolor[HTML]{175290}} \color[HTML]{F1F1F1} 1239.825 & {\cellcolor[HTML]{FBE6DA}} \color[HTML]{000000} 143.967 & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 259.391 & {\cellcolor[HTML]{FBD0B9}} \color[HTML]{000000} 110.858 & {\cellcolor[HTML]{F19E7D}} \color[HTML]{000000} -0.215 & {\cellcolor[HTML]{FCE2D2}} \color[HTML]{000000} 213.264 \\
\bottomrule
\end{tabular}
}
\caption{Scoring Rules. The best score in each column is marked \textbf{bold}. Note that the $LEAR-\mathcal{N}(0, \sigma)$, the ODR-IND and the MODR-IND models do not model the dependence structure.}
\end{center} \label{tab:results_scores}
\end{table}
The $p$-values of the DM-test {\color{black} (see Figure \ref{fig:dm_tests} in \ref{app:figures})} 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 \cite{buchweitz2025asymmetric} and the different behaviors observed by \cite{pinson2013discrimination, marcotte2023regions} and \cite{alexander2024evaluating}.}
{\color{black}
Looking at Figure \ref{fig:skill_scores}, 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.
}
\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{figures/skill_scores_timeseries.pdf}
\caption{Rolling 182-day average of skill scores. The skill score is defined as $\operatorname{SS} = 1 - \operatorname{LS}_\text{model} / \operatorname{LS}_\text{baseline}$, where the baseline is the LEAR+N(0,$\Sigma$) model. A positive skill score indicates an improvement over the baseline. The secondary $y$-axis shows the daily average day-ahead price.}
\label{fig:skill_scores}
\end{figure}
{\color{black}
A common theme for distributional regression models is the trade-off between accuracy in the point predictions and the distributional error metrics. \cite{marcjasz2023distributional} and \cite{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 \citep[see e.g.][{\color{black} who focus on the DSS and LS as evaluation metrics}]{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{sec:overfitting}. 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: \begin{align*}
\operatorname{Spike} &=
\{t \mid \min(\ensuremath{\bm{\mathrm{y}}}_t) < \operatorname{Q}_{0.05}(\min(\ensuremath{\bm{\mathrm{y}}}_t)) \lor \max(\ensuremath{\bm{\mathrm{y}}}_t) > \operatorname{Q}_{0.05}(\max(\ensuremath{\bm{\mathrm{y}}}_t))\}, \quad \text{and}\\
\operatorname{Spread} &= \{t \mid \max(\ensuremath{\bm{\mathrm{y}}}_t) - \min(\ensuremath{\bm{\mathrm{y}}}_t) > \operatorname{Q}_{0.90}(\max(\ensuremath{\bm{\mathrm{y}}}_t) - \min(\ensuremath{\bm{\mathrm{y}}}_t))\}.
\end{align*}
The results for the CRPS, LS and the miscoverage are given in Table \ref{tab:results_spike}. 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{tab:results_spike_app} and mirror roughly the results for the CRPS and LS. Figures \ref{fig:calibration_bands_spikes} and \ref{fig:calibration_bands_spreads} 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.
}}
\begin{table}[htb!]
\begin{center}
\resizebox{\textwidth}{!}{\begin{tabular}{lrrrrrrrrrrrrrrrrrrrrrrrr}
\toprule
& \multicolumn{4}{c}{CRPS} & \multicolumn{4}{c}{LS} & \multicolumn{4}{c}{$\text{MC}_{0.95}$} \\
& \multicolumn{2}{c}{Spread} & \multicolumn{2}{c}{Spike} & \multicolumn{2}{c}{Spread} & \multicolumn{2}{c}{Spike} & \multicolumn{2}{c}{Spread} & \multicolumn{2}{c}{Spike} \\
& Low & High & Low & High & Low & High & Low & High & Low & High & Low & High \\
\midrule
LEAR-N(0, $\sigma$) & {\cellcolor[HTML]{C5DFEC}} \color[HTML]{000000} 10.99 (10) & {\cellcolor[HTML]{1D5FA2}} \color[HTML]{F1F1F1} 37.85 (02) & {\cellcolor[HTML]{7EB8D7}} \color[HTML]{000000} 11.96 (07) & {\cellcolor[HTML]{195696}} \color[HTML]{F1F1F1} 29.02 (03) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 103.42 (15) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 169.63 (15) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 104.74 (15) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 157.46 (15) & {\cellcolor[HTML]{EA8E70}} \color[HTML]{F1F1F1} -0.11 (14) & {\cellcolor[HTML]{D6604D}} \color[HTML]{F1F1F1} -0.37 (13) & {\cellcolor[HTML]{EE9677}} \color[HTML]{000000} -0.10 (14) & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} -0.45 (12) \\
LEAR-N(0, $\Sigma$) & {\cellcolor[HTML]{C5DFEC}} \color[HTML]{000000} 10.99 (09) & {\cellcolor[HTML]{1D5FA2}} \color[HTML]{F1F1F1} 37.85 (03) & {\cellcolor[HTML]{7EB8D7}} \color[HTML]{000000} 11.96 (08) & {\cellcolor[HTML]{195696}} \color[HTML]{F1F1F1} 29.02 (02) & {\cellcolor[HTML]{9DCBE1}} \color[HTML]{000000} 79.69 (07) & {\cellcolor[HTML]{6EAED2}} \color[HTML]{F1F1F1} 123.39 (07) & {\cellcolor[HTML]{87BEDA}} \color[HTML]{000000} 80.07 (07) & {\cellcolor[HTML]{D4E6F1}} \color[HTML]{000000} 119.74 (07) & {\cellcolor[HTML]{E27B62}} \color[HTML]{F1F1F1} -0.13 (15) & {\cellcolor[HTML]{C6413E}} \color[HTML]{F1F1F1} -0.43 (14) & {\cellcolor[HTML]{E48066}} \color[HTML]{F1F1F1} -0.12 (15) & {\cellcolor[HTML]{E27B62}} \color[HTML]{F1F1F1} -0.52 (14) \\
LEAR-CP & {\cellcolor[HTML]{FBD0B9}} \color[HTML]{000000} 11.75 (13) & {\cellcolor[HTML]{2065AB}} \color[HTML]{F1F1F1} 38.27 (04) & {\cellcolor[HTML]{E6EFF4}} \color[HTML]{000000} 12.64 (12) & {\cellcolor[HTML]{266CAF}} \color[HTML]{F1F1F1} 30.16 (04) & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} (17) & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} (17) & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} (17) & {\cellcolor[HTML]{FFFFFF}} \color[HTML]{000000} (17) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.64 (16) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.79 (16) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.63 (16) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} -0.86 (16) \\
LEAR-GARCH & {\cellcolor[HTML]{59A1CA}} \color[HTML]{F1F1F1} 10.54 (06) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 35.11 (01) & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 11.40 (04) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 27.27 (01) & {\cellcolor[HTML]{BF3338}} \color[HTML]{F1F1F1} 98.06 (14) & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} 135.37 (08) & {\cellcolor[HTML]{D7634F}} \color[HTML]{F1F1F1} 97.02 (14) & {\cellcolor[HTML]{D05548}} \color[HTML]{F1F1F1} 144.52 (14) & {\cellcolor[HTML]{FAE8DE}} \color[HTML]{000000} -0.03 (08) & {\cellcolor[HTML]{EAF1F5}} \color[HTML]{000000} -0.10 (05) & {\cellcolor[HTML]{D7E8F1}} \color[HTML]{000000} -0.01 (05) & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} -0.28 (07) \\
ODR-IND & {\cellcolor[HTML]{144E8A}} \color[HTML]{F1F1F1} 10.03 (03) & {\cellcolor[HTML]{D1E5F0}} \color[HTML]{000000} 49.54 (07) & {\cellcolor[HTML]{124984}} \color[HTML]{F1F1F1} 11.16 (03) & {\cellcolor[HTML]{D4E6F1}} \color[HTML]{000000} 39.18 (07) & {\cellcolor[HTML]{E27B62}} \color[HTML]{F1F1F1} 94.39 (13) & {\cellcolor[HTML]{FDD9C4}} \color[HTML]{000000} 143.14 (14) & {\cellcolor[HTML]{E37E64}} \color[HTML]{F1F1F1} 95.57 (13) & {\cellcolor[HTML]{FBCEB7}} \color[HTML]{000000} 132.34 (13) & {\cellcolor[HTML]{D5E7F1}} \color[HTML]{000000} 0.01 (04) & {\cellcolor[HTML]{3A87BD}} \color[HTML]{F1F1F1} 0.03 (03) & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} 0.03 (10) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} -0.10 (01) \\
ODR-GC & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 9.87 (01) & {\cellcolor[HTML]{9DCBE1}} \color[HTML]{000000} 46.18 (05) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 10.99 (01) & {\cellcolor[HTML]{98C8E0}} \color[HTML]{000000} 35.96 (05) & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 73.87 (05) & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 114.50 (05) & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 75.17 (05) & {\cellcolor[HTML]{195696}} \color[HTML]{F1F1F1} 102.66 (06) & {\cellcolor[HTML]{4997C5}} \color[HTML]{F1F1F1} -0.00 (03) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.02 (01) & {\cellcolor[HTML]{DEEBF2}} \color[HTML]{000000} 0.01 (06) & {\cellcolor[HTML]{3783BB}} \color[HTML]{F1F1F1} -0.15 (03) \\
ODR-SPGC & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 9.88 (02) & {\cellcolor[HTML]{9DCBE1}} \color[HTML]{000000} 46.22 (06) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 10.99 (02) & {\cellcolor[HTML]{9BC9E0}} \color[HTML]{000000} 36.01 (06) & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 73.90 (06) & {\cellcolor[HTML]{1C5C9F}} \color[HTML]{F1F1F1} 114.55 (06) & {\cellcolor[HTML]{276EB0}} \color[HTML]{F1F1F1} 75.21 (06) & {\cellcolor[HTML]{195696}} \color[HTML]{F1F1F1} 102.65 (05) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.00 (01) & {\cellcolor[HTML]{3A87BD}} \color[HTML]{F1F1F1} 0.03 (03) & {\cellcolor[HTML]{EFF3F5}} \color[HTML]{000000} 0.02 (09) & {\cellcolor[HTML]{2B73B3}} \color[HTML]{F1F1F1} -0.13 (02) \\
MODR-OLS-CD-IND & {\cellcolor[HTML]{4393C3}} \color[HTML]{F1F1F1} 10.45 (05) & {\cellcolor[HTML]{FBCEB7}} \color[HTML]{000000} 60.09 (09) & {\cellcolor[HTML]{4695C4}} \color[HTML]{F1F1F1} 11.71 (06) & {\cellcolor[HTML]{F8BFA4}} \color[HTML]{000000} 48.49 (09) & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} 90.59 (11) & {\cellcolor[HTML]{F9EDE5}} \color[HTML]{000000} 139.14 (10) & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} 91.93 (11) & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} 126.91 (11) & {\cellcolor[HTML]{FAE7DC}} \color[HTML]{000000} -0.03 (09) & {\cellcolor[HTML]{FBE6DA}} \color[HTML]{000000} -0.15 (06) & {\cellcolor[HTML]{DEEBF2}} \color[HTML]{000000} -0.01 (07) & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} -0.32 (10) \\
MODR-OLS-CD-$\Sigma$ & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} 11.88 (14) & {\cellcolor[HTML]{E17860}} \color[HTML]{F1F1F1} 67.75 (12) & {\cellcolor[HTML]{FCD3BC}} \color[HTML]{000000} 13.29 (14) & {\cellcolor[HTML]{DA6853}} \color[HTML]{F1F1F1} 54.82 (13) & {\cellcolor[HTML]{134C87}} \color[HTML]{F1F1F1} 72.06 (04) & {\cellcolor[HTML]{0C3D73}} \color[HTML]{F1F1F1} 111.61 (03) & {\cellcolor[HTML]{134C87}} \color[HTML]{F1F1F1} 73.26 (04) & {\cellcolor[HTML]{0D3F76}} \color[HTML]{F1F1F1} 100.64 (03) & {\cellcolor[HTML]{E0ECF3}} \color[HTML]{000000} -0.02 (05) & {\cellcolor[HTML]{3681BA}} \color[HTML]{F1F1F1} -0.03 (02) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 0.00 (01) & {\cellcolor[HTML]{71B0D3}} \color[HTML]{F1F1F1} -0.18 (04) \\
MODR-LASSO-CD-$\Sigma$ & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 13.09 (16) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 83.18 (16) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 14.95 (16) & {\cellcolor[HTML]{67001F}} \color[HTML]{F1F1F1} 66.14 (16) & {\cellcolor[HTML]{0E4179}} \color[HTML]{F1F1F1} 71.49 (03) & {\cellcolor[HTML]{0E4179}} \color[HTML]{F1F1F1} 112.08 (04) & {\cellcolor[HTML]{0E4179}} \color[HTML]{F1F1F1} 72.72 (03) & {\cellcolor[HTML]{0E4179}} \color[HTML]{F1F1F1} 100.79 (04) & {\cellcolor[HTML]{0A3B70}} \color[HTML]{F1F1F1} -0.00 (02) & {\cellcolor[HTML]{F7B799}} \color[HTML]{000000} -0.22 (10) & {\cellcolor[HTML]{BBDAEA}} \color[HTML]{000000} 0.01 (03) & {\cellcolor[HTML]{F8F2EF}} \color[HTML]{000000} -0.31 (08) \\
MODR-OLS-MCD-IND & {\cellcolor[HTML]{3F8EC0}} \color[HTML]{F1F1F1} 10.42 (04) & {\cellcolor[HTML]{FBCEB7}} \color[HTML]{000000} 60.06 (08) & {\cellcolor[HTML]{4393C3}} \color[HTML]{F1F1F1} 11.69 (05) & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 48.35 (08) & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} 90.37 (10) & {\cellcolor[HTML]{FAEAE1}} \color[HTML]{000000} 139.46 (11) & {\cellcolor[HTML]{F8BFA4}} \color[HTML]{000000} 91.79 (10) & {\cellcolor[HTML]{F9F0EB}} \color[HTML]{000000} 126.46 (09) & {\cellcolor[HTML]{E4EEF4}} \color[HTML]{000000} -0.02 (06) & {\cellcolor[HTML]{FBCEB7}} \color[HTML]{000000} -0.18 (09) & {\cellcolor[HTML]{A2CDE3}} \color[HTML]{000000} -0.01 (02) & {\cellcolor[HTML]{E6EFF4}} \color[HTML]{000000} -0.27 (06) \\
MODR-OLS-MCD-$\Sigma$ & {\cellcolor[HTML]{FAE7DC}} \color[HTML]{000000} 11.55 (12) & {\cellcolor[HTML]{DD7059}} \color[HTML]{F1F1F1} 68.64 (14) & {\cellcolor[HTML]{F9F0EB}} \color[HTML]{000000} 12.93 (13) & {\cellcolor[HTML]{D25849}} \color[HTML]{F1F1F1} 56.01 (14) & {\cellcolor[HTML]{08366A}} \color[HTML]{F1F1F1} 70.96 (02) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 110.49 (01) & {\cellcolor[HTML]{08366A}} \color[HTML]{F1F1F1} 72.19 (02) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 99.20 (01) & {\cellcolor[HTML]{F7F5F4}} \color[HTML]{000000} -0.02 (07) & {\cellcolor[HTML]{FBE3D4}} \color[HTML]{000000} -0.15 (07) & {\cellcolor[HTML]{D2E6F0}} \color[HTML]{000000} -0.01 (04) & {\cellcolor[HTML]{DEEBF2}} \color[HTML]{000000} -0.26 (05) \\
MODR-LASSO-MCD-$\Sigma$ & {\cellcolor[HTML]{DD7059}} \color[HTML]{F1F1F1} 12.28 (15) & {\cellcolor[HTML]{960F27}} \color[HTML]{F1F1F1} 78.74 (15) & {\cellcolor[HTML]{D86551}} \color[HTML]{F1F1F1} 14.01 (15) & {\cellcolor[HTML]{900D26}} \color[HTML]{F1F1F1} 62.86 (15) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 70.62 (01) & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 110.91 (02) & {\cellcolor[HTML]{053061}} \color[HTML]{F1F1F1} 71.85 (01) & {\cellcolor[HTML]{073467}} \color[HTML]{F1F1F1} 99.66 (02) & {\cellcolor[HTML]{F9C2A7}} \color[HTML]{000000} -0.06 (12) & {\cellcolor[HTML]{F5AA89}} \color[HTML]{000000} -0.24 (11) & {\cellcolor[HTML]{FDDCC9}} \color[HTML]{000000} -0.04 (12) & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} -0.41 (11) \\
MODR-OLS-LRA-IND & {\cellcolor[HTML]{96C7DF}} \color[HTML]{000000} 10.76 (07) & {\cellcolor[HTML]{FBCCB4}} \color[HTML]{000000} 60.27 (10) & {\cellcolor[HTML]{87BEDA}} \color[HTML]{000000} 12.01 (09) & {\cellcolor[HTML]{F8BB9E}} \color[HTML]{000000} 48.83 (10) & {\cellcolor[HTML]{F8BDA1}} \color[HTML]{000000} 90.66 (12) & {\cellcolor[HTML]{FAE8DE}} \color[HTML]{000000} 140.08 (13) & {\cellcolor[HTML]{F8BB9E}} \color[HTML]{000000} 92.10 (12) & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} 126.99 (12) & {\cellcolor[HTML]{FBE6DA}} \color[HTML]{000000} -0.03 (10) & {\cellcolor[HTML]{FDD9C4}} \color[HTML]{000000} -0.17 (08) & {\cellcolor[HTML]{EDF2F5}} \color[HTML]{000000} -0.02 (08) & {\cellcolor[HTML]{F9F0EB}} \color[HTML]{000000} -0.31 (09) \\
MODR-OLS-LRA-$\Sigma$ & {\cellcolor[HTML]{BBDAEA}} \color[HTML]{000000} 10.95 (08) & {\cellcolor[HTML]{E6866A}} \color[HTML]{F1F1F1} 66.68 (11) & {\cellcolor[HTML]{B3D6E8}} \color[HTML]{000000} 12.26 (10) & {\cellcolor[HTML]{DB6B55}} \color[HTML]{F1F1F1} 54.65 (12) & {\cellcolor[HTML]{FCD7C2}} \color[HTML]{000000} 89.00 (08) & {\cellcolor[HTML]{F9EEE7}} \color[HTML]{000000} 138.95 (09) & {\cellcolor[HTML]{FCD5BF}} \color[HTML]{000000} 90.48 (08) & {\cellcolor[HTML]{F7F5F4}} \color[HTML]{000000} 125.42 (08) & {\cellcolor[HTML]{FDDCC9}} \color[HTML]{000000} -0.04 (11) & {\cellcolor[HTML]{DE735C}} \color[HTML]{F1F1F1} -0.33 (12) & {\cellcolor[HTML]{F9EBE3}} \color[HTML]{000000} -0.03 (11) & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} -0.45 (12) \\
MODR-LASSO-LRA-$\Sigma$ & {\cellcolor[HTML]{C7E0ED}} \color[HTML]{000000} 11.00 (11) & {\cellcolor[HTML]{DF765E}} \color[HTML]{F1F1F1} 68.00 (13) & {\cellcolor[HTML]{D8E9F1}} \color[HTML]{000000} 12.51 (11) & {\cellcolor[HTML]{DD7059}} \color[HTML]{F1F1F1} 54.21 (11) & {\cellcolor[HTML]{FCD3BC}} \color[HTML]{000000} 89.29 (09) & {\cellcolor[HTML]{FAEAE1}} \color[HTML]{000000} 139.65 (12) & {\cellcolor[HTML]{FBD0B9}} \color[HTML]{000000} 90.72 (09) & {\cellcolor[HTML]{F9EFE9}} \color[HTML]{000000} 126.57 (10) & {\cellcolor[HTML]{F5AA89}} \color[HTML]{000000} -0.08 (13) & {\cellcolor[HTML]{A51429}} \color[HTML]{F1F1F1} -0.57 (15) & {\cellcolor[HTML]{F5A886}} \color[HTML]{000000} -0.08 (13) & {\cellcolor[HTML]{D6604D}} \color[HTML]{F1F1F1} -0.56 (15) \\
\bottomrule
\end{tabular}
}
\caption{{\color{black} Price spike and spread analysis for the CRPS and the miscoverage $\text{MC}_{0.95}$. A spike event is defined if the min/max price of a day exceeds the 5\%{} resp. 95\%-quantile of all min/max prices. A large spread event is defined as the min-to-max spread exceeding the 90\%-quantile of all min-to-max spreads. Numbers in brackets give the model ranking.}}\label{tab:results_spike}
\end{center}
\end{table}
{\color{black}
Before concluding the results, we want to emphasize some limitations of this study. \begin{itemize}
\item We acknowledge that the proposed models yield a trade-off between accuracy in the marginal distributions and correct modeling of the dependence structure. This trade-off is reflected also in the disconnect between marginally dominated scoring rules (CRPS and ES) and the LS, DSS, and VS, and therefore, we yield somewhat ambiguous results.
\item While the results of the spike analysis are encouraging, the parametric assumption on the distributional family might be too restrictive in some cases and could be complemented with non-parametric approaches, or with multivariate extreme value models to capture the tail behavior. This could also help to improve the forecasting performance in the marginal distributions and hence CRPS and ES.
\end{itemize}
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}
\subsection{Computation Time}
}\label{sec:results_comptime}
\begin{table}[htb]
\begin{minipage}[t]{.7\linewidth}
\vspace{0pt}
\resizebox{\linewidth}{!}{
\begin{tabular}{lrrrrr}
\toprule
Model & Initial Fit & Avg. Update & Total Time & Est. Speedup \\
\midrule
LEAR-N(0, $\sigma$) & {\cellcolor[HTML]{F7FBFF}} \color[HTML]{000000} 1.23 & {\cellcolor[HTML]{E7F1FA}} \color[HTML]{000000} 0.06 & {\cellcolor[HTML]{F2F8FD}} \color[HTML]{000000} 103.80 & {\cellcolor[HTML]{B5D4E9}} \color[HTML]{000000} $\times$ 21 \\
LEAR-N(0, $\Sigma$) & {\cellcolor[HTML]{F7FBFF}} \color[HTML]{000000} 1.23 & {\cellcolor[HTML]{E7F1FA}} \color[HTML]{000000} 0.06 & {\cellcolor[HTML]{F2F8FD}} \color[HTML]{000000} 103.80 & {\cellcolor[HTML]{B5D4E9}} \color[HTML]{000000} $\times$ 21 \\
LEAR-CP & {\cellcolor[HTML]{F7FBFF}} \color[HTML]{000000} 1.23 & {\cellcolor[HTML]{E7F1FA}} \color[HTML]{000000} 0.06 & {\cellcolor[HTML]{F2F8FD}} \color[HTML]{000000} 103.80 & {\cellcolor[HTML]{B5D4E9}} \color[HTML]{000000} $\times$ 21 \\
LEAR-GARCH & {\cellcolor[HTML]{F3F8FE}} \color[HTML]{000000} 1.44 & {\cellcolor[HTML]{89BEDC}} \color[HTML]{000000} 0.27 & {\cellcolor[HTML]{A1CBE2}} \color[HTML]{000000} 501.31 & {\cellcolor[HTML]{F7FBFF}} \color[HTML]{000000} $\times$ 5 \\
ODR-IND & {\cellcolor[HTML]{6AAED6}} \color[HTML]{F1F1F1} 55.31 & {\cellcolor[HTML]{82BBDB}} \color[HTML]{000000} 0.30 & {\cellcolor[HTML]{92C4DE}} \color[HTML]{000000} 598.78 & {\cellcolor[HTML]{2272B6}} \color[HTML]{F1F1F1} $\times$ 169 \\
ODR-GC & {\cellcolor[HTML]{6AAED6}} \color[HTML]{F1F1F1} 55.32 & {\cellcolor[HTML]{82BBDB}} \color[HTML]{000000} 0.30 & {\cellcolor[HTML]{92C4DE}} \color[HTML]{000000} 599.61 & {\cellcolor[HTML]{2272B6}} \color[HTML]{F1F1F1} $\times$ 168 \\
ODR-SPGC & {\cellcolor[HTML]{6AAED6}} \color[HTML]{F1F1F1} 55.46 & {\cellcolor[HTML]{60A7D2}} \color[HTML]{F1F1F1} 0.45 & {\cellcolor[HTML]{71B1D7}} \color[HTML]{F1F1F1} 870.45 & {\cellcolor[HTML]{3787C0}} \color[HTML]{F1F1F1} $\times$ 116 \\
MODR-OLS-CD-IND & {\cellcolor[HTML]{7DB8DA}} \color[HTML]{000000} 39.08 & {\cellcolor[HTML]{D6E6F4}} \color[HTML]{000000} 0.08 & {\cellcolor[HTML]{D8E7F5}} \color[HTML]{000000} 190.40 & {\cellcolor[HTML]{08478D}} \color[HTML]{F1F1F1} $\times$ 375 \\
MODR-OLS-CD-$\Sigma$ & {\cellcolor[HTML]{4594C7}} \color[HTML]{F1F1F1} 131.51 & {\cellcolor[HTML]{A4CCE3}} \color[HTML]{000000} 0.20 & {\cellcolor[HTML]{A3CCE3}} \color[HTML]{000000} 494.41 & {\cellcolor[HTML]{083877}} \color[HTML]{F1F1F1} $\times$ 487 \\
MODR-LASSO-CD-$\Sigma$ & {\cellcolor[HTML]{0E59A2}} \color[HTML]{F1F1F1} 752.57 & {\cellcolor[HTML]{084A91}} \color[HTML]{F1F1F1} 2.41 & {\cellcolor[HTML]{09529D}} \color[HTML]{F1F1F1} 5156.91 & {\cellcolor[HTML]{0E59A2}} \color[HTML]{F1F1F1} $\times$ 267 \\
MODR-OLS-MCD-IND & {\cellcolor[HTML]{77B5D9}} \color[HTML]{000000} 43.97 & {\cellcolor[HTML]{E0ECF8}} \color[HTML]{000000} 0.07 & {\cellcolor[HTML]{DFEBF7}} \color[HTML]{000000} 163.73 & {\cellcolor[HTML]{083776}} \color[HTML]{F1F1F1} $\times$ 491 \\
MODR-OLS-MCD-$\Sigma$ & {\cellcolor[HTML]{4493C7}} \color[HTML]{F1F1F1} 134.60 & {\cellcolor[HTML]{B2D2E8}} \color[HTML]{000000} 0.16 & {\cellcolor[HTML]{ABD0E6}} \color[HTML]{000000} 434.25 & {\cellcolor[HTML]{08306B}} \color[HTML]{F1F1F1} $\times$ 567 \\
MODR-LASSO-MCD-$\Sigma$ & {\cellcolor[HTML]{084184}} \color[HTML]{F1F1F1} 1530.21 & {\cellcolor[HTML]{083370}} \color[HTML]{F1F1F1} 3.57 & {\cellcolor[HTML]{083877}} \color[HTML]{F1F1F1} 8066.77 & {\cellcolor[HTML]{084B93}} \color[HTML]{F1F1F1} $\times$ 347 \\
MODR-OLS-LRA-IND & {\cellcolor[HTML]{9DCAE1}} \color[HTML]{000000} 21.62 & {\cellcolor[HTML]{F7FBFF}} \color[HTML]{000000} 0.04 & {\cellcolor[HTML]{F7FBFF}} \color[HTML]{000000} 92.11 & {\cellcolor[HTML]{084082}} \color[HTML]{F1F1F1} $\times$ 429 \\
MODR-OLS-LRA-$\Sigma$ & {\cellcolor[HTML]{2777B8}} \color[HTML]{F1F1F1} 301.48 & {\cellcolor[HTML]{3A8AC2}} \color[HTML]{F1F1F1} 0.78 & {\cellcolor[HTML]{3F8FC5}} \color[HTML]{F1F1F1} 1735.92 & {\cellcolor[HTML]{08509B}} \color[HTML]{F1F1F1} $\times$ 317 \\
MODR-LASSO-LRA-$\Sigma$ & {\cellcolor[HTML]{08306B}} \color[HTML]{F1F1F1} 2469.66 & {\cellcolor[HTML]{08306B}} \color[HTML]{F1F1F1} 3.81 & {\cellcolor[HTML]{08306B}} \color[HTML]{F1F1F1} 9440.61 & {\cellcolor[HTML]{083979}} \color[HTML]{F1F1F1} $\times$ 478 \\
\bottomrule
\end{tabular}
}
\end{minipage}
\begin{minipage}[t]{.29\linewidth}
\vspace{0pt}
\includegraphics[width=\linewidth]{figures/computation_times.png}
\end{minipage}
\caption{Computation times. All timings are in seconds. The out-of-sample data for the forecasting study consists of 1831 days. We update the 24-dimensional distributional regression model on each day. All experiments are run on a standard laptop (Intel Core i7 (16 Threads, 4.9 GHz), 32GB RAM). Estimated speed-ups are calculated by taking $\text{Speedup} = (\text{Initial Fit} \times T) / \text{Total Time}$. {\color{black} The figure on the right shows the share of the time spent on the initial fit and on the out-of-sample updating}.}
\label{tab:computation_times}
\end{table}
Table \ref{tab:computation_times} 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 \cite{hirsch2024online} for the univariate online distributional regression case in {\color{black} direct} comparison {\color{black} between online estimation and repeated batch estimation}.
\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{figures/eff_frontier.pdf}
\caption{Efficient frontier of computation time against average log-score for the online multivariate distributional regression model for different update frequencies in the batch and online case. Note that the $x$-axis is in log-scale.}
\label{fig:results_comptime_eff}
\end{figure}
{\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{fig:results_comptime_eff}.}} We compare the results for online learning, rolling and expanding window batch estimation in Figure~\ref{fig:results_comptime_eff}. {\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{tab:computation_times}.
}
{\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
\section{Discussion and Conclusion}\label{sec:conclusion}
Distributional learning algorithms such as GAMLSS and deep distributional networks have been used successfully for probabilistic electricity price forecasting \cite[PEPF, see e.g.][]{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 \citep[JPB, see][]{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, \texttt{sklearn}-like API to facilitate the usage by other researchers and data scientists \citep{pedregosa2011scikit} and contributed the implementation to the \texttt{ondil} package \citep{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 \citep{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{sec:online_estimation}. 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
\section*{Acknowledgments}
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.
\section*{Declaration of Interest}
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.
\section*{Data Statement}
The data used in this paper has been provided by \cite{lipiecki2024postprocessing} and cannot be shared further by the author. Therefore, for any requests with respect to the data, please contact \cite{lipiecki2024postprocessing} directly. A shorter data set of the same explanatory variables is available at the \texttt{GitHub} repository of \cite{marcjasz2023distributional}.\footnote{See: \url{https://github.com/gmarcjasz/distributionalnn}.} For the differences between both data sets please consult \cite[][Section 2 and Section 5.1.3]{lipiecki2024postprocessing}.
\section*{Generative AI Statement}
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
\bibliographystyle{elsarticle-harv}
\bibliography{references}
\newpage