Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
66,662 characters · 17 sections · 48 citation commands
Hierarchical forecasting for aggregated curves with an application to day-ahead electricity price auctions
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Aggregation, coherent forecasts, reconciliation, hierarchical time series, forecast combinations.
\spacingset{1.45}
Curves are frequently encountered structures in \color{black} various scientific disciplines, and especially in economics and finance. Prominent examples are \color{black} supply and demand curves marshall2009principles, mankiw2014principles, pindyck1995microeconomics. Other well-known examples are forwards and futures price curves hull2003options, yield curves gurkaynak2007us, Engel curves aitchison1954synthesis, banks1997quadratic, Philips curve phillips, and various types of cost curves eiteman1952shape, pindyck1995microeconomics. \color{black} The accurate estimation of these curves is important because other measures such as equilibrium values can be derived from them, and modeling a curve as a whole has the benefit of preserving additional information such as its shape or slope, which could be very useful in deriving the strategies of market participants, equilibrium values, or probabilistic forecasts, rather than modeling the equilibrium prices directly xmodel. \color{black}
\color{black} It is important to note that we refer to curves as a single observation at a point in time, as is the case in functional data analysis (FDA) shang2017grouped, hyndman2007robust, shang2011nonparametric. \color{black} \color{black} In this study, we consider finite-dimensional representations of curves given by a finite grid of $x$-values and their corresponding $y$-values, often referred to as $f(x)$. \color{black} Every curve with well-defined increments can be readily aggregated or disaggregated into marginal and cumulative values. For example, Figure (ref) shows a simulated curve and its corresponding marginal values at each step. The cumulative values can be obtained from the marginal values by cumulatively summing them up, and the marginal values can be obtained from the cumulative values by \color{black} discrete \color{black} differentiation. \color{black} Hence, from a computational point of view, a curve is simply a vector of aggregated values. \color{black} This inherent structure of aggregated curves can be interpreted as a hierarchy, and thus hierarchical time series theory can be applied to potentially improve the forecasting accuracy hyndman2018forecasting.
An hierarchical time series is a set of time series that naturally relate to each other, i.e. sum up or break down, according to a certain logic \color{black} wickramasuriya2019optimal, hyndman2018forecasting, spiliotis2021hierarchical. This property is referred to as coherency and it is generally not given for individual forecasts. Approaches that allow independently forecasted time series at certain or all levels of the hierarchy to be made coherent are referred to as reconciliation methods. \color{black}
\color{black} Notable advances have been recently made in optimal reconciliation approaches as well as in improving classic approaches, such as top-down and bottom-up approaches, in terms of their forecasting accuracy hyndman2018forecasting. Most notable is the optimal minimum trace reconciling method wickramasuriya2019optimal. Other methods have also been proposed, such as a game-theoretically optimal reconciliation approach van2015game, averaging approaches called level conditional coherent (LCC) and combined conditional coherent (CCC) point forecasts di2021forecast,hollyman2021understanding and machine-learning based reconciliation spiliotis2021hierarchical, bregere2022online, huard2020hierarchical. \color{black}
\color{black} In this study, we introduce four novel features into the field of hierarchical forecasting. First, we exploit the fact that all aggregated curves have an implicit hierarchical structure and that reconciliation approaches can be used to potentially improve their forecast accuracy. To our knowledge, this is the first time that aggregated curves have been treated as hierarchical structures. We consider various methods of constructing and deconstruction the curves, including different representations of aggregated curves. Second, we introduce a new, simple reconciliation approach called aggregated-down with similar complexity to the top-down approach, which we recommend using as a benchmark method alongside bottom-up and top-down approaches. Third, we study minimum trace optimal reconciliation approaches for aggregated curves in detail. In particular, we provide a result to show that under some assumptions, the reconciliation approach is independent of the representation of the curve. Finally, we applied all of these approaches in an empirical setting to forecast the supply and demand curves for day-ahead electricity price auctions. We conclude that forecast accuracy can be improved through hierarchical reconciliation. \color{black}
\color{black} To assess the effect of hierarchical reconciliation approaches on aggregated curves in an empirical setting we consider an application to the German day-ahead electricity markets. In this auction-based market, submitted bids are aggregated to form supply and demand curves narajewski2022optimal. The intersection of the curves yields the day-ahead electricity price. Hence, the market mechanism itself fits our idea of aggregated curves, and thus the application of hierarchical forecasting is a natural extension. \color{black}
\color{black} xmodel, ziel2018probabilistic, haben2021probabilistic forecast the day-ahead electricity price by forecasting the corresponding bids, aggregating them to form supply and demand curves, and computing the resulting intersection of these curves to yield the final price forecasts. This implicit application of the simple bottom-up reconciliation approach improved both point and probabilistic price forecasts. \color{black} Recent related papers kulakov2020x, mestre2020forecasting, soloviova2021efficient which forecast electricity supply and demand curves show similar findings.
The remainder of this paper is organized as follows. In \hyperref[sec:hierarchy]{Section 2}, we introduce the notation for the inherent hierarchical structure of aggregated curves. \color{black} We also describe different ways to construct and represent curves via aggregation and/or disaggregation. In \hyperref[sec:reconciliation]{Section 3}, we present the reconciliation approaches that we consider in this study, including the new aggregated-down approach and their simplified formulas according to the specific hierarchical structure of aggregated curves. Furthermore, we show that under some assumptions, the reconciliation approach is independent of the previously introduced curve representations. In \hyperref[sec:simulation]{Section 4}, we present a simulation study of the new aggregated-down approach where we analyze and compare with the other reconciliation methods considered in this study. \color{black} In \hyperref[sec:Application]{Section 5}, we introduce the market clearing mechanism for the German day-ahead electricity market and the model used to forecast the day-ahead supply and demand curves. We also present our empirical study, the data used, and the empirical results. We summarize the results and provide our conclusions in \hyperref[sec:conclusions]{Section 6}.
\color{black} The hierarchical structure of a curve can be represented in multiple ways, i.e. the relationship between the bottom-level or marginal values and the curve itself, or the aggregated values. The curve can be obtained from an aggregation procedure with the marginals. Alternatively, we can start with the curve and describe a disaggregation relationship to obtain the bottom values.
However, a natural method of representation that we refer to as the canonical structure of an aggregated curve is shown in Figure (ref). At the end of this section, we also discuss other representations. \color{black}
\color{black} The canonical representation of an aggregated curve at we illustrate in Figure (ref) is regarded as the standard representation. This representation also appears naturally when modeling auction curves in economics xmodel. It can be obtained by considering the bottom values together with an aggregation procedure to receive the aggregated values, but also the other way around. \color{black}
\color{black} Based on the first approach above, we introduce the following notation. Let $\boldsymbol b = ( b_{1} , \dots , b_{n} )'$ be the $n$-dimensional vector of bottom-level or marginal values with $n>1$, as shown in Figure (ref). Now, we introduce their aggregation by the vector $\boldsymbol a = ( a_{1} , \dots , a_{n} )'$, which is an $n$-dimensional vector of aggregated or cumulative values starting from the top level, where $a_i$ for $i \in \{1,\dots,n\}$ represents the curve at point $i$. Formally, this is defined as
i.e. the cumulative sum of $\boldsymbol b$. In addition, the recursive relationship holds
We note that Equation ((ref)) is simply the mathematical representation of Figure (ref). The top level $a_{n}$ is the most aggregated level.
Alternatively, we can introduce the canonical representation starting from the aggregated values $\boldsymbol a$. Then, we can obtain the bottom values $\boldsymbol b$ by disaggregation or differencing with an initial value for $b_1$:
Furthermore, we can express $\boldsymbol b$ as $\boldsymbol b = \boldsymbol D \boldsymbol a$ where $\boldsymbol D$ is an invertible $n$-dimensional quadratic matrix defined as
Consequently, the matrix representation of equation (ref) is $\boldsymbol a = \boldsymbol D_n^{-1}\boldsymbol b$ where $\boldsymbol D_n^{-1}$ is a lower-left triangular matrix with ones on the diagonal and lower triangle. \color{black} s \color{black} Using the definitions of aggregated values $\boldsymbol a$ and bottom-level values $\boldsymbol b$ we can compactly write all of the values shown in Figure (ref) as the $(2n-1)$-dimensional vector
which contains the values of $\boldsymbol a$ in inverse order except for its first value $a_1$, which matches $b_1$. Thus, $(y_1,\ldots, y_n)'$ may be regarded as aggregated values and $(y_n,\ldots, y_{2n-1})'$ as bottom-level values. \color{black}
The aggregation relationship within $\boldsymbol y$ can be represented using the $(2n-1)\times n$-dimensional summation matrix $\boldsymbol S$, which describes the hierarchical structure of the considered data such that $\boldsymbol y = \boldsymbol S \boldsymbol b$ holds. \color{black} For the considered canonical representation we have
where $\boldsymbol U_{n-1}$ is a unit anti-diagonal upper-left triangular matrix of dimension $(n-1)$ and $\boldsymbol I_n$ is an identity matrix of dimension $n$. \color{black}
\color{black} Other situations can be considered where even more parts of the hierarchical curve are available. Using the notation $b_{i:j} = \sum_{l=i}^j b_l$ the canonical setting uses $b_{1:j}$ for $j>1$ for the aggregated values, i.e. $\boldsymbol a = (b_{1:1}, b_{1:2},\ldots, b_{1:n})'$. We could also consider a setting where, e.g. $b_{2:3} = b_2+b_3$, is also available, which would enrich the hierarchical structure with this additional information. However, the resulting structures that include general $b_{i:j}$ combinations could include up to $n(n-1)/2$ values in the corresponding $\boldsymbol y$ vector. This would increase computational costs substantially, so we only study the specific structure where forecasts $\widehat{\boldsymbol y}$ (or equivalently $\widehat{\boldsymbol a}$ and $\widehat{\boldsymbol b}$) are available to the forecaster. \color{black}
\color{black} In this subsection, we discuss representations other than the canonical one for aggregated curves. We generalize the representation scheme so the canonical representation is embedded. In particular, we consider a situation where the aggregated values $\boldsymbol a$ are given and the scope involves defining a disaggregation rule to receive the bottom values $\boldsymbol b$.
We introduce $k$ as the element of $\boldsymbol a$ where from we start the disaggregation rule. In the canonical representation we disaggregate $\boldsymbol a$ starting from the first value $k=1$ of the curve, i.e. for $k=1$ it holds that $b_k=a_k$, $b_i= a_{i} - a_{i-1}$ for $i>k$ (see Equation (ref)). This starting point of aggregation/disaggregation $k$ could be changed to achieve an alternative bottom-values vector $\boldsymbol b_{[k]}$. Clearly, $\boldsymbol b = \boldsymbol b_{[1]}$ and if we start disaggregating from the end, i.e. $k=n$, we obtain $b_{[n],1}= a_n$, $b_{[n],2}= a_{n-1}-a_{n}$, and in general, $b_{[n],i}=a_{n-i+1} - a_{n-i+2}$. Thus, $b_{[n],i}$ has the opposite sign to $b_i$ for $i<n$. However, this approach can be embedded in the canonical representation, which we receive by defining $b_i = b_{[n],n-i+1}$. Thus, the theory of canonical representations can also be applied.
If we consider the disaggregation procedure with an initial value at $k$ where $1<k<n$ with the corresponding value $b_{[k],k} = a_{k}$, then we obtain a representation that is substantially different from the $\boldsymbol b_{[1]}$ and $\boldsymbol b_{[n]}$ situations. The reason for this difference is that we require two directions of aggregation, with one for bottom values larger than $k$ and the other one for smaller values. Thus, we have
We observe that the special cases of $\boldsymbol b_{[1]}$ (the canonical representation) and $\boldsymbol b_{[n]}$ can be defined using the definition (ref). For illustrative purposes and to easier understand, we provide an $n=6$-dimensional example for $k=1,3,6$ in Table (ref).
In the general $\boldsymbol b_{[k]}$ setting, the summation matrix $\boldsymbol S_{[k]}$ is given by
where $\boldsymbol L_{k-1}$ is a $(k-1)$-dimensional matrix that contains $1$ on the lower anti-diagonal. The general summation matrix $\boldsymbol S_{[k]}$ in (ref) also nests the canonical case $\boldsymbol S$ in (ref) for $k=1$. The vector $\boldsymbol y_{[k]}$ in the corresponding hierarchy that satisfies $\boldsymbol y_{[k]} = \boldsymbol S_{[k]} \boldsymbol b$ is
where we define $\boldsymbol a_{[-k]}$ as the reversed vector $\boldsymbol a$ without the $k$th element, i.e. $\boldsymbol a_{[-k]}=(a_n,\ldots, a_{k+1},a_{k-1},\ldots, a_1)'$. Clearly, for $k=1$ we have $\boldsymbol y = \boldsymbol y_{[k]}$ (see definition (ref)).
However, the different representations of the hierarchical structure are actually only formal representations and do not automatically provide different reconciled forecasts by themselves, as shown at the end of the next section. To show that this is the case, we observe that the following relations hold:
and
These relations help us define Equation (ref) in matrix form as $\boldsymbol b_{[k]} = \boldsymbol A_{[k]}\boldsymbol a$ with matrix $\boldsymbol A_{[k]}$, which yields $\boldsymbol y_{[k]} = \boldsymbol S_{[k]}\boldsymbol A_{[k]}\boldsymbol a$. For $k=1$, we have $\boldsymbol A_{[k]} = \boldsymbol D_n$ from Equation (ref). $\boldsymbol A_{[k]}$ has the form:
In addition, according to definition (ref), a matrix $\boldsymbol B_{[k]}$ exists that satisfies $\boldsymbol y = \boldsymbol B_{[k]} \boldsymbol y_{[k]}$. We note that $\boldsymbol B_{[k]}$ has the structure
Furthermore, it is easy to check that $\boldsymbol B_{[k]}$ is orthogonal, i.e. it holds $\boldsymbol B_{[k]}^{-1} = \boldsymbol B_{[k]}'$. We observe that $\boldsymbol B_{[k]}$ is a generalized permutation matrix that contains permutations and reflection components.
Finally, using $\boldsymbol y = \boldsymbol B_{[k]} \boldsymbol y_{[k]}$, $\boldsymbol y_{[k]} = \boldsymbol S_{[k]}\boldsymbol A_{[k]}\boldsymbol a$ and $\boldsymbol a = \boldsymbol D_n^{-1}\boldsymbol b$ we obtain that $\boldsymbol y = \boldsymbol B_{[k]}\boldsymbol S_{[k]}\boldsymbol A_{[k]} \boldsymbol D_n^{-1}\boldsymbol b $. Thus, it holds that $\boldsymbol S = \boldsymbol B_{[k]}\boldsymbol S_{[k]}\boldsymbol A_{[k]} \boldsymbol D_n^{-1}$ for all $k$, which is valuable for further analysis. \color{black}
In general, the hierarchical structure does not hold if each time series is forecasted individually. As stated by wickramasuriya2019optimal, we call these individual forecasts incoherent or base forecasts, and denote them by $\widehat \boldsymbol y$. \color{black} The coherent or reconciled forecasts for which the structure in Figure (ref) holds are denoted by $\widetilde{\boldsymbol y}$, and they formally satisfy: $\widetilde{\boldsymbol y} = \boldsymbol S \widetilde{\boldsymbol b}$, where $\widetilde{\boldsymbol b}$ are the bottom values in $\widetilde{\boldsymbol y}$. \color{black}
The coherency of forecasts is a desired property of aggregated curves for further analysis, e.g. for trading applications, but the main application of reconciliation is generating accurate values for the aggregated level $\boldsymbol a$. All linear reconciliation approaches for any hierarchical structure can be compactly written in matrix notation as
\color{black} The $ (2n-1) \times n$-dimensional mapping matrix $\boldsymbol P$ maps the base forecasts $\boldsymbol y$ to the bottom level, hyndman2018forecasting. $\boldsymbol P$ is different for each reconciliation approach. The product $\boldsymbol S\boldsymbol P$ is sometimes referred to as the reconciliation or projection matrix.
In the previous section, we described the canonical representation and other representations. We note that they imply the same hierarchical structure and it holds that $\boldsymbol S = \boldsymbol B_{[k]}\boldsymbol S_{[k]}\boldsymbol A_{[k]} \boldsymbol D_n^{-1}$. Furthermore, the results are actually equivalent for specific reconciliation approaches considered in the simulation study and application, i.e. they are not affected by the choice of $k$ (see Section (ref) for more details). Therefore, in the following subsections, we only consider the canonical representation. \color{black}
\color{black} Table (ref) summarizes the formulas for our three benchmark reconciliation approaches: bottom-up, top-down and aggregated-down. The first two approaches are well-established procedures, so we will not provide their details but instead we refer to previous studies (gross1990disaggregation,hyndman2018forecasting). For the top-down and aggregated down methods, we provide three options for calculating the disaggregating proportions, i.e., using the average ratio (ar), ratio of averages (ra), and forecasted values (fo). Many more options could be considered for the top-down approach (gross1990disaggregation). We describe the new aggregated-down approach in the next subsection. \color{black}
We use bu as an abbreviation for bottom-up, and tdar, tdra, tdfo for the top-down approaches using the average ratio, the ratio of averages, and forecasted values respectively.
\color{black} The aggregated-down approach is essentially a localized top-down approach, where the disaggregating proportions are calculated based on the node above and not based on the top-most aggregated level. The motivation for this approach is based on the simple assumption that proportions calculated using values that are closer in the hierarchical structure are expected to be more accurate than when using values that are further away. We provide support for this assumption based on our simulation study in the next section. \color{black}
\color{black} We denote the corresponding disaggregating proportions by $q_j$ where $\boldsymbol q = ( q_1 , \dots , q_{n} )'$ is the vector of proportions. By design (see Fig. (ref)), it holds that
Thus, $q_j$ is the disaggregation proportion for $b_j$ given $a_{n-j+1}$ at the corresponding connection of the tree. \color{black} The mapping matrix is defined as $$\boldsymbol P_{ad} =
,$$ where $\boldsymbol Q = \text{Antidiag}( \boldsymbol q)$ is a $n \times n$ anti-diagonal matrix with the elements on the anti-diagonal starting from top to bottom, i.e., $$\boldsymbol Q =
.$$
We require that $q_1=1$ every time, which corresponds to $\tilde{b}_1 = \hat{b}_1$. \color{black} Analogous to the top-down approach, the proportions can be estimated using various methods, which are summarized in Table (ref). It should be noted that the methods for calculating the average ratio and ratio of averages proportion calculation methods have similar complexity for the top-down and aggregated-down approaches. However, the formula for calculating the aggregated-down approach using the forecasted values $\widehat{\boldsymbol y}$ is much simpler than that for the top-down approach. \color{black}
We use adar, adra, adfo as abbreviations for the aggregated-down approaches using the average ratio, the ratio of averages, and forecasted values, respectively (see Table (ref)).
hyndman2018forecasting and wickramasuriya2019optimal introduced the optimal reconciliation approach by showing that the optimal mapping matrix that obtains the best, unbiased coherent forecasts is given by
where $\boldsymbol W=\text{Var}[\boldsymbol y-\widehat \boldsymbol y]$ is the variance-covariance matrix of the base forecast errors. Equation ((ref)) is the result obtained by minimizing the variance of the coherent forecasts. $\boldsymbol W$ is not known \color{black} and must be estimated. We consider multiple estimators for $\boldsymbol W$, which are listed below. \color{black} The first five estimators were proposed by hyndman2018forecasting and wickramasuriya2019optimal.
\color{black}
Next, we briefly describe the optimal reconciliation approach for other representations. In this situation, for some $k$, we have $$\widetilde{\boldsymbol P}_{[k]} = (\boldsymbol S_{[k]}'\boldsymbol W_{[k]}^{-1}\boldsymbol S_{[k]})^{-1}\boldsymbol S_{[k]}' \boldsymbol W_{[k]}^{-1}.$$ Furthermore, we can show that under mild assumptions regarding the forecast method and the reconciling matrix $\boldsymbol W_{[k]}$, the reconciliation approach preserves the hierarchical structure, which holds in the sense that the result does not depend on the choice of $k$:
The proof is presented in the \hyperref[sec:appendix]{Appendix}. We recall that $\boldsymbol B_{[k]}$ is orthogonal, and thus the assumption $\widehat{\boldsymbol y} = \boldsymbol B_{[k]} \widehat{\boldsymbol y}_{[k]}$ is satisfied if the forecasting algorithm that provides $\widehat{\boldsymbol y}$ is invariant to orthogonal transformations, e.g., this holds for linear regressions. In addition, we note that $\boldsymbol W^{-1}_{[k]} = \boldsymbol B_{[k]}'\boldsymbol W^{-1}\boldsymbol B_{[k]} $ is trivially satisfied if $\boldsymbol W = \boldsymbol I_{2n-1}$ due to orthogonality of $\boldsymbol B_{[k]}$.
\color{black}
\color{black} The new aggregated-down approach is motivated based on the simple idea of closeness, i.e. proportions calculated using values that are closer in the hierarchical structure are expected to be more accurate than those calculated when using values that are further away due to avoidance of error aggregation. To support this claim, we conducted a simulation study, which also allowed us to compare the performance with other approaches. \color{black}
\color{black} We simulated Vector Auto-Regressive VAR(1) processes for the bottom-level values. We replicated the simulation 1000 times for each combination of the parameters specified in Table (ref). \color{black}
\color{black} We also considered the case of correlated errors with the variance-covariance matrix $A= 0.3 \boldsymbol I_n + 0.7 \boldsymbol 1 \boldsymbol 1'$ in combination with the coefficient matrix $\boldsymbol \Phi = 0.7 \boldsymbol I_n$. We generated the aggregated values from the bottom-level simulations and fitted an Auto-Regressive AR(1) process without intercept for each level of the hierarchy. We calculated the forecast accuracy from 1-step-ahead forecasts for each level. Next, we computed the root mean square errors (RMSE) value of each reconciling approach considered. \color{black}
\color{black} The results are shown in Table (ref) for the setup with $\boldsymbol \Phi = 0.7 \boldsymbol I_n$ and $A= \boldsymbol I_n$. We left out the results for the average ratio and ratio of averages approaches for top-down and aggregated-down since they were always greatly inferior to those obtained using the base case. We also omit some optimal methods for brevity because these tended to produce very similar results. \color{black}
\color{black} The key finding from this simulation study is that the aggregated-down approach using forecasted values was very similar to the other methods, and even optimal ones. By contrast, the top-down approach using forecasted values was markedly inferior in the setups considered. As expected, the accuracy generally decreases as more bottom-levels ($n$) were added. However, this was disproportionally the case with the top-down approach, which was caused by extreme sensitivity to outliers because the proportions were multiplicatively connected and common factors were present across the hierarchy. We refer to this issue as the error inheritance of proportions.
This situation is clear when we consider the formulas. For illustrative purposes, we consider the simplest proportions $\hat p_{fo,n}$, and $\hat p_{fo,n-1}$: $$\hat p_{fo,n} = \frac{\hat b_{n}}{\hat a_{n-1}+ \hat b_n}, \ \ \hat p_{fo,n-1} = \frac{\hat b_{n-1}}{\hat a_{n-2}+ \hat b_{n-1}} \underbrace{\frac{\hat a_{n-1}}{ \hat a_{n-1}+ \hat b_n}}_{=1-\hat p_{fo,n}} $$ Clearly, the formulas are sensitive to outliers when the denominator $\hat a_{n-1}+ \hat b_n$ is very small, i.e. when $\hat a_{n-1} \approx -\hat b_n$. If this is the case, then all of the subsequent proportions will generally be outliers because $\hat p_{fo,n}$ is part of $\hat p_{fo,n-1}$, and so on. It should be noted that this occurs if we allow negative values. If we only consider monotonously increasing curves, then these outliers do not occur, as shown in Table (ref), where we simulated a strictly positive AR(1) process leading to errors close to the other methods. To understand the sensitivity to outliers, Table (ref) shows the results for the same setup with all outliers defined as $|\hat p_{fo,1}|>50$ removed. The errors improved for greatly tdfo, especially with increasing values of $N$. \color{black}
\color{black} This issue is essentially error aggregation and is important considering the specific hierarchical structure of aggregated curves because the hierarchy considered is completely vertical and each additional point on the curve deepens the hierarchy instead of widening it. Thus, if we allow for negative values, the outliers become larger for tdfo when the hierarchy is deeper. The aggregated-down approach is not sensitive to this distance issue since proportions are not multiplicative, but instead they are always calculated using values immediately above the hierarchy level considered, and non-monotonic curves can be easily accommodated. \color{black} The other simulation setups yielded similar results, as shown in Appendix Tables (ref), (ref), (ref), and (ref).
It is important to note that we consider the aggregated-down approach as a simple benchmark approach. It is not equivalent to more sophisticated optimal reconciliation methods. Nevertheless, aggregated-down has the advantages of being intuitive and using very simple formulas for calculating proportion, especially compared to the top-down approaches, and is programmatically simple as well as efficient. Furthermore, it is not affected by the error aggregating issue, unlike the top-down approach.
\color{black} The European day-ahead electricity market is a daily, blind auction market where electricity prices are determined for each hour of the next day. The auction closes at 12:00 CET each day until when the market participants, i.e. the buyers and sellers, can submit buy (bid) and sell (ask) orders for each hour of the next day. These bids are the volumes that buyers are willing to buy and that producers are willing to sell at certain prices. For the time range of data considered, bid volumes can be specified for prices ranging from -500 EUR/MWh up to 3000 EUR/MWh with an increment of 0.1 EUR/MWh\footnote{Since 2022-05-11 the upper price limit was increased to 4000 EUR/MWh.}. After gate closure, the volumes bid across all participating countries are aggregated and unique prices and volumes are generated for each separate market area. One of these areas is Germany-Luxembourg, which we modeled and forecasted in this study. The market clearing results as well as the 24 supply and demand curves are published shortly after gate closure at 12:00 CET. All major relevant aspects of the German day-ahead electricity market are illustrated in Figure (ref) (see petropoulos2022forecasting for more details). Only selected prices can be bid on out of a total of 35001 possible prices, so the curves have a characteristic step-like appearance due to the many zero-volume prices, i.e. prices that were not bid on xmodel. \color{black}
Figure (ref) shows the German volumes bid for each price for a selected day and hour. In total, there are 35001 prices that can be bid on for each hour, where most usually have a volume of zero, \color{black} as can be observed by the many empty intervals between the bars \color{black}. Certain prices are generally preferred by market participants, such as multiples of 10 and those in the vicinity of 40 EUR/MWh, as shown by the clusters formed around these prices in Figure (ref). The hourly bids are aggregated to produce the day-ahead supply and demand curves for each hour, as shown in Figure (ref).
The X-Model introduced by xmodel is an approach for forecasting day-ahead electricity prices as the intersection of supply and demand curves. The electricity prices are not modeled directly but instead they are obtained from the forecasts of the whole day-ahead curves as their intersection.
xmodel,ziel2018probabilistic,kulakov2020x, haben2021probabilistic modeled the day-ahead supply and demand curves by grouping the prices into price classes, thus drastically reducing the dimensionality, e.g., from 350001 prices to a few price classes (40 in our study), to make the problem computationally feasible. The prices are split into price classes by inverting the supply and demand curves at a pre-specified grid of equidistant volumes xmodel. The result is shown in Figure (ref), where the curves are approximated by ca. 20 points each. \color{black} This form of dimensionality reduction via binning can be interpreted as a special case of applying FDA where the basis functions are constants, i.e. dummy variables, that do not overlap. This works very well in our case because the original curves are not smooth. An additional advantage is computability because it would be quite costly to smoothen the curves first, model them, and then un-smooth them again for the final representation. xmodel use a simple but very effective method for reconstructing the forecasted curves from the price class approximations, i.e., by using historical proportions for each price. This method was shown to capture the shape of the original curves very well without requiring vast computational resources. However, this reconstruction step is only important for calculating the price as the intersection of the supply and demand curves and does not affect the forecasting accuracy of the curves themselves, hence we refer to the original X-model paper for more information xmodel as well as other papers that use more classical FDA approaches to implement the X-model soloviova2019modeling, soloviova2021efficient. \color{black}
Thus, the response variable is the sum of the volume within a price class, which is modeled separately for every hour and class. For $40$ price classes comprising both supply and demand, a total of $40 \times 24$ regressions must be conducted to forecast curves for all hours in the next day. xmodel modeled each price class volume marginally, and thus the supply and demand curves were generated by cumulatively summing up the forecasted values, which is equivalent to forecasting only the bottom-level values and using the bottom-up reconciliation approach described earlier. In this study, we modeled both the marginal and cumulative responses, and used them to compare and contrast the different reconciliation approaches. It should be noted that the method using marginal values will always be equivalent to the bottom-up approach.
To model the responses, we used a combination of autoregressive and external regressors. Let then $X_{S,d,h}^{(c)}$ and $X_{D,d,h}^{(c)}$ be the supply and demand volumes at day $d$ and hour $h$ of price class $c$ among the price classes generated for the supply and demand curves, respectively. These constitute the bulk of the regressor matrix because each price class volume will depend on its lags according to a specific lag structure. The external regressors are denoted by $X_{X,d,h}^{(1)}, \dots, X_{X,d,h}^{(M_X)}$ for a total of $M_X$ external regressors, and they comprise the prices for coal, gas, oil and CO$_2$ emissions (EUAs), the day-ahead prices and volumes on the previous day, as well as the day-ahead forecasts for the country-wide load, solar, onshore wind, and offshore wind production. \color{black} The day-ahead data were taken from www.epexspot.com and the forecasts from www.entsoe.eu. \color{black}
Let $M_S$ and $M_D$ be the number of price classes for the supply and demand curves, respectively. Then. all regressors can be compactly written as $$\textbf{X}_{d,h} = \left ( X_{1,d,h}, \dots, X_{M,d,h} \right )' = \left ( \left (X_{S,d,h}^{(c \in C_S)} \right ), \left ( X_{D,d,h}^{(c \in C_D)} \right ), \left ( X_{X,d,h}^{(c \in C_X)} \right ) \right )',$$ where $C_S$ and $C_D$ are the sets of price classes for the supply and demand sides, respectively, $C_X$ is the set of external regressors $C_X = \left \{ X_{X,d,h}^{(1)} 1, \dots, X_{X,d,h}^{(M_X)} \right \}$, and $M=M_S+M_D+M_X$.
To capture the weekly seasonality, we also included dummy regressors for every day of the week which we denoted by a function $W_k(d)$ that returns the day of the week $d$. The full model can be written as
$$X_{m,d,h} = \sum_{k=1}^{M} \sum_{j=1}^{24} \sum_{k \in \mathcal{I}_{m,h}(l,j)} \phi_{m,h,l,j,k} X_{l,d-k,j} + \sum_{k=2}^{7} \psi_{m,h,k} W_k(d) + \varepsilon_{m,d,h},$$
for $m \in \{1, \dots, M_S + M_D\}$ and $\mathcal{I}_{m,h}(l,j)$ represents the sets of possible lags, which we defined as \[ \mathcal{I}_{m,h}(l,j) =
. \] We fitted the models using lasso xmodel as implemented in the glmnet R package.
We conducted a day-ahead rolling window forecasting study for each day between 2019-01-01 and 2019-06-30. We used a rolling window length of $730 = 2 \times 365$ days to forecast each day. The price classes were generated only once for the first day-ahead forecast, i.e. using the 2017 and 2018 data, and they were kept constant throughout the study. All data were hourly except for coal, gas, oil, and EUA prices, which were daily. We used marginal values as regressors to forecast the marginal values $\widehat b_i$ and cumulative values as regressors to forecast the cumulative values $\widehat a_i$. We chose $\widehat a_n \leftarrow \widehat b_1$ for each day and hour, as in Figure (ref).
To measure the forecasting accuracy, we used two popular error measures comprising the mean absolute errors (MAE) defined as $$\textrm{MAE}_m^{\textrm{test}} = \frac{1}{24\cdot \#(\mathcal{D})} \sum_{d\in \mathcal{D}} \sum_{h=0}^{23} |X_{m,d,h}-\widehat X_{m,d,h}|$$ and the RMSE defined as $$\textrm{RMSE}_m^{\textrm{test}} = \frac{1}{\#(\mathcal{D})} \sum_{d\in \mathcal{D}} \sqrt{\frac{1}{24} \sum_{h=0}^{23} (X_{m,d,h}-\widehat X_{m,d,h})^2}$$ where $\mathcal{D}$ is a set containing all 181 forecasted days, $\#(\cdot)$ is a function that returns the number of elements in a set, and $\widehat X_{m,d,h}$ represents the respective forecasted value.
The averages over the MAEs and RMSEs for the overall price classes are shown in Figure (ref). The top-down and aggregated-down approaches using historical proportions yielded worse results than the simple bottom-up method. Therefore, we do not include these results but they can be inspected in the comprehensive Tables (ref) and (ref). For the supply curve, the optimal lambda reconciliation approach obtained the smallest MAE on average, approximately 33 MWh lower than that with the marginal model. For the demand curve, the optimal shrinkage reconciliation approach returned the smallest MAE on average, approximatively 97 MWh lower than that of with the marginal model. Similar results were obtained in terms of RMSEs, with differences of approx. 39 MWh and 109 MWh for the supply and demand curves respectively.
\color{black} Table (ref) shows the absolute MAE and percentage differences for some selected classes and the mean over all classes. Again, we do not include the results for the top-down and aggregated-down approaches using historical proportions for the same reasons stated for Figure (ref). When using historical values to calculate the proportions, the top-down approach was superior to aggregated-down. However, using forecasted values to calculate the proportions for the aggregated-down approach yielded better results than the bottom-up case on average. The top-down approach still yielded less accurate values than the bottom-up approach even when using forecasted values to calculate the proportions. Hence, aggregated-down was superior to bottom-up and top-down only when using forecasted values for the proportions. \color{black}
Tables (ref) and (ref) in the Appendix show the detailed MAEs and RMSEs over each price class. For the supply curve, no approach consistently yielded the lowest errors across all price classes. For example, the first three price classes had the lowest MAEs and RMSEs using the simple bottom-up approach. The cumulative and marginal models in the tables refer to simply forecasting all bottom-level values with the corresponding marginal or cumulative regressors and simply cumulatively summing them up, which corresponds to a bottom-up approach. It should be noted that in the results tables, the explicit bottom-up approach is equal to the cumulative approach per definition. The optimal weighted least squares (WLS) reconciliation approach yielded the lowest errors for the largest number of classes, followed by the shrinkage and lambda approaches. For the demand curve, the results were consistent for all classes, where the lowest errors were achieved by the optimal shrinkage reconciliation approach (option 5 among the Section \hyperref[sec:opt]{3.3}).
\color{black} In this study, we considered the hierarchical structure of aggregated curves and different representations. \color{black} We presented several reconciliation methods comprising established bottom-up, top-down, and minimum-trace optimal reconciliation approaches in an aggregated curves setting wickramasuriya2019optimal. In addition, we introduced a new aggregated-down approach with comparable methodological complexity to the bottom-up and top-down approaches. \color{black} We provided a theoretical insight that under mild assumptions regarding the forecasting and reconciling method, the reconciling result is independent of the representation of the curve. These approaches were then applied in a simulation study, and for forecasting the supply and demand curves \color{black} of the German day-ahead electricity market. In the latter case, we showed that the reconciling approaches for aggregated curves can improve forecasting accuracy compared with the standard approaches.
The results showed that there is a single reconciliation method did not outperform the others every time. However, the Section \hyperref[sec:opt]{3.3} (optimal approaches), particularly the shrinkage, WLS, and lambda approaches, obtained the best results in most cases, where they obtained considerable improvements compared with the bottom-up base case. The reconciliation method that improves the forecast by the greates amount will probably be specific to the data.
We suggest that it may be useful to consider multiple methods because even simple approaches such as Section \hyperref[sec:ad]{3.2} (aggregated-down) yielded improvements at certain points on the curve compared with the bottom-up approach, which is the current state-of-the-art method xmodel, ziel2018probabilistic, forecast2021haben. \color{black} Based on our finding that the top-down approach did not improve the forecasts on average whereas the aggregated-down approach using forecasted values to calculate proportions obtained improvements, we conclude that the aggregated-down approach can be applied as a simple benchmark method that is superior to top-down. In addition, we conclude that it is important to have access to all base forecasts to calculate the proportions because they can lead to substantial improvements in the forecasts compared with only using historical values. \color{black}
Our study could be extended further by including more recent approaches, such as machine-learning-based spiliotis2021hierarchical and conditional coherency reconciliation methods di2021forecast. Averaging or using different reconciliation approaches for each point on the aggregated curve could also be considered, especially if coherency is not necessary.