Multivariate Probabilistic CRPS Learning with an Application to Day-Ahead Electricity Prices

Multivariate Probabilistic CRPS Learning with an Application to Day-Ahead Electricity Prices

Introduction

Forecast combination (sometimes referred to as expert aggregation or ensembling) has recently gained much traction. We know from theory that combination methods work well to combine different but well-performing model classes cesa2006prediction. As gaillard2016additive pointed out, it is always recommended to use different classes of models, e.g., regression and time series type models, neural network models, decision tree learning models, and other machine learning and artificial intelligence methods.

{This paper proposes a novel online updating scheme for combining the marginals of the corresponding multivariate distribution across quantiles (also referred to as horizontal aggregation). We know from Sklar's theorem that we can decompose any multivariate distribution into the marginals and a copula. That is, we can improve the marginals (i.e., by using a strictly proper scoring rule like the CRPS) while leaving the copula untouched. In consequence, we require only the reporting of the forecasted marginal distribution. The proposed method considers dependencies between the combination weights across quantiles and marginals through a simple but flexible smoothing procedure. We assume a basic metric or spatial structure in the multivariate dimension. Such a metric structure is present when forecasting a univariate time series several steps ahead or predicting one-dimensional spatial data.}

Online learning algorithms are particularly attractive for forecasting where frequent short-term forecasts are essential for the application domain (e.g., energy, weather, finance, retail). The proposed algorithm generalizes the probabilistic CRPS learning framework presented in berrisch2021crps. It is based on exponential weighted averaging (EWA) and yields optimal asymptotic convergence rates with respect to the best individual forecast and the best convex combination of all forecasts wintenberger2017optimal.

Considerable research on forecasting combination already exists. bordignon2013combining, nowotarski2014empirical, avci2018managing combine point-forecasts using various batch methods. marcjasz2020probabilistic, Serafin2019averaging apply batch methods to probabilistic forecast combination. Some authors also applied online learning algorithms for point forecasting nowotarski2016improving and probabilistic forecasting gaillard2016additive, gonzalez2021new. The work above focuses on developing distinct forecasting models and on combination methods. gaillard2015forecasting discuss how model development can be optimized in the framework of aggregation of experts.

In electricity price forecasting, dynamic aggregation techniques, where the combination weights are adjusted based on past performance, tend to perform better than simple constant weight techniques gaillard2015forecasting, marcjasz2018selection, maciejowska2020pca. However, they consider multivariate updating schemes that use the same weight for all time series. Most other work in energy forecasting considers all time series to be independent and therefore combines forecasts separately bordignon2013combining, nowotarski2016improving, nitka2023combining. {Neither approach considers possible dependencies of combination weights between marginals. Consequently, we can expect potential improvements by exploiting this metric structure of electricity prices by considering updating schemes that assign different weights to all neighboring price forecasts of the day and considering possible dependencies between combination weights. Of course, the same logic applies to other areas of application.}

The contributions of this manuscript are manifold:

The remainder of this paper is structured as follows. Section (ref) discusses the general multivariate probabilistic combination setting and discusses CRPS learning using quantile regression. Section (ref) presents the proposed multivariate generalization of online CRPS learning and summarizes its asymptotic properties. Additionally, we discuss possible extensions of the proposed method. Those extensions to the core algorithm add hyperparameters that have to be specified. Therefore, we elaborate on two possible strategies for hyperparameter tuning in Section (ref). Section (ref) continues with an empirical application of the proposed algorithm. We apply the methodology to multivariate probabilistic forecasts of Day-Ahead power prices. We discuss the data, elaborate on the specific algorithms we consider, and present a detailed analysis of the obtained results. Section (ref) discusses limitations, introduces potential enhancements, and concludes.

Multivariate CRPS Learning

The combination setting

In this paper, we consider the combination of multivariate probabilistic forecasts. {In particular, we consider a setting where the forecasts are given as quantiles of all marginals of a multivariate distribution.} berrisch2021crps show that pointwise forecast combinations potentially outperform standard methods where weights are constant over all distribution quantiles. We apply this idea to a multivariate setting by computing weights depending on the quantile and the marginals. First, we discuss batch learning methods and propose a dimension reduction technique that bridges the gap between flexible pointwise and robust constant procedures. Afterward, we show how the proposed online learning algorithm of berrisch2021crps can be extended for combining {the marginals of} multivariate probabilistic forecasts.

{Let $\widehat{\boldsymbol F}_{t} = (\widehat{F}_{t,1},\ldots, \widehat{F}_{t,K})$ be a vector of $K$ univariate distributions representing the marginal distribution of the corresponding multivariate distribution, resp. the set of experts that we want to combine.} We consider the combination across quantiles (also known as horizontal aggregation):

We evaluate the performance using the cumulative CRPS over all marginals. Therefore, the weights shall be chosen to minimize the cumulative CRPS of all marginals. We can approximate the CRPS by the sum over Quantile Losses ($\operatorname*{QL}$)

for an equidistant dense grid $\boldsymbol{\mathcal{P}} = ( p_1,\ldots, p_P )$ with $p_i<p_{i+1}$ and $p_{i+1}-p_i = h$ for all $p$. Clearly, $P\to \infty$ induces $h \to0$, $p_1\to0$, $p_P\to1$ and the approximation converges to the CRPS gneiting2011making, gneiting2011quantiles. marcjasz2022distributional omitted the scaling factor of 2 in equation (ref) as it does not affect the optimization, and there is no natural interpretation of the CRPS. We follow this approach to ensure comparability of the results.

This relationship enables us to compute pointwise weights based on quantile losses. We can extend this idea by optimizing weights not only depending on the quantile $p$ but also on the marginal $d$:

We are interested in setting $w_{t,k}$ such that the CRPS of $\widetilde{F}_{t,d}$ is minimized.

CRPS learning using quantile regression

Pointwise CRPS learning has the potential to outperform standard CRPS learning methods. However, the best pointwise weights in (ref) must be estimated. Theoretically, a pointwise approach has to be applied to all probabilities $p\in(0,1)$ and all marginals $\boldsymbol{\mathcal{D}} = (1, 2, \ldots, D)$ such that the bivariate weight function $\boldsymbol w_{t,k}$ can be specified. However, we can never evaluate infinitely many values for $p$. On the same page, the computation may be infeasible if $D$ is very large. Therefore, we must consider some finite-dimensional representation for the weight functions $\boldsymbol w_{t,k}$. A suitable option is representing the weight functions $\boldsymbol w_{t,k}$ using a finite-dimensional representation using splines. Bivariate splines are a suitable option in this scenario. We can express them as follows:

This is essentially the same as univariate splines with $L$-dimensional parameter vector $(\beta_1,\ldots,\beta_L)'$. However, the support of $f$ is 2-dimensional. Thus, we need many more basis functions $L$ to have a suitable description of $f$. A popular way to describe the bivariate basis function $\boldsymbol \varphi_l$ in (ref) is to assume a tensor structure mclean2014functional, wood2017gen. In the bivariate case, the spline function is a product of two univariate ones. In addition, $\varphi_{1,l}$ and $\varphi_{2,l}$ are usually chosen such that $\varphi_{1,l_1}$ interacts with each of the considered basis functions $\varphi_{2,l_2}$. This, allows to renumerate the problem such that $l=(l_1,l_2)$, and yields

This can be used to express the bivariate weight function as a product of the $\widetilde{P} \times \widetilde{D}$ parameter matrix $\boldsymbol \beta_{t,k}$ and the bivariate basis represented by $\boldsymbol{\varphi}^{\text{mv}}$ and $\boldsymbol{\varphi}^{\text{pr}}$:

Given $T$ historic forecasts $\widehat{F}^{-1}_{t,d,k}(p)$, the cooresponding realizations $Y_{t,d}$, and the index of marginals $\boldsymbol{\mathcal{D}} = (1,2, \ldots, D)$ we can estimate the $\widetilde{D} \times \widetilde{P} \times K$-dimensional parameter tensor $\boldsymbol \beta_{t}$ by minimizing the corresponding CRPS using (ref):

The second line uses the shift-invariance of the quantile loss and quantile regression notation $\rho_p(z) = \operatorname*{QL}_p(0,z) = z(p-\mathbb{1}\{z< 0\})$ koenker2017handbook.

Still, computing (ref) requires the evaluation of all distribution forecasts. As discussed, this is often not possible in practice. If we restrict the evaluation to a grid of probabilities $\boldsymbol{\mathcal{P}}$ problem (ref) simplifies with (ref) to

In general, quantile regression problems can be solved efficiently using linear programming solvers koenker2017handbook. However, (ref) is not a simple quantile regression problem, but a joint quantile regression sangnier2016joint, chun2016graphical. The parameters $\beta_{t,j,l,k}$ are active for multiple quantiles. Thus, adequate estimation requires solving the optimization problem for $\widetilde{D} \times \widetilde{P} \times K$ parameters, which can be computationally costly if $\widetilde{D}$, $\widetilde{P}$, and $K$ are large.

However, if we choose both basis $\boldsymbol \varphi^\text{pr}$ and $\boldsymbol \varphi^\text{mv}$ so that $\boldsymbol \varphi_i^\text{mv}=\mathbb{1}{\{d_i\}}$ on $\boldsymbol{\mathcal{D}}$ for $d_i\in \boldsymbol{\mathcal{D}}=(d_1,\ldots, d_D)$ and $\boldsymbol \varphi_i^\text{pr}=\mathbb{1}{\{p_i\}}$ on $\boldsymbol{\mathcal{P}}$ for $p_i\in \boldsymbol{\mathcal{P}}=(p_1,\ldots, p_P)$ then (ref) can be disentangled into $\widetilde{D} \times \widetilde{P}$ separate quantile regression problems. This is

for $p\in \boldsymbol{\mathcal{P}}$ and $d \in \boldsymbol{\mathcal{D}}$ where $\widehat{F}^{-1}_{t,h,k}(p)$ are the experts for the $p$-quantile and the $d$-marginal.

Quantile regression (ref) will lead to linear optimality on $\boldsymbol{\mathcal{D}}$ and $\boldsymbol{\mathcal{P}}$, as long as standard regularity conditions required for the quantile regression are satisfied koenker2001quantile. However, we might assume further restrictions to reduce the estimation risk, e.g., the solution is a convex combination. taylor1998combining discussed many related plausible restrictions for quantile combination concerning bias correction, positivity, and affinity, among others.

A potential issue of pointwise algorithms is quantile crossing. This problem occurs if we have $\widetilde{F}_{t,d}^{-1}(p_i) > \widetilde{F}_{t,h}^{-1}(p_j)$ for some $p_i,p_j\in(0,1)$ with $p_i<p_j$. {In this case, we recommend rearranging the predictions as sorting is known to improve the forecasting performance chernozhukov2010quantile.}

Multivariate Online CRPS Learning

Batch-learning approaches, like quantile regression, evaluate the entire history for estimating new combination weights, which is computationally costly. Therefore, we suggest to use online learning methods instead.

Online learning is often called prediction under expert advice. In this context, experts refer to the models producing the predictions (or predictive distributions). The person or model that combines the experts' predictions is called forecaster. A key element of online learning methods is (cumulative) regret. It is defined as:

i.e., the cumulative difference between the loss of the expert's predictions $\widehat{X}_{t,k}$ and the prediction of the forecaster $\widetilde{X}_{t}$ for a loss function $\ell$. $\widetilde{X}_{t,k}$ might be a predicted quantile $\widehat{F}^{-1}_{t,k}(p)$ of expert $k$ as discussed in the previous section. $R_{t,k}$ is called regret because it indicates how much the forecaster regrets not following the experts' advice cesa2006prediction. With (ref), we can formulate the EWA:

where $K$ refers to the number of experts, $w_{0,k}$ to the initial weights of an expert $k$, and $\eta$ to the learning rate, which defines how fast the weights adjust to changes in the regret cesa2006prediction. We can express this aggregation rule in terms of past weights and the loss suffered by the experts (right-hand side of (ref)). This highlights that there is no need for evaluating the entire history when adjusting weights.

EWA yields optimal convergence rates of ${\mathcal{O}}(T)$ towards the best expert for exp-concave loss functions cesa2006prediction. It means that the algorithm's performance (in terms of risk) is asymptotically not worse than the performance of the best expert. A more ambitious property that can also be satisfied is the convex aggregation property. It ensures that the risk of the algorithm is not worse than the risk of the best convex combination of the experts. For an algorithm to satisfy this property, the gradient trick is needed devaine2013forecasting. For exp-concave losses, this gives optimal convergence rate ${\mathcal{O}}(\sqrt{T})$ with respect to the best convex combination of the experts cesa2006prediction. This property also holds for losses that satisfy some Bernstein condition, such as the MAE, when algorithms like Bernstein Online Aggregation (BOA) are used. These algorithms use regularized updating techniques to improve convergence and stability properties wintenberger2017optimal.

berrisch2021crps adapted BOA to probabilistic problems. The new algorithm is called CRPS learning because it optimizes the CRPS of the target distribution using pointwise optimization on a grid of quantiles. The weights can vary over time and in different parts of the distribution. CRPS learning still maintains the fast convergence of BOA. This algorithm for combining $\widehat{\boldsymbol X}_{t}=(\widehat{X}_{t,1},\ldots, \widehat{X}_{t,K})$ to $\widetilde{X}_{t}=\boldsymbol w_{t-1}'\widehat{\boldsymbol X}_{t}$ can be summarized as follows:

where $\boldsymbol x^+$ and $\boldsymbol x^-$ denote the elementwise positive and negative parts of $\boldsymbol x$ and $\odot$ the elementwise product (Hadamard product). The learning rate $\boldsymbol \eta_{t}$ determines the weight adjustment speed. It depends on the bound estimator $\boldsymbol E_{t}$ and $\boldsymbol V_t$, which is an estimator for the variance. The algorithm describes how weights are calculated on a full quantile grid $\boldsymbol{\mathcal{P}}$. First, the instantaneous regret is calculated (ref). Then the learning rate ((ref)-(ref)) is adjusted. In (ref) the cumulative regret is calculated. Afterward, we calculate the weights (ref). Finally, the forecaster uses $\boldsymbol w_{t}$ to calculate $\widetilde{X}_{t+1}$ and starts over with (ref).

Several extensions of online learning algorithms were proposed in the literature. However, they can also be applied in standard Batch-Learning settings. The benefits of these extensions have been confirmed in empirical studies. Some extensions, like shrinkage operators, are also valuable to nest specific weighting strategies into the learning algorithm.

Smoothing

As mentioned, we apply the general CRPS learning idea to multivariate data. Therefore, we adopt the two weight-smoothing methods of the original CRPS Learning algorithm. The first consists of the dimension reduction method using basis matrices. The approach is analogous to the idea discussed in Section (ref). Using a bivariate basis, we can reduce the dimensionality of the instantaneous regret from $D \times P$ to $\widetilde{D} \times \widetilde{P}$:

As usual, we can use this reduced regret to carry out the online learning algorithm. After obtaining weights (we refer to them as $\boldsymbol \beta_{t,k}$) on this condensed version of the regret, we can utilize the basis matrices once again to obtain weights in our original dimensions of interest $\boldsymbol w_{t,k} = {\boldsymbol B^\text{mv}} \boldsymbol \beta_{t,k} {\boldsymbol B^\text{pr}}'$.

{This relatively simple method yields a powerful property: It bridges the gap between purely pointwise weight optimization based on quantiles and the constant approach where a single weight is optimized with respect to the CRPS. That means we can move from a setting with low flexibility (i.e., a few parameters) and low estimation risk to a very flexible one (with many parameters) at the price of high estimation risk.}

Another option is to smooth the weights using penalized smoothing. This method can be applied after the estimation, i.e., after the updating step. Hereby we consider two sets of bounded basis functions $\boldsymbol \psi^{\text{\text{mv}}}=(\psi_1,\ldots, \psi_{D})$ and $\boldsymbol \psi^{\text{\text{pr}}}=(\psi_1,\ldots, \psi_{P})$ on $(0,1)$ that we will use for penalized smoothing.

Then the weights can be represented by

with parameter matix $\boldsymbol b_{t,k}$. We estimate $\boldsymbol b_{t,k}$ by penalized $L_1$- and $L_2$-smoothing which minimizes

for each $k$ given $\boldsymbol \beta_{t,k}$ with differential operator $\mathcal{D}_q$ of order $q$. The differential order characterizes the smoothing penalty, and $\lambda\geq 0$ characterizes the roughness penalty. Typically, $q=2$ is considered along with cubic B-Splines to penalize for roughness wang2011smoothing, wood2017generalized. However, we prefer using $q<2$ here. This smoothes towards constant weights over $\boldsymbol{\mathcal{P}}$ for $\lambda\to \infty$ and not towards a linear relationship between weights and probabilities as for $q=2$. As pointed out in berrisch2021crps no argument supports shrinkage towards a linear relationship. In contrast, shrinkage towards constant weights yields the non-pointwise CRPS-learning theory of constant weight functions. However, let us remark that the penalized smoothing approach with $\lambda\to \infty$ yields a different result than the simple basis smoothing approach mentioned before with $\boldsymbol \varphi = \varphi_1 \equiv 1$.

In applications, we only apply this function bases approach on finite grids of probabilities $\boldsymbol{\mathcal{P}}=(p_1, \ldots, p_P)$ and a finite number of marginals $\boldsymbol{\mathcal{D}}=(1, \ldots, D)$. If we consider B-Spline basis functions $\boldsymbol \psi^\text{mv}$ and $\boldsymbol \psi^\text{pr}$, then an explicit solution based on ordinary least squares exists for (ref). This explicit solution has a ridge regression representation. The smoothed weights matrix $\boldsymbol w_{t,k}$ is then given by

with basis matrices $\boldsymbol B^\text{mv} = \boldsymbol \psi^\text{mv}\left( \boldsymbol{\mathcal{D}} \right)$ and $\boldsymbol B^\text{pr} = \boldsymbol \psi^\text{pr}\left( \boldsymbol{\mathcal{P}} \right)$, penalty matrices $\boldsymbol D^\text{mv}_q$ and $\boldsymbol D^\text{pr}_q$. We can easily compute the penalty matrix if the b-spline basis has equidistant knots. Hereinafter, we distinguish $\boldsymbol D^S_q$ and $\boldsymbol D^G_q$, which refer to the equidistant case and the general case where knot placement does not have to be equidistant, respectively. Let $\Delta$ denote the matrix difference operator:

Now, $\boldsymbol D^{S}_q$ can be easily computed as $\boldsymbol D^{S}_q= \Delta^q \boldsymbol I$. The computation of $\boldsymbol D^{G}_q$ is more intricate since non-equidistant knots are permitted. The calculation involves an additional weighting step with respect to the non-equidistant distribution of the knots. We elaborate on this topic briefly since the literature is surprisingly scarce li2022general. The required difference matrix can be computed as

where $W_q$ are weighting matrices that depend on the order of the B-Spline basis, denoted as $o$, and the knots. Let $J$ denote the number of inner knots. The dimension of the difference matrices $\Delta$ in (ref) depend on $W_q$ are $(J+o-q)\times(J+o-q+1)$. The total number of knots will be $J + 2o$. We can specify the weighting matrices as:

The quantity $o-q$ represents the lag used to differentiate the knots. If the knots are equidistant, then $\bm{W}_q$ will be proportional to the identity matrix $\boldsymbol I$. Therefore, it nets the standard P-Spline, which uses $\boldsymbol D^{S}_q= \Delta^q \boldsymbol I$. This gives rise to the general P-Spline estimator. However, $\bm{W}_q$ is only proportional to $\boldsymbol I$ rather than equal to it. Therefore, scaling needs to be applied for the results of both estimators to coincide. The scaling can be applied to lambda, the penalty matrix, or the difference matrix. To state this formally, let $\bm{P}_q^S = {\bm{D}_q^S}'\bm{D}_q^S$ and $\bm{P}_q^G = {\bm{D}_q^G}'\bm{D}_q^G$ denote the penalty terms of the standard and general P-Spline estimators. The scaling factor with respect to the penalty $\bm{P}_q^G$ is $\left(\operatorname{Tr}\left(\bm{W}_q\right)/(J+o-q)\right)^{2q}$ so the following relation holds:

This is only valid for equidistant B-Splines. For non-equidistant B-Splines, the generalized P-Spline is the only appropriate estimator. However, $\bm{P}_q^G$ should always be scaled to ensure the comparability between lambda values in equidistant and non-equidistant situations.

For notational brevity, we denote the first and last part of (ref) as $\boldsymbol{\mathcal{H}}^\text{mv}$ and $\boldsymbol{\mathcal{H}}^\text{pr}$, respectively, the so-called hat matrices. Fortunately, they do not depend on time-varying components; therefore, we can compute them prior to the main online learning task, which yields a great reduction in the algorithm's computational complexity.

Knot placement for Smoothing Splines

For both smoothing approaches discussed above, the knots of the B-Spline basis must be placed. A well-established approach is placing plenty of equidistant knots. However, as discussed above, non-equidistant knot placement is valid if the penalty is defined accordingly. We consider equidistant and non-equidistant knots. Thereby, the non-central beta distribution with the following parameterization is used for distributing the knots:

Where $I_x$ is the incomplete beta function, $a$ and $b$ are shape parameters, and $c$ is the non-centrality parameter johnson1995continuous. Algorithm (ref) describes the knot placement in detail. It returns equidistant knots if $\mu = 0.5$, $\sigma = 1$, $c = 0$ and the tailweight parameter $\tau = 1$.

Figure (ref) shows B-Spline Basis' for different knot placements for the inputs of Algorithm (ref).

Shrinkage operators and Forgetting

Shrinkage operators are well-known in statistical learning theory. They help to reduce the overfitting problem by shrinking a solution. The P-Spline smoothing discussed above can also be interpreted as a shrinkage operator. However, simple shrinkage operators can also be applied to $\boldsymbol \beta_{t}$. We consider three additional shrinkage operators: the fixed share operator $\mathcal{F}$, the soft-thresholding operator $\mathcal{S}$, and the hard-thresholding operator $\mathcal{H}$. They are defined as

for some $\phi\in[0,1]$, $\nu\geq0$ and $\kappa\geq0$. The fixed share operator shrinks towards the naive combination. This is preferable if no prior information on the experts' performance is available. For some shrinkage problems, there are theoretical guarantees for improvements tu2011markowitz, cesa2012mirror. Applications in the context of online learning include, e.g., cesa2012mirror, gonzales2021new. The thresholding operators $\mathcal{S}$ and $\mathcal{H}$ were also considered in online learning contexts previously dalalyan2012sharp, gaillard2017sparse. Applying thresholds leads to sparse solutions. Both appear in several situations for specific linear model estimators. Most notably, soft-thresholding is the key operator in the coordinate descent algorithm for estimating the lasso friedman2007pathwise. Applying any threshold operator potentially violates affinity constraints (incl. the convexity constraint). Therefore, projections to the desired solution space should be applied.

As mentioned, cumulative regret is a key element in online learning. However, in settings with a long history, there might be structural breaks in the data. These breaks motivate the introduction of the forgetting factor. It means that only a limited amount of the old cumulative regret is considered for adjusting the weights. In other words, the algorithm forgets about some part of the past performance. In online learning, usually, exponential forgetting is chosen guo2018online, messner2019online,ziel2021smoothed. The Regret with a forgetting factor $\theta\in[0,1]$ is formally defined as

where $\theta = 0$ correspons to no forgetting. Optimal values for the forgetting factor $\theta$ are usually close to $0$. The forget should be applied to all hidden state variables in sophisticated online learning procedures like BOA.

Full Model and Hyperparameter Optimization

Algorithm (ref) shows the multivariate online CRPS-Learning algorithm. This includes all extensions discussed in Subsections (ref) and (ref). Considering all extensions, this algorithm contains five general hyperparameters (the forget rate $\phi$, the parameters of the shrinkage operators $\theta$, $\kappa$, $\nu$) as well as 30 hyperparameters concerning the design of basis and hat matrices.

The algorithm is versatile as it handles several special cases discussed in the literature. One such case is the uniform combination, also known as the naive combination. This can be calculated using the Fixed-Share operator $\mathcal{F}_\phi$ with $\phi=1$, resulting in uniform weights. There are more efficient ways of calculating uniform weights. However, adjusting the value of $\phi$ allows a smooth transition from the uniform solution to the solution computed by BOA. Another typical particular case is constant weights, where each expert receives a specific weight. This can be calculated by setting both basis matrices, $B^\text{mv}$, and $B^\text{pr}$, to the unity Vector of length $D$ and $P$, respectively. This leads to weights without variation across marginals and probabilities (Constant). Setting only one of the basis matrices to the unity vector will result in weights that are constant over either marginals (Constant Mv) or probabilities (Constant Pr). Additionally, setting all smoothing matrices to identity produces pointwise weights, and optimizing $\lambda$ in the hat matrices concerning the predictive CRPS results in possibly smoothed weights. These cases are illustrated in Figures (ref)-(ref).

The extensions discussed in Subsections (ref) and (ref) require the specification of various hyperparameters. There are many possible hyperparameters to choose from, and we do not have any prior information on the best values. This means that it is impossible to test all combinations of these parameters in each iteration of the forecasting task. The latter would be ideal, but it is impractical due to the required computational resources. As a result, we need to use other, less demanding methods for tuning these hyperparameters. In this paper, we will utilize three approaches for tuning.

The first approach to hyperparameter tuning is using a sophisticated search algorithm based on random forest and optimizing towards the lowest CRPS on a subset of our observations (i.e., a training set). We utilize the R-Package mlrMBO to execute this optimization mlrMBO. This approach brings one significant advantage: the search algorithm can efficiently search the considered space by repeatedly reevaluating the objective function. However, once the final set of hyperparameters is selected, it will remain constant throughout the forecasting task. Additionally, the tuning is only executed using a small subset of the dataset. This could be a problem as the chosen set of parameters may not be optimal for the rest of the dataset, especially if there are structural breaks. Hereinafter, we will refer to this approach as Bayesian fix as it fixes the hyperparameters after utilizing a Bayesian search algorithm.

The second approach uses the online function, which is included in the profoc R-Package profoc_package. It implements the proposed algorithm and an online tuning strategy for the hyperparameters. This strategy considers a random sample of all possible hyperparameter sets, and for each iteration, the combination with the lowest aggregate CRPS is chosen. In contrast to the Bayesian fix approach, we define all possible parameter sets before the learning task. However, this method dynamically selects the parameter set based on past performance, allowing for dynamic adjustments if underlying properties change. The most significant drawback of this approach is that only a random sample of the hyperparameter space is considered. However, this approach has the advantage of adjusting dynamically to changes in the data. Therefore, we will refer to this approach as Sampling Online.

It is also possible to combine both approaches. In this case, mlrMBO optimizes on a subset of the data. Afterward, online uses the parameter combinations that got proposed in the mlrMBO optimization. This has the potential to profit from efficient exploration of the hyperparameter space and the ability to adjust to changes in the data dynamically. After this, we will refer to this approach as Bayesian Online.

Application to Multivariate Probabilistic Day-Ahead Power Prices

In day-ahead electricity price forecasting, we consider the price $Y_{t,h}$ at day $t$ and product $h=1, \ldots, H$ of the day. For hourly electricity prices, we have $H=24$, and therefore $h$ is often referred to as hour, see ziel2018day. We consider the forecasts of marcjasz2022distributional, which covers the period from December 27, 2018 to December 31, 2020. These forecasts are based on German electricity market data starting in January 2015. barunik2023learning also use that data in their probabilistic forecasting study with the same design. They are hourly forecasts of eight models, i.e., neural network specifications. The forecasts are given as distributional parameters for each hour (i.e., $\boldsymbol{\mathcal{D}} = (1, \ldots, 24)$) of all 736 Days. We use those distributional forecasts for calculating quantiles on the equidistant grid of percentiles $\boldsymbol{\mathcal{P}}=(0.01,\dots,0.99)$.

The performance of combinations is mainly determined by two factors: the performance of the considered experts and the diversity between them. The first should naturally be high as an expert can only be beneficial if it provides valuable information; the latter is equally important since there is close to no benefit in combining very similar forecasts. Figure (ref) shows the correlation between the experts. We show Pearson's correlation on the lower triangle, which takes values in $[-1,1]$. In the upper triangle, we show the distance correlation. The distance correlation is a non-linear dependency measure that takes values in $[0,1]$ and characterizes stochastic independence szekely2007measuring. Unsurprisingly, we observe positive values for both dependence measures as all time series forecast the same target. However, all values are clearly below 1. This indicates diversity between experts, which is beneficial for the combination task.

The simulations of \citet*{berrisch2021crps} show superior performance for penalized smoothing compared to the basis smoothing approach. Therefore, we solely use the penalized smoothing approach for our learning task. We use 99 knots, i.e., one on each quantile.

We consider the knot placement and the other extensions discussed in Subsections (ref) and (ref). Table (ref) summarizes the considered hyperparameters. That is, we have a total of 15 tuning parameters to optimize.