{'fraction_non_alphanumeric': 0.038428142598427154, 'fraction_numerical': 0.018085291126051253, 'mean_word_length': 4.627530364372469, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From point forecasts to multivariate probabilistic forecasts: The Schaake shuffle for day-ahead electricity price forecasting

Modeling price risks is crucial for economic decision making in energy markets. Besides the risk of a single price, the dependence structure of multiple prices is often relevant. We therefore propose a generic and easy-to-implement method for creating multivariate probabilistic forecasts based on univariate point forecasts of day-ahead electricity prices. While each univariate point forecast refers to one of the day's 24 hours, the multivariate forecast distribution models dependencies across hours. The proposed method is based on simple copula techniques and an optional time series component. We illustrate the method for five benchmark data sets recently provided by Lago et al. (2020). Furthermore, we demonstrate an example for constructing realistic prediction intervals for the weighted sum of consecutive electricity prices, as, e.g., needed for pricing individual load profiles.


Introduction
Day-ahead electricity price forecasting is a critical element in the decisionmaking of energy companies. Accordingly, an active applied research literature is concerned with developing and comparing price forecasting methods. The majority of this literature addresses point forecasts, for which a wide range of methods has been proposed. These efforts, which have recently been reviewed by Hong et al. (2020) and Lago et al. (2020), have given rise to a rich toolkit While univariate point forecasts are an important first step, they are limited in two respects. First, they do not address forecast uncertainty, which is an important concern for economic decision making. The studies surveyed by Nowotarski and Weron (2018) and Ziel and Steinert (2018) address this aspect by considering probabilistic forecasting in energy markets. While growing, the corresponding literature and range of approaches is still more limited than for point forecasting. Second, univariate forecasts are concerned with a single measurement unit (such as a given time period or a given location), whereas economic decisions often require joint forecasts for several time periods or locations. In applications, the two aspects of forecast uncertainty and multivariate dependence may well interact. For example, managing daily price risks requires information on both the uncertainty in hourly prices and information on dependencies of prices across hours. Consider the sum of the prices from two consecutive hours. If the price of the first hour is higher than predicted, the price of the second predicted hour is likely to be higher as well. Assuming independence of these prices would underestimate the uncertainty of the sum, leading to suboptimal economic decisions. Dependence modeling thus plays a vital role in uncertainty estimation. This paper considers a generic and easy-to-implement method for constructing multivariate probabilistic forecasts based on univariate point forecasts. The method uses past forecast errors to learn about forecast uncertainty and dependencies across measurement units. If necessary, we first fit simple time series models to the univariate point forecast errors, thus accounting for possible predictability in the errors' conditional mean and variance. For dependence, we apply straightforward empirical copula methods. Being based on forecast errors, our approach leverages the rich literature on univariate point forecasting (thus avoiding to re-invent the wheel) while at the same time addressing forecast uncertainty and multivariate dependence. It is inspired by the literature on post-processing multivariate ensemble forecasts in meteorology (Clark et al., 2004;Schefzik et al., 2013;Vannitsem et al., 2021). Here, the point forecast errors of numerical weather prediction models are used to learn about forecast uncertainty and (possibly) about multivariate dependencies across variables, locations or time points. Related post-processing approaches have been considered in various other areas of application, such as macroeconomics (Krüger and Nolte, 2016;Clark et al., 2020).
Post-processing seems particularly promising for forecasting energy prices or demand, given the rich literature on point forecasting and the wide availability of data for learning structures in past point forecasting errors. Nevertheless, the literature on post-processing of energy forecasts is relatively sparse. In a univariate context, studies such as Marcjasz et al. (2020) and Kath and Ziel (2021) discuss approaches to estimate forecast uncertainty from past point forecast errors. Phipps et al. (2020a,b) and Ludwig et al. (2020) consider post-processing univariate or multivariate weather forecasts, in a situation where these are used as an input to energy forecasting models. Muniain and Ziel (2020) deal with bivariate probabilististic price forecasting, using a residual based approach that is conceptually similar to post-processing. Janke and Steinke (2020) apply an implicit generative model to generate a multivariate forecast distribution for energy prices from an ensemble of univariate point forecasts. Furthermore, Chai et al. (2019) employ a Gaussian copula to generate scenarios on the basis of an ensemble from extreme learning machines. Compared to the latter two studies, our use of empirical copula methods appears considerably simpler to implement and more easily comparable to the large meteorological literature. Other uses of copula methods to model multivariate dependence in energy contexts include Toubeau et al. (2019), Pinson and Girard (2012), Manner et al. (2016) and Pircalabu and Benth (2017).
We next present a more specific summary of our approach, full details of which are given in Section 3. We consider a vector of 24 day-ahead point forecasts, each of which refers to one hour h of the following day t (see Lago et al. (2020) for a description of day-ahead markets). From now on, we stick to this hourly example to simplify the presentation, but quarter-hourly or minute-byminute forecasts could be handled analogously. The prediction is subject to an error defined by t,h = y t,h −ŷ t,h , where y t,h denotes the true price andŷ t,h the prediction for each day t and hour h = 1, . . . , 24. Conditional on the information set F t−1 , available at day t − 1, t,h is a random variable with unknown, conditional distribution function that we model as G t,h (x) := F h ( x−µ t,h σ t,h ). Here, F h reflects an hour specific shape of the distribution while σ t,h and µ t,h are scale and location parameters.
In what we call the Error Learning Phase, the error distributions as well as the scale and location parameters are estimated. We allow a potentially timevarying conditional mean µ t,h = E[ t,h |F t−1 ] and variance σ 2 t,h = V[ t,h |F t−1 ]. The nonzero conditional mean µ t,h is motivated by possible misspecification of the point forecasting model, leading to predictable forecast errors. Time variation in the conditional variance σ 2 t,h accommodates for differences in price uncertainty over time. If predictable structure in the data is apparent, µ t,h and σ 2 t,h can be tracked by time series models, e.g., by estimating the models' parameters in rolling windows of past forecast errors, separately for each hour. Alternatively, with no predictable structure apparent, we simply set G t,h (x) := F h (x), so that the error distribution remains constant over time. We call the latter the 'raw'-error approach. In the Dependence Learning Phase, the procedure aims to learn the dependence of the point forecast errors. Therefore, we estimate the joint distribution of the 24-dimensional vector of day-ahead forecasting errors t,h . We disentangle univariate error distributions and their dependence by inferring the copula from the corresponding uniformly standardized values, i.e., û t,h :=Ĝ t,h ( t,h ), whereĜ t,h has been obtained from the Error Learning Phase. Finally, the learned error distribution and dependence structure are combined in the Forecasting Phase to construct a forecast distribution consisting of m simulated time paths of day-ahead electricity prices, with each path consisting of 24 hourly prices. This simulated distribution encodes both forecast uncertainty for each individual hour, and the temporal dependence across the day's 24 hours. Among others, it can be used to assess the likelihood of price-related events (such as the event of a price spike), and to measure the uncertainty of summary random variables (such as the sum of hourly prices). We consider an example for the price of a Standard Load Profile, reflecting a typical daily pattern of energy consumption, and show that the proposed method is able to capture the dependence structure of hourly prices, leading to a realistic assessment of price risks. Furthermore, we highlight the importance of the dependence structure in a multivariate probabilistic forecast by comparing our results to their naïve counterpart that assumes temporal independence of the prices.
The remainder of the paper is structured as follows. Section 2 presents the data sets of Lago et al. (2020) and day-ahead load data that we consider in addition. In Section 3 we explain our proposed method in detail and give a toy example. Section 4 provides a case study for multivariate probabilistic electricity price forecasting in five markets and an example for the price of a Standard Load Profile in Germany. Further, we apply our method on a highly seasonal load forecasting time series to show the benefit of the time series component. Last, Section 5 summarizes our work and points towards future research.

Data Sets
The method we consider is based on a data set of past point forecasts and corresponding realizations. In our illustrations, we utilize point forecasts of energy prices and the corresponding realized prices. As discussed in the introduction, a wide range of point forecasting models and data sources are available for this purpose. Here we consider the five benchmark data sets recently made available by Lago et al. (2020). The data covers five day-ahead electricity markets, each spanning a history of six years. In addition to observations of day-ahead prices, two time series of influential exogenous variables are included, which differ for each market. To illustrate the data, Figure 1 displays the realized day-ahead electricity prices in EUR per MW from 4 January 2012 to 31 December 2017 in the data set for the German electricity market (EPEX-DE). We generate point forecasts for the five data sets from the implementation of the LEAR model by Lago et al. (2020) with a two year calibration window. 1 The LEAR model is a parameter-rich ARX model (that is, an autoregressive model with a large number of 'exogenous' regressors or features) estimated using LASSO for feature selection. In total, 247 features are considered, including day-ahead prices, day-ahead forecasts of the two exogenous variables, and historical day-ahead forecasts of the exogenous variables, all stemming from previous days and weeks. Additionally, a dummy variable for the day of the week is included. Finally, all variables are preprocessed with the arc hyperbolic sine (asinh) variance stabilizing transformation. Details on the data and the LEAR model can be found in Lago et al. (2020, Section 4.2).  In our forecast evaluation analysis, we further consider a Standard Load Profile (SLP) that provides a practically relevant example of intra-day seasonality in energy demand. We use this profile to construct a single daily energy price from 24 hourly prices. This setup allows us to assess whether our multivariate probabilistic predictions of hourly prices enable a realistic estimation of the uncertainty in daily prices. In practice, German energy suppliers use SLPs to model the consumption patterns of electricity cus-tomers without registered power metering. The SLPs replace the non-existent load profile of these customers with a forecast of electricity consumption during every quarter-hour. We utilize SLP G0 (see Figure 2), which represents the average of all industrial SLPs in Germany and is provided by the German Association of Energy and Water Industries (Bundesverband der Energie-und Wasserwirtschaft e.V.).
We also demonstrate the application of our method for the related task of load forecasting. To this end, we use day-ahead load data from the entso-e transparency platform (https://transparency.entsoe.eu/) for the German market from 2016 to 2020. Load forecasting is essential for planning and operation of energy suppliers, system operators and other market participants (Weron, 2013). Figure 3 shows the time series of errors that results from a 'black box' forecast provided on the platform. Upon closer inspection, the error series displays strong seasonal patterns and therefore calls for appropriate time series modeling techniques for postprocessing.

Sample Choices for Forecasting Scheme
Training of our model is based on two phases. First, the univariate error distribution is learned (Error Learning Phase). Second, the joint behavior of realized and standardized errors is modeled by their copula (Dependence Learning Phase). Both phases require calibration windows for learning and we implement the proposed method using backward-looking rolling windows. The windows should be long enough to provide stable estimates but also short enough to adapt quickly to changing market developments. Here, we use about one year (364 days) of raw errors t,h for the Error Learning Phase. Then, the raw errors are standardized with the help of an appropriate time series model, and the most recent 90 standardized errors are used to infer the distribution of standardized errors F h . Alternatively, when the raw errors do not show temporal structures, the 24 distributions F h , h = 1, . . . , 24 are inferred directly from the most recent 90 raw errors without filtering them through a time series model. We call the latter 'raw-error' approach. Finally, the Dependence Learning Phase estimates multivariate dependence by computing the copula between the 24 distributions F h based on a windows of the last 90 filtered or raw errors, respectively. While the trade-off between too short or too long window sizes needs to be resolved individually for each use case, our experiments indicate that the lengths proposed above are a good baseline setting for hourly prices.
After the next day's distribution of observations is predicted, the windows are rolled over by one day and the estimates of the error distributions and dependency structures are updated. The initial rolling windows used for each data set are shown in Table 1 along with the overall resulting evaluation periods.  Table 1: Sample choices in the empirical analysis. We show the initial backward looking windows for the learning phases which are then rolled over until the end of the evaluation periods (rightmost column).

Methods
The method we propose is inspired by Clark et al. (2004) and Schefzik et al. (2013) who consider multivariate dependencies in weather forecasting models, e.g., across locations or across different weather variables. Intuitively, we learn the univariate forecast distributions of 24 hours of the day (Error Learning Phase) and their dependence structure (Dependence Learning Phase) in history and use this information to construct multivariate forecasts in the future (Forecasting Phase). We next detail the three steps of the method. An algorithmic description is available in the appendix (Algorithms 1 and 2).

Error Learning Phase
In the Error Learning Phase, we estimate the conditional distributions where the indices represent the day (t) and hour (h). The conditional mean and variance of the error, µ t,h and σ 2 t,h , can be estimated via a time series model. This is recommended if, e.g., the point prediction has high autocorrelation, seasonal effects or heterogeneous variance patterns. The estimatesμ t,h andσ t,h then result from the forecast of the time series model. The required degree of complexity of the model partly depends on the point forecasting model that generated this sequence of forecast errors. In an idealized setup, one would not require a model for the expected forecast error µ t,h , as the latter is equal to zero if the forecast model is correctly specified (see e.g. Pesaran and Weale, 2006, Section 2.3). However, due to complex seasonality and persistence patterns, this goal is hard to achieve for practical energy price forecasts (see e.g. Maciejowska et al., 2021). Time variation in σ 2 t,h reflects different market phases and seems unrelated to the accuracy of the point forecasting model. 2 If no time series model is applied, we set G t,h (x) = F h (x), i.e., we assume that the distribution of y t,h |F t−1 does not change over time.
To complete the model specification, we must estimate or fix the distribution . We consider two versions: A parametric variant, where we suppose thatẑ t,h follows a standard normal distribution, and a non-parametric variant, where we estimate F h via the empirical distribution function. We hence set where z ∈ R, and n is the length of the error learning phase.

Dependence Learning Phase
The Dependence Learning Phase aims to learn the joint distribution of all t,h . Since we have already estimated the marginal distributions G t,h (x) of each hour h, we next estimate their dependence structure using copula techniques.
To motivate our use of copulas, note that we need the joint distribution of the 24 day ahead errors t,h of the point forecast model for h = 1, . . . , 24. This joint distribution is represented by the multivariate distribution function with t = ( t,1 , . . . , t,24 ), where θ = (θ 1 , . . . , θ 24 ) ∈ R 24 is a vector of threshold values. For example, for a vector of zero thresholds, θ 0 = (0, . . . , 0), G t (θ 0 ) yields the probability that all 24 forecast errors are negative (or zero). In general, the distribution function G t is a complicated object that depends on both the marginal distribution of each element t,h as well as the dependence structure across the 24 elements. However, copula methods allow us to specify the 24 elements' marginal distributions and their dependence separately. Nelson (2006), Joe (2014), and Durante and Sempi (2015) provide textbook introductions to copulas. By Sklar's theorem, there exists a copula function C : where the inputs are the elements' marginal CDFs, (Sklar, 1959). That is, the copula C constructs a multivariate CDF from 24 marginal CDFs. A simple example of a copula is given in the case of independence, where C(G t,1 (θ 1 ), . . . , G t,24 (θ 24 )) = 24 h=1 G t,h (θ h ), i.e., the copula function simply multiplies all of its arguments. In practice, independence of the hourly prices is clearly unrealistic and more flexible models of dependence are required. We can achieve this by choosing and estimating an appropriate copula function C. Both parametric and non-parametric methods are available for this purpose. In the parametric case, we require additional parameters that determine the function C. As a prominent example, the Gaussian copula sets C(G t,1 (θ 1 ), . . . , G t,24 (θ 24 )) = Φ Σ (Φ −1 (G t,1 (θ 1 )), . . . , Φ −1 (G t,24 (θ 24 ))), (4) where Φ is the CDF of the univariate standard normal distribution, and Φ Σ is the multivariate CDF of a Gaussian random vector with mean zero and correlation matrix Σ. The latter matrix has 24 × 23/2 = 276 unique elements, so that the Gaussian copula allows for considerable flexibility in modeling dependence, at the cost of parameter estimation uncertainty. In order to estimate the Gaussian copula, we use the rank correlation estimator (see, e.g., Genest and Favre 2007).
In contrast to parametric copulas, non-parametric copulas do not assume a specific functional form for C. The use of nonparametric copulas in the context of multivariate forecast preprocessing was introduced by Clark et al. (2004), who named the reordering idea after Dr. J. Schaake, member of the National Weather Service Office of Hydrologic Development. Schefzik et al. (2013) established the connection to empirical copulas. Our exposition of non-parametric copulas loosely follows Schefzik et al. (2013, Section 3.3). Consider a training sample covering m days of past point forecast errors for every hour h, i.e., {(ε t,1 , . . . , ε t,24 ) : t = 1, . . . , m}.
We standardize these forecast errors by considering the corresponding quantile Assume for simplicity that there are no ties, i.e., the training sample contains m×24 unique values ofû t,h . We denote the rank ofû t,h withinû 1,h , . . . ,û m,h by rk(û t,h ). That is, for the day t at which the smallest quantile level is observed, we have rk(û t,h ) = 1. Similarly, the day t with the largest quantile level yields rk(û t,h ) = m. The empirical copula based on a training sample of length m is then given by for integers 0 ≤ i 1 , . . . , i 24 ≤ m. For example, suppose that i 1 = . . . = i 24 = 10. In this case, the above function yields the empirical frequency of the event that all 24 quantile levels are simultaneously smaller or equal than their tenth smallest observed value. Clearly, setting i 1 = . . . = i 24 = 0 yields a function value of zero, and setting i 1 = . . . = i 24 = m yields a function value of one.

Forecasting Phase
We next describe how to draw multivariate forecasts for the next day t. From the previous two steps, we have obtained estimates of the univariate distributions G t,h (x) and the copula of standardized residualsẑ t,h . We now combine these two pieces to construct the desired multivariate distribution.
Specifically, we first construct a probabilistic univariate forecast distribution of size m, for each hour h in the next day, given by predicted quantiles of forecast error for i = 1, . . . , m, whereF −1 h is the inverse cumulative distribution function of the standardized residuals. In the case of parametric (Gaussian) marginal distributions, we have thatF −1 h = Φ −1 is the quantile function of the standard normal CDF. Thus, Equation (7) generates a stylized sample by computing m equally spaced quantiles of the forecast distribution. In the 'raw-error' case without a time series model, we simply setμ t,h ≡ 0 andσ t,h ≡ 1 in the formula above.
Up to now, the 24 ensembles represent the 24 marginal forecast densities for the hours, but lack the correct dependence structure. Therefore, the resulting m univariate ensembles per hour have to be paired according to the learned copula to encode the desired dependence structure and a discrete representation of the copula is needed. Such a representation is given by a suitable rank matrixR, containing the pairing of univariate ensembles. For the parametric copula approach, we derive the rank matrix from a random sample of size m from the fitted parametric copula model. I.e., we assign ranks to m random samples from the learned copula in each dimension starting with one for the smallest value. For the non-parametric approach, the rank-matrixR is calculated from the points of the empirical copula in the Dependence Learning Phase (see Equation 6). We then sort the univariate ensemble forecasts for every hour according to the rank matrixR to adopt the dependence structure given by the copula. We illustrate this approach in Example 3.1 which complements the meteorological example of Clark et al. (2004, Section 3b).
Example 3.1. To illustrate the method (and in particular, the Forecasting Phase), we next consider a toy example based on a Dependence Learning Phase of length m = 7, and four (instead of 24) hours of the day. First, Table 2 illustrates the four univariate forecast distributions. Equation (7) together with m = 7 implies that we represent each distribution by a set of quantiles at levels (1/8, . . . , 7/8). The differences in the four univariate distributions (e.g., the lower median price forecast at 00:00 compared to 12:00) reflect intra-daily seasonality in energy prices. This seasonality is typically already reflected by the point forecasting model that we consider as an input to our method.   Table 3 illustrates an estimated rank matrixR. For example, the table's first row shows the ranks for the first day of the Dependence Learning Phase. Prices were fairly low on that day: Among others, the price at 00:00 was the lowest recorded in the training sample, as compared to the prices at 00:00 for the other m − 1 = 6 days in the training sample. The table also shows strong positive dependence of the prices across hours, as reflected by positive correlation of the ranks.
Finally, Table 4 shows the multivariate forecast distribution that results from combining the univariate distributions in Table 2 with the rank structure from  Table 3. For example, the first draw of the multivariate forecast distribution consists of the lowest quantile for the price at 00:00 (given by 6.1), the secondlowest quantile for the price at 06:00 (given by 31.6), the lowest quantile for the price at 12:00 (given by 27.2), and the second-lowest quantile for the price at 18:00 (given by 36.5). 00:00 06:00 12:00 18:00 t = 1 1 2 1 2 t = 2 4 3 3 5 t = 3 5 4 7 7 t = 4 2 1 2 1 t = 5 3 5 5 6 t = 6 7 7 6 4 t = 7 6 6 4 3  Observe that the number of elements is the same in Tables 2, 3 and 4. At the same time, the interpretation differs across the three tables: While Table 2 represents quantile levels per time point in each row, rows in Table 3 represent ranks, and forecast draws in Table 4. This observation mirrors the simple yet clever construction behind the Schaake shuffle: Based on a dependence training sample of a given size (here, m = 7), it constructs a multivariate forecast sample of the same size. As a consequence, larger training samples yield larger and potentially more informative forecast distributions. On the other hand, short samples are quicker to adapt to possible structural breaks in dependence patters. In our empirical analysis below, we consider m = 90.

Forecast Evaluation
For the evaluation of our forecast distributions, we follow the principle of , according to which a probabilistic forecast should maximize sharpness subject to calibration. Calibration means that, ideally, each observation should resemble a random draw from the predictive distribution. On the other hand, sharpness requires that the predictive distribution be as concentrated as possible. Proper scoring rules (Gneiting and Katzfuss, 2014; assess sharpness and calibration simultaneously by assigning numerical predictive performance measures. A scoring rule is a function S(F, y) → R, where F is the forecast distribution and y the observed outcome. We consider scoring rules in negative orientation, i.e., a smaller score value indicates a better forecast. A scoring rule is called proper if stating the true forecast distribution yields the best expected score. We consider the Continuous Ranked Probability Score (CRPS) to evaluate the marginal forecast distributions. Specifically, we evaluate the forecast distribution for each of the 24 horizons, and then compute the average score. Furthermore, we use the Energy Score (ES) to evaluate the multivariate forecast distribution. We further apply the Diebold and Mariano (1995) test to assess the statistical significance of differences in forecast performance. Finally, we check the calibration of multivariate forecasts graphically with Average Rank Histograms (Thorarinsdottir et al., 2016). We provide details on these evaluation methods in Appendix A.

Electricity Price Forecasting
Our probabilistic forecasts are based on the rolling window scheme explained in Section 2.2. We examine several variants of the proposed method to investigate the impact of different modeling choices on forecast performance. Additionally, we compare all settings to their simplified (independence) counterparts ignoring the dependence structure.

Settings
The variants we consider differ in their specification of the margin (standardizing forecast errors via a time series model versus using raw forecast errors; using a non-parametric versus Gaussian marginal distribution) and the copula (non-parametric versus Gaussian). For ease of presentation, we do not consider all possible configurations but focus on the settings listed in Table 5: The Schaake-NP variant uses non-parametric margins and a non-parametric copula, as well as an AR(1)-GARCH(1,1) time series model for error standardization. The Schaake-P variant uses Gaussian margins, a Gaussian copula, and an AR(1)-GARCH(1,1) for error standardization. Last, the Schaake-Raw variant is based on raw forecast errors, i.e., without time series standardization, but is otherwise identical to the Shaake-NP variant. As a simple evaluation of copula modeling performance, we compare each variant with its independence counterpart that uses the same marginal distributions of forecast errors for each hour h, but assumes them to be independent across hours. The AR-GARCH calculation is done using the 'rugarch' package (Ghalanos, 2020) in R (R Core Team, 2013), while everything else is implemented in Python 3 (Van Rossum and Drake, 2009). Figures 4 (a), 4 (b) and 4 (c) provide illustrative calibration checks for the EPEX-DE data set. The univariate verification rank histograms for all 24 hours are displayed on the diagonal and the bivariate scatterplots of realized quantiles  on the off-diagonals. Each point in the figure corresponds to one day in the forecasting phase. Figure 4 (a) displays the result for the Schaake-NP, while panel (b) displays the case of its parametric counterpart Schaake-P, and panel (c) refers to the case of raw errors (Schaake-Raw). The off-diagonal elements in Figures 4 (a) to (c) indicate a decreasing dependence for hours that are further apart, which seems plausible. The verification rank histograms (diagonal elements) suggest good calibration of the univariate forecast distribution for the non-parametric case (Figure 4 a), as indicated by uniform histograms. Univariate calibration seems to be slightly worse in the parametric case (diagonal elements of Figure 4 (b), where most of the histograms are hump-shaped. The margins of the non-parametric raw-error approach in Figure 4 (c) also seem to be well calibrated. Thus, the results suggest better calibration for the two non-parametric versions. Furthermore, satisfactory calibration of the raw error approach suggests that error post-processing via a time series model is not of much importance in this case study. Figure 5 shows the Average Rank Histograms for various forecast distributions. The nonparametric Schaake methods in panels (a) and (c) display slight signs of miscalibration (i.e., non-uniform histograms) which however seem rather unsystematic. The histogram for the parametric case in panel (b) is humpshaped, which indicates under-dispersion and is in line with our observation in Figure 4 (b). The three naïve independence counterparts in panels (d)-(f) clearly lack calibration, as indicated by distinct U-shaped patterns of the Average Rank Histograms.   Finally, Table 6 summarises the Energy Score and CRPS for all settings and the corresponding baselines in the five markets. We naturally observe the same values for the copula based approaches and their independence counterparts for the CRPS as they share the same marginal distributions. In terms of the Energy Score, the three copula based approaches consistently outperform their independence counterparts. Furthermore, we note a better performance for both non-parametric versions (Schaake-NP/Schaake-Raw) than for the parametric error distributions (Schaake-P). This result is in line with the lack of calibration of the Schaake-P method documented earlier. Remarkably, the performance of the raw error approach (Schaake-Raw) is not systematically worse than that of the time series based non-parametric variant (Schaake-NP). For the DE, PJM and FR data sets, both variants perform similarly, and differences in score performance (Diebold-Mariano tests) are typically insignificant at the 5% level. For the BE data set, which is challenging due to several positive spikes in prices (Lago et al., 2020, Section 3.4), Schaake-Raw outperforms Schaake-NP. This may be due to increased robustness of Schaake-NP in an unstable data environment. On the contrary, in the Nord Pole data, Schaake-NP outperforms Schaake-Raw. Overall, the competitive performance of the raw error approach indicates that there is little need of systematic error postprocessing. This suggests that the forecast models provided by Lago et al. (2020) yield a good fit to the data, leaving little unmodeled heterogeneity in the models' forecast errors. The need for error postprocessing may be more pronounced in a setup with less sophisticated forecasting models. We consider such a case in the context of load forecasting below.

Load Profile Price Prediction Intervals
We illustrate the importance of the correct dependence structure in forecasting with the example of the Standard Load Profile (SLP) G0 from the Bundesverband der Energie-und Wasserwirtschaft e.V.. To do so, we compare daily price forecasts for the SLP G0 from the Schaake-NP and Schaake-Raw variants with their independence counterparts over the time period from 04/01/2016 to 31/12/2017. We average the quarter-hourly SLP across seasons and aggregate it to full hours. The true price for the SLP of each day is then calculated with realized prices from the EPEX-DE dataset. Finally, we evaluate how often this price realizes within the 93.33% quantile of the respective forecast distribution. The positive dependence between consecutive hours leads to wider prediction 2 0 1 7 -0 9 -2 9 2 0 1 7 -1 0 -0 9 2 0 1 7 -1 0 -1 9 2 0 1 7 -1 0 -2 9 2 0 1 7 -1 1 -0 8 2 0 1 7 -1 1 -1 8 2 0 1 7 -1 1 -2 8 2 0 1 7 -1 2 -0 8 2 0 1 7 -1 2 -1 8 2 0 1 7 -1 2 -2 8 intervals for the SLP. Figure 6 illustrates this aspect for a three-month subsample of the data (for better visibility), whereas Figure C.7 in the Appendix presents the entire sample period.  a simple alternative method, one could compute a prediction interval from past price realizations. In the present case, using this method with a sample of 90 past prices (in line with the choice for the Dependence Learning Phase) also yields a satisfactory coverage rate of 88.95%. Note, however, that the simple method is specific to the particular load profile considered here. By contrast, the Schaake methods can be adapted to any given load profile.

Forecasting intraday load data
We finally give a short example for other operational scenarios where our proposed method can be applied. Instead of hourly prices, this example uses intra-day load data of the entso-e transparency platform for Germany from 2016 to 2020 (https://transparency.entsoe.eu/dashboard/show). We average the day-ahead quarter-hourly load forecasts ([6.1.B]) and actual loads ([6.1.A]) over each hour, hence resulting in 24 hourly forecast error observations, in line with the setup of the price data considered earlier. We choose this data set to substantiate the use of time series models as an optional building block within the proposed method. In contrast to the point forecast errors for the price data above, point forecast errors for the load data contain time series structure with a marked seasonal component (see, e.g., the analysis in Maciejowska et al. (2021)). Without claiming to be optimal, we tackle this stylized fact by using a SARIMA(1,0,0)(1,0,0,7) time series with one AR parameter and one seasonal AR parameter with a lag of one week, i.e., 7 days for each hourly error time series (Schaake-NP). See Hyndman and Athanasopoulos 2021, Section 9.9 for an introduction to SARIMA models. We display the results for the Energy Score and CRPS for the years 2019 and 2020 in Table 8  The table indicates that raw error approaches are clearly inferior to SARIMA based error postprocessing. Thus, in contrast to the case of price data considered abover, there appears to be a clear need for error postprocessing in this example. We conjecture that this result is driven by the lower sophistication of the (unknown) load forecasting model, as compared to the forecasting models provided by Lago et al. (2020) that we considered for the energy price data.

Conclusion
This paper proposes a postprocessing method to create multivariate forecast distributions for day-ahead electricity prices out of any point forecasting model. The method is motivated by two main aspects: First, many sophisticated point forecasting models have been introduced in the literature, so that we may take the availability of a 'good' point forecasting model as given. Second, many economic decision require probabilistic multivariate forecasts which are hardly available at present.
Our method exploits several time series of univariate point forecast errors and creates a multivariate forecast distribution which inherits their dependence structure. In a case study on energy price forecasting based on benchmark models of Lago et al. (2020), a simple raw-error variant of the method performs well. This result indicates that the point forecast errors contain little exploitable structure, so that using a time series model for these errors is not necessary. In an additional case study on load forecasting, the point forecast errors contain pronounced seasonal patterns, and removing these patterns via a time series model leads to clear gains in forecasting performance. Throughout our empirical analysis, simple nonparametric techniques (for estimating the marginal distribution and the copula) outperform parametric ones based on normality assumptions.

Appendix A. Algorithmic Description of Proposed Method
Here we present algorithmic descriptions of the proposed method, in the version with a time series model (Algorithm 1) and in the raw-error version without a time series model (Algorithm 2). For ease of presentation, we omit time indices in this section; for example, we denote the price for hour h of day t by y h instead of y t,h .
Algorithm 1: Schaake shuffle with time series model (Schaake-NP/Schaake-P) if parametric margins==True then use specified distribution, e.g., N (0, 1) and its inverse CDF Φ −1 : for each ensemble member i do y h,ens (i) = (ŷ h +μ h ) + Φ −1 (i/(m + 1)) ×σ h end else use empirical CDFF h and its inverse CDFF −1 h : for each ensemble member i do y h,ens (i) = (ŷ h +μ h ) +F −1 h (i/(m + 1)) ×σ h end end end if parametric dependence==True then Fit specified parametric copula model to C Sample m times from parametric copula model of C and create rank-matrix else Create rank matrix out of saved values from C end Pair up univariate forecast ensembles according to the rank-matrix of C. end univariate forecast distribution F (represented by its cumulative distribution function, or CDF) and an outcome y, the CRPS is defined by where Y and Y are independent random variables distributed according to F , and 1(x ≥ y) denotes the indicator function which equals one if x ≥ y and zero otherwise. We aggregate the T × 24 score values (for each day and hour) by computing their mean.

Appendix B.3. Energy Score
The multivariate extension of the CRPS, the Energy Score (ES). For a multivariate forecast distribution F , the ES is given by For an ensemble forecast F ens with membersŷ ( see Gneiting et al. (2008). We average the ES across the T days of the evaluation sample.

Appendix B.4. Diebold-Mariano-Test
We apply the Diebold and Mariano (1995) test for testing whether differences in forecast performance are statistically significant. The test considers the null hypothesis that E(S 1,t − S 2,t ) = 0, where S 1t and S 2t denote the scores of forecasting methods 1 and 2 on day t. In our case, the scores correspond to either the CRPS or the Energy Score. The Diebold-Mariano test is a t-test with test statistic given byδ 12 σ(δ 12 ) whereδ 12 = T −1 T t=1 (S 1,t − S 2,t ) is the mean score difference, andσ(δ 12 ) is its estimated standard deviation. The test statistic is standard normally distributed under the null hypothesis. Since we focus on one-day ahead forecasts, we setσ(δ 12 ) equal to the sample standard deviation.