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.
105,456 characters · 13 sections · 79 citation commands
Ensemble Forecasting for Intraday Electricity Prices: Simulating Trajectories
\setcitestyle{square,sort&compress}
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} electricity price forecasting, power markets, intraday market, continuous-trade markets, XBID, ensemble forecasting, probabilistic forecasting, short-term forecasting, trajectories, generalized additive models, lasso, logistic regression, zero-inflated distribution, scenario simulation
\spacingset{1.45}
Intraday continuous electricity markets gain on importance every day goodarzi2019impact. Their primary purpose is to handle the uncertainty in electricity generation and load arisen since the day-ahead markets kath2018value. A number of events can cause the uncertainty, e.g. unexpected power plant outage or changing weather conditions. The latter one is the result of the global trend of investing in weather-dependent renewable power sources and is a subject of modelling and forecasting maciejowska2020assessing. The need of intraday continuous trading is fulfilled by the power exchanges and transmission system operators (TSO) Viehmann2017. They allow the market participants to trade the energy continuously up to 5 minutes before the delivery, e.g. in France or Germany, and to trade it cross-border, e.g. using the cross-border intraday (XBID) market kath2019modeling. Even though there is a clear evidence of the importance of this kind of markets, the researchers do not investigate them in terms of forecasting as willingly as the day-ahead market.
The day-ahead market is the main electricity spot market with a long history of research on electricity price forecasting weron2014electricity. Recent studies on the electricity price forecasting (EPF) in day-ahead markets consider i.a. the probabilistic forecasting and forecasting combination. nowotarski2018recent present a review of probabilistic EPF and muniain2020probabilistic use it to simulate peak and off-peak prices. marcjasz2018selection combine point forecasts achieved using different calibration windows while uniejewski2019importance and serafin2019averaging do it for probabilistic forecasts. A very big part of the recent EPF literature are also hybrid models yang2017electricity, wang2017multi, zhang2020adaptive and neural networks xiao2017research,bento2018bat, keles2016extended. Also the market integration plays an important role in price formation in both day-ahead and intraday markets what is elaborated by lago2018forecasting and kath2019modeling.
The role of the intraday markets in the balancing of electricity systems was emphasized and explained by ocker2017german and koch2019short on the basis of the German electricity market. They observed that the introduction of the intraday continuous market in Germany partially led to a substantial decrease in the demand for balancing energy while the wind and solar energy generation increased. karanfil2017role clarify the reason for the spread between day-ahead and intraday prices in Denmark, while maciejowska2019day forecast the price spread between the day-ahead and intraday markets based on the Polish and German data. The continuity of the intraday market has encouraged the researchers to investigate the transaction arrival process narajewski2019estimation, bidding behaviour kiesel2017econometric, graf2020modeling, rintamaki2020strategic and optimal trading strategies aid2016optimal, ayon2017aggregators, glas2020intraday. The impact of fundamental regressors on the price formation in the intraday market was examined by pape2016fundamentals, gurtler2018effect, and kremer2019fundamental.
The literature on the EPF in the intraday markets is not that broad as in the day-ahead markets or as the one regarding other aspects of the intraday markets. monteiro2016short and andrade2017probabilistic conducted the EPF for the Iberian intraday market, however it is not a continuous market, and thus their studies are more similar to these on day-ahead markets. uniejewski2019understanding, narajewski2019econometric and marcjasz2020beating performed the EPF in the German Intraday Continuous Market, while oksuz2019neural in the Turkish Intraday market. An outcome of the second one was an indication of the weak-form efficiency of the investigated market. This was partially confirmed by janke2019forecasting, who forecasted the distribution of prices during the last three hours of trading and concluded that forecasting of the central quantiles yields marginal improvement to the naive benchmark. However, marcjasz2020beating managed to outperform the most recent price by using an ensemble of it and a lasso-estimated model.
The only four papers on EPF in the German intraday market considered the ID$_3$-Price (a volume-weighted average price of transactions in the last three hours before delivery) as the most important price index in the German intraday market and conducted forecasting of it. This paper focuses on the ID$_3$ index as well, but not directly. Instead of forecasting its price we simulate the paths of 5-minute volume weighted average price during the time-frame of the index. This way we obtain a distribution forecast of the prices in every 5 min window during the last three hours before the delivery. An example of this approach can be seen in Figure (ref). We motivate our research with results on the weak-form efficiency of the market concluded by narajewski2019econometric and a possible application of the methodology to trading of the electricity and optimal redispatch management.
In purpose of modelling and forecasting of the trajectories, we utilize the generalized additive models for location scale and shape (GAMLSS) rigby2005generalized which extends the generalized additive models (GAM) hastie1990generalized. This methodology found applications to the electricity load pierrot2011short, gaillard2016additive and day-ahead price serinaldi2011distributional, gianfreda2018stochastic, abramova2020forecasting forecasting, but never to the intraday electricity markets. The model for price difference $\Delta P$ is fitted to the Student's t-distribution and mixed with the Dirac distribution, i.e. $\Delta P \sim (1-\alpha)\delta_0 + \alpha \text{t}$. $\alpha$ is assumed to be a Bernoulli variable with probability $\pi$ and is modelled using the logistic regression. We estimate it with the lasso method Tibshirani1996. A broader description of the modelling exercise can be found in Section (ref).
The forecasting part utilizes a rolling window study. This is a very common study type in the EPF and is widely utilized by researchers uniejewski2019understanding, narajewski2019econometric. We analyse both in-sample characteristics and evaluate the out-of-sample forecasting performance.
The major contributions of this paper are as follows:
The remainder of this paper has the following structure. In the next section, we describe the market. The third section consists of the data description and descriptive statistics. Then, a broader explanation of the estimation methods is presented, followed by the description of the considered models and benchmarks. In the fifth section, the forecasting study and evaluation measures are introduced and discussed in detail. In the sixth section, we present the results which consist of the in-sample analysis with relevant model interpretations and the out-of-sample evaluation. The final section concludes this paper. The methodology used in the paper is very innovative, especially in regard to the intraday electricity markets. We present it with an application to the German Intraday Continuous Market, but it can be easily used with any other intraday electricity continuous market.
The German Intraday Continuous Market allows to trade hourly, half-hourly and quarter-hourly products. We conduct the study using the most liquid part of the market -- the hourly one. This is in line with other EPF studies in intraday markets. Trading of hourly products in the German Intraday Continuous begins every day at 15:00 for the 24 products of the following day. It is possible to trade the electricity until 30 minutes (in the whole market) and up to 5 minutes (within respective control zones) before the delivery. In the meantime, between hour 22:00 and 60 minutes before the delivery the cross-border trading within XBID system is possible kath2019modeling. This system went live on 18th June 2018. A visualization of the trading timeline can be seen in Figure (ref). For more details on the German electricity market, we recommend the paper of Viehmann2017.
The most important price measure in the German intraday market is the volume-weighted average price of transactions in the last three hours of trading, called ID$_3$. The index takes into account only these transactions that happen until the gate closure 30 minutes before the delivery, so in fact it measures the last two and a half hours of trading before the gate closure. The relevance of ID$_3$ is an outcome of the behaviour of traders in the intraday market -- most of the transactions are held in this time period making it very liquid. This results in a high interest of practitioners and researchers in the ID$_3$-Price. For more details on the index visit the webpage of EPEX SPOT or see e.g. narajewski2019econometric.
To measure the prices during the trading period, we use the $_x$ID$_y$ defined by narajewski2019econometric. Let us recall the definition of $_x$ID$_y$. Let $b(d,s)$ be the start of the delivery of a product $s$ on day $d$. By $\mathbb{T}_{x,y}^{d,s} = \left[b(d,s) - x - y, b(d,s) - x \right)$, $x \ge 0$ and $y > 0$, we denote the time interval between $x+y$ and $x$ minutes before the delivery, and by $\mathcal{T}^{d,s}$ we denote a set of timestamps of transactions on the product. The $_x$ID$_y$ is defined by
where $V_k^{d,s}$ and $P_k^{d,s}$ are the volume and the price of $k$-th trade within the transaction set $\mathbb{T}_{x,y}^{d,s}\cap \mathcal{T}^{d,s}$ respectively. Let us note that the $_x$ID$_y$ is simply a volume-weighted average price of transactions in the time interval of length $y$ hours and ending $x$ hours before the delivery.
In the case of $\mathbb{T}_{x,y}^{d,s}\cap \mathcal{T}^{d,s} = \emptyset$ we use the value of $_{x+y}$ID$_y$, that is to say the previous observed volume-weighted average price measured on the time period of the same length.\footnote{In narajewski2019econometric this value is set to the price of the last transaction. This adjustment is caused by the fact that in this paper we work with 5-minute time intervals, leading to a significant number of the events of no trade in the time interval. This would often result in an artificial change of the price, compared to the previously observed $_x$ID$_y$.} In the case of no trades appearing since the start of trading, the price is set to the price of the corresponding Day-Ahead Auction.
The data used in purpose of this study consists of all transactions on hourly products in the German Intraday Continuous Market between 16th July 2015 and 1st October 2019. A more general descriptive statistics were presented by narajewski2019econometric. As mentioned in the previous section, the XBID system started to function on 18th June 2018. This means that XBID trades were possible only on around 30% of the days in the data. In the forecasting study, we use $D = 365$ days of the data as in-sample, and therefore the analysis in this section is based only on the initial in-sample, i.e. the data between 16th July 2015 and 14th July 2016. The start of the data is set to the first day of lead change in Germany from 45 min to 30 min in order to avoid this structural break. In this paper, we aggregate the transactions using the $_x\text{ID}_y^{d,s}$ with $y = 5$ min, and this way we obtain dense time series data. As said before, we are particularly interested in the evolution of prices during the last 2.5 hours of trading before the gate closure, so we use $x \in \mathcal{J} = \{180, 175, \dots, 35, 30 \}$, where $x$ is denoted in minutes. This way we observe $T = 31$ price points a day, what results in $T D = 31 \times 365 = 11315$ in-sample observations and $T$-dimensional simulated trajectories. Subsequently, we use a very specific setting, but it can be applied to any other continuous intraday market with other input variables.
As the market shows strong indications of weak-form efficiency, we focus on modelling of the price differences $\Delta P_t^{d,s} = {}_{(T-t) y + 30}\text{ID}_y^{d,s} -{}_{(T-(t-1)) y + 30}\text{ID}_y^{d,s} $ instead of pure prices $P_t^{d,s} = {}_{(T-t) y + 30}\text{ID}_y^{d,s}$. We also introduce the $P_t^{d,s}$ notation for simplicity. Due to the usage of price differences and to the fact that the data is aggregated using 5-minutes grid, we observe a high frequency of observations with no trade, and thus price differences equal to 0. This is depicted in Figure (ref). One can see that lack of transactions happens more often to the night and morning hours. In Figure (ref), we zoom in the tails of the histograms from Figure (ref). We also plot there densities of 4 distributions fitted to the data: the normal distribution $\mathcal{N}(0, \widehat{\sigma})$ and the t-distribution $\text{t}(0, \widehat{\sigma}, \nu)$ with fixed $\nu \in \{2.5,3,4\}$ and estimated $\widehat{\sigma}$ using maximum likelihood estimation ignoring the no-trade observations. Based on Figure (ref) it is clear that the price differences $\Delta P_t^{d,s}$ are heavy-tailed. One can see that even the t-distribution with $\nu = 4$ seems to be not heavy-tailed enough for the data. This indicates that the tail-index of the price differences may be lower than 4 which would mean that the fourth moment of the $\Delta P_t^{d,s}$ might not exist what is a strong indication for heavy tails.
Figure (ref) shows the frequency of the no-trade event over time to delivery. We see that the overall behaviour is very similar across all products -- the closer to the delivery, the less observations without transactions. What is different among the products is the level of the frequency. It is clear that the frequency decreases as the product time increases and the reason for it may be the time distance from the Day-Ahead and Intraday Auctions. It is intuitive that since these auctions the uncertainty could be smaller for the first products and higher for the last ones, but the smallest values of frequency are achieved not for the evening, but for the day-peak hours. This can be explained by higher activity in the market due to higher expected demand.
Figure (ref) shows the in-sample standard deviation of price differences $\Delta P_t^{d,s}$ over time to delivery. The dashed lines depict the standard deviation of the whole samples, independent of time. If the price processes would be similar to random walk, the sample standard deviation over time should be oscillating around these dashed lines. The behaviour in Figure (ref) is clearly different, with a spike in the last 30 minutes before gate closure. This suggests that the variance should be a subject of modelling. Figure (ref) presents the partial autocorrelation function of the absolute price differences $\left|\Delta P_t^{d,s}\right|$ to explore potential conditional heteroscedasticity in the heavy-tailed data. Figure (ref) shows that the most significant are the first three lags. Also, lags up to 6 may contain some information. Surprisingly, lags around 31 seem to be significant too, but this is most likely some daily dependence.
We assume the price differences $\Delta P_t^{d,s}$ to follow a 4-parametric distribution -- a mixture of the Dirac $\delta_0$ distribution and the 3-parametric t-distribution, sometimes referred as zero-inflated t-distribution:
where $\alpha_t^{d,s} = \mathds{1}(V_t^{d,s} \neq 0)$ is a Bernoulli variable of the event that there is non-zero volume of energy traded on product $s$ on day $d$ at time $t$ with probability $\pi_t^{d,s}$ and $F_t^{d,s}$ is the 3-parametric t-distribution $t(\mu_t^{d,s}, \sigma_t^{d,s}, \nu_t^{d,s})$ where $\mu_t^{d,s} \in \mathbb{R}$ is the mean, $\sigma_t^{d,s} > 0$ the standard deviation and $\nu_t^{d,s} > 2$ the degrees of freedom. The t-distribution is estimated with GAMLSS framework rigby2005generalized.
The GAMLSS is an expansion of the GAM hastie1990generalized and it allows to model not only the expected value of a response variable, but also potentially the higher moments, represented by scale and shape parameters. Namely, let $Y$ be a random variable with a density function $f(y|\Theta)$, where $\Theta$ is a set of up to four distribution parameters. Then each $\theta_i \in \Theta$ may be modelled by
where $g_i$ is some link function, $J_i$ is a number of explanatory variables and $h_{ji}$ is a smooth function of explanatory variable $x_{ji}$. Note that function $h_{ji}$ does not have to be a parametric function. In our exercise, we use the following link functions
$g_1$ is a standard link function for the expected value. $g_2$ is a link function that we call "logident" and we introduce it in order to avoid exponential inverse function for high values of estimates. The third link function is simply a natural logarithm shifted to 2 for preserving the condition that $\nu > 2$. The three link functions are plotted in Figure (ref). The models for $F_t^{d,s}$ are estimated using the \verb+gamlss+ package in R stasinopoulos2007generalized.
Due to the novelty of the exercise, we cannot use any literature benchmarks, as well as any standard approaches to the modelling of volatility, e.g. GARCH. Even though the data looks like time series, the biggest problem lies in the gap between days. We model each product separately, and for each product we have 31 observations every day. In the corresponding time series, the observations on day $d$ appear in 5-minute breaks, while the time difference between the last observation on day $d$ and the first on day $d+1$ is around 21 hours.
Furthermore, there is no direct link between the prices on day $d$ and day $d+1$ as they are for different delivery periods with potentially different fundamental market situations. Thus, the usage of GARCH-type components to address conditional heteroscedasticity is not straight-forward. Instead, as simple benchmarks we use models that assume the distribution of $\Delta \mathbf{P}^{d,s} = (\Delta P_1^{d,s}, \Delta P_2^{d,s}, \dots, \Delta P_T^{d,s})$ to be multivariate, random walk models, and a model that uses in-sample price differences to create an ensemble forecast. Also, as advanced benchmarks linear quantile regression with copula models are considered.
In the following subsection, the more complicated models are considered. We model explicitly the probability of non-zero number of transactions, the mean, and the variance of fitted distribution. We present the models from the least to the most complex and show the results similarly. This allows us to observe the gain caused by every new part of the model.
We introduce a dependency structure between the first three parameters of the $G_t^{d,s}$ distribution, i.e. $\pi_t^{d,s}$, $\mu_t^{d,s}$ and $\sigma_t^{d,s}$, and the data. For the fourth parameter, the degrees of freedom $\nu_t^{d,s}$, we assume the constancy. The $G_t^{d,s}$ distribution is estimated in a 2-step approach. First, the $\pi_t^{d,s}$ parameter is estimated, and then the $F_t^{d,s}$ distribution is fitted to the in-sample price differences $\Delta P_t^{d,s}$ for which the value of $\alpha_t^{d,s}$ is 1.
In the first step, we build a logistic model for $\pi_t^{d,s}$
The model explains the logit function with 4 main components: price difference impact, time dummies, fundamental regressors and regression on $\alpha_t^{d,s}$. Price difference impact consists of 3 most recent price differences, 6 most recent absolute price differences and a sum of absolute prices differences lagged by 7 to 12. This component addresses the overall impact of price volatility on $\pi_t^{d,s}$. We expect to observe more trades when the prices are more volatile. Time dummies consist of three weekday dummies and time to maturity dummies. The weekday dummies for Monday, Saturday and Sunday are chosen literature-based. A number of studies misiorek2006point, uniejewski2016automated, ziel2018day have proven that usage of these dummies in EPF substantially improves the forecasting performance. These three dummies indicate the end of the week with Monday being a transition day. The use of time to maturity dummies is clear when we take a look again at Figure (ref). It is expected that $\pi_t^{d,s}$ rises as we approach the gate closure. Fundamental regressors consist of day-ahead forecasts of total load, solar generation, wind onshore generation and wind offshore generation. It is expected that higher load and share of renewables should rise the uncertainty in the market, and encourage market participants to trade more. The last, but not the least is the regression on $\alpha_t^{d,s}$. We do not use the regression directly, but instead we use the average of last $j$ observed values of $\alpha_t^{d,s}$ which we denote by $\bar{\alpha}_{t-j}^{d,s}$. We expect these values to have a significant impact on the prediction of $\pi_t^{d,s}$. Intuitively, the higher these averages, the higher the value of $\pi_t^{d,s}$.
Model (ref) consists of 61 regressors in total. To avoid overfitting problems, we estimate the model using the least absolute shrinkage and selection operator (lasso) of Tibshirani1996. Let us recall that if we possess a logistic model $ \log\left(\frac{\pi}{1-\pi}\right) = \mathbf{X}' \boldsymbol{\beta}$ for the Bernoulli variable $\alpha$ with $P(\alpha = 1) = \pi$, then the lasso estimator $\widehat{\boldsymbol{\beta}}^{\text{lasso}}$ is given by
where $l$ is the corresponding log-likelihood
$\widetilde{\mathbf{X}}$ is a standardization of $\mathbf{X}$ and $\lambda$ is a tunable shrinkage parameter. This method found already many successful applications to the EPF and intraday markets ziel2016forecasting,uniejewski2019understanding, narajewski2019econometric. In this exercise, we utilize the \verb+glmnet+ package in R by friedman2010regularization. The estimation is conducted using a BIC-tuned $\lambda$ value chosen from an exponential grid of 100 values.
Let us now take a look at the $F_t^{d,s}$ distribution in equation (ref). We consider four versions of it. In the first one, we assume that $F_t^{d,s}$ follows $ t(0, \sigma_t^{d,s}, \nu_t^{d,s})$ with constant $\sigma_t^{d,s}$ and $\nu_t^{d,s}$. We denote it simply by Mix.RW.t. The $F_t^{d,s}$ distribution is fitted to the in-sample price differences with non-zero transaction number using the GAMLSS. With this model we can observe the gain of using a complex model for the $\pi_t^{d,s}$ parameter. Figure (ref) shows fitted densities to the histograms presented in Figure (ref). They were obtained with model Mix.RW.t.
The second model utilizes $F_t^{d,s}$ with modelled $\mu_t^{d,s}$ and constant $\sigma_t^{d,s}$ and $\nu_t^{d,s}$, and we denote it by Mix.t.mu. This model helps us understand the outcome of modelling of the expected value of $\Delta P_t$. However, a preliminary analysis has shown that most of the regressors used in model (ref) were not significant for modelling of $\mu_t^{d,s}$. The only significant were the three most recent price differences. Therefore, we model the expected value with
The next model uses $F_t^{d,s}$ with $\mu_t^{d,s} \equiv 0 $, modelled $\sigma_t^{d,s}$ and constant $\nu_t^{d,s}$. We denote it by Mix.t.sigma. The formula for the standard deviation is as follows
where $h_1$ and $h_2$ are smooth non-linear P-spline functions. The P-splines simply combine equally-spaced B-splines and discrete penalties. More information on P-splines can be found in eilers2015twenty. Let us note that the model described by equation (ref) uses much more regressors than in equation (ref). The explanation of the choice of the variables is very similar to the one of the model described by equation (ref). We explain the standard deviation of price differences with: lagged absolute price differences, weekday dummies, fundamental regressors, lagged values of $\alpha_t^{d,s}$ and non-linearities in most recent price and time to maturity variables. We expect that the absolute price changes are a suitable explanatory variable for the standard deviation as motivated through Figure (ref). The fundamental regressors are supposed to have a positive linear correlation with the $\sigma_t^{d,s}$. For the Saturday and Sunday dummies we might expect a negative impact due to lower trading activity on weekends, but also a positive impact due to the fact that higher bid-ask spreads are plausible. The lagged values of $\alpha_t^{d,s}$ indicate if the market participants traded lately, and thus we believe that it could identify higher price difference's variance. The last two regressors are expected to have a non-linear impact on the formation of $\sigma_t^{d,s}$, and therefore they are estimated using P-splines. Figure (ref) provides already an evidence that the standard deviation varies over time to maturity. Moreover, we suspect that extreme values of most recent price $P_{t-1}^{d,s}$ result in a higher variance due to a relatively inelastic supply curve in extreme price areas.
The last and at the same point the most complicated model uses $F_t^{d,s}$ with $\mu_t^{d,s}$ and $\sigma_t^{d,s}$ modelled and constant $\nu_t^{d,s}$. We denote it by Mix.t.mu.sigma. The $\mu_t^{d,s}$ is modelled using the formula from equation (ref) and the $\sigma_t^{d,s}$ using the formula from equation (ref). Let us mention that we could make the mixture model even more complex by modelling the degrees of freedom parameter $\nu_t^{d,s}$. However, a preliminary analysis has shown that it does not yield any significant improvement while increasing heavily the computational cost. Thus, in the forecasting study we analyse the performance of 8 models described in this section.
The first benchmark model uses one of $D = 365$ historical trajectories to model the price difference vector $\Delta \mathbf{P}^{d,s} = (\Delta P_1^{d,s}, \Delta P_2^{d,s}, \dots, \Delta P_T^{d,s})$. We denote it by Naive and its formula is given by
where $d' \sim \mathcal{U}(\{d-1,\ldots, d-D\})$ is a uniform random variable indicating the day used to model the price difference. Let us note that a fixed $d'$ index is used to model the whole price trajectory, i.e. for every $t \in \{1,2,\dots, T\}.$ This model assumes that the future trajectories can be forecasted using simply the past ones.
The second and the third benchmark models assume that the price difference vector $\Delta \mathbf{P}^{d,s}$ follows a multivariate normal and t-distributions, respectively. They are denoted by MV.N and MV.t and are given by
where $\boldsymbol{\varepsilon}^{d,s} \sim \mathcal{N}\left(\mathbf{0}, \boldsymbol{\Sigma}^{d,s}\right)$ in the case of $\mathbf{MV.N}$ and $\boldsymbol{\varepsilon}^{d,s} \sim t \left(\mathbf{0}, \boldsymbol{\Sigma}^{d,s}, {\nu}^{d,s} \right)$ in the case of MV.t. Let us note that the covariance matrix ${\boldsymbol{\Sigma}}^{d,s}$ and degrees of freedom $\nu^{d,s}$ are estimated by fitting the respective distributions to the in-sample observations. Moreover, the degrees of freedom $\nu^{d,s}$ is assumed to be constant for all $t \in \{1,2,\dots, T\}.$
The next benchmark model is the random walk version of the mixture model described by equation (ref), and we denote it by RW.t.mix.D. The formula is as follows
where $\varepsilon_t^{d,s} \sim G_t^{d,s}$ with $\widehat{\pi}_t^{d,s} = \frac{1}{D T} \sum_{i = d - D}^{d-1} \sum_{t=1}^{T} \mathds{1}(V_t^{i,s} \neq 0)$, ${\mu}_t^{d,s} \equiv 0$ and constant ${\sigma}_t^{d,s}$ and ${\nu}_t^{d,s}$. These values are estimated based on the in-sample data.
The fifth benchmark model is a modification of the RW.t.mix.D. We denote it by RW.t, and we simply set ${\pi}_t^{d,s} \equiv 1$ which means that we do not incorporate the mixing part and assume that the price differences follow the t-distribution. The last and the simplest of the random walk models assumes the price differences to follow a Gaussian distribution $\mathcal{N}(0, ({\sigma}_t^{d,s})^2)$ and is denoted by RW.N. In terms of the $G_t^{d,s}$ distribution, we simply modify the RW.t model by taking ${\nu}_t^{d,s} \to \infty$.
Later, we consider the random walk models from the simplest RW.N to the most complex RW.t.mix.D. This allows us to observe the gain of introducing more complex structure of the distribution. Let us note that model RW.N assumes exponentially decaying tails of the price differences $\Delta P_t^{d,s}$. Comparing it to model RW.t we measure the gain of assuming heavier, polynomially decaying tails. Based on the number of outlier observations in the German intraday market and on Figure (ref), we expect it to perform better than the Gaussian random walk. Then, considering the RW.t.mix.D helps us to understand the gain of the introduction of the mixture.
As mentioned, we are unable to use any literature-based benchmarks as this is the first paper on ensemble forecasting in intraday electricity markets. However, it is possible to implement scenario generating methods that are utilized in other research areas. Thus, as advanced benchmark we utilize a smoothed linear quantile regression model with two copulas: Gaussian and independence, and we denote them by LQR.Gauss and LQR.ind, respectively. A very similar approach was applied recently in the purpose of generating density forecasts of significant sea wave height and peak wave period gilbert2020probabilistic.
First, we build the linear quantile regression (LQR) model using the same set of regressors as for the mixture models. The formula is as follows
for $\tau \in \{0.01, 0.02, \dots, 0.99 \}$ and $t = 1, 2 \dots, T$. That is to say, we build separate models for each quantile $\tau$ and each time point $t$. Let us note that due to the design of the model, we can use only the regressor values available at the time of forecasting $t = 0$ (i.e. 3 h 5 min before the delivery). This results in the fact that here we model all $T$ time points using the same data, what is contrary to the mixture models where we can use autoregressive variables due to the recursive character of the models. We estimate the LQR models using the \verb|quantreg| package in R quantreg.
In the next step, a spline interpolation is applied over all fitted $Q^{\tau}\left(\Delta P_t^{d,s}\right)$ for $\tau \in \{0.01, 0.02, \dots, 0.99 \}$ and for every $t = 1, 2 \dots, T$. In order to preserve the monotonicity of the estimated cumulative distribution function (CDF) we compute a monotonic cubic spline using Hyman filtering hyman1983accurate. This way we obtain a smooth and monotonic $T$-dimensional CDF function
where $\Psi_t^{d,s}\left( \min_{\{d-1, \dots, d-D \}}\left( \Delta P_t^{d,s} \right) \right) = 0$ and $\Psi_t^{d,s}\left( \max_{\{d-1, \dots, d-D \}}\left( \Delta P_t^{d,s} \right) \right) = 1$. To assess the dependency structure of the price differences $\Delta P_t^{d,s}$ over $t$ we use two copulas: Gaussian and independence. The Gaussian copula for a given correlation matrix $\mathbf{R}$ can be written as
where $\Phi^{-1}$ is the inverse CDF of a standard normal distribution and $\boldsymbol{\Phi}_\mathbf{R}$ is a joint CDF of a multivariate normal distribution with mean vector zero and covariance matrix equal to the correlation matrix $\mathbf{R}$. We estimate the correlation matrix $\mathbf{R}$ simply by calculating the in-sample correlation matrix.
We use a rolling window forecasting study approach with $D = 365$ days in-sample size and $N = 1173$ days out-of-sample size. The in-sample data consists of $D T$ data points where $T=31$ in this study. We model each of the $S = 24$ hourly products separately and our forecasting time is 185 minutes before the delivery of product $s$ on day $D+1$. That is to say, we can utilize all the information from the in-sample data and from the day $D+1$ until 185 minutes before the delivery. At this time we forecast $M = 1000$ times the first price difference $\Delta P_1^{d,s} = {}_{180}\text{ID}_5^{d,s} - {}_{185}\text{ID}_5^{d,s}$. Based on these forecasts and explanatory data we simulate $M$ second price differences $\Delta P_2^{d,s}$ and we continue this recursive process until we reach the gate closure. Figure (ref) provides an outline of the exercise. This gives us $M$ simulated trajectories, each consisting of $T=31$ points. After that, we move the window forward by one day and repeat the exercise until the end of out-of-sample data. However, in the case of benchmark models we do not use the recursive algorithm as there is no recursion in their formulas.
Before we discuss the evaluation design in detail, we recall that we are mainly interested in the evaluation of the forecasted $T$-dimensional distribution of the price vector $\mathbf{P}^{d,s} = \left(P_{1}^{d,s},\ldots,P_{T}^{d,s}\right)$ which is represented by the predicted ensemble. Indeed, the multivariate cumulative distribution function of the ensemble coincides with the underlying cumulative distribution function if the ensemble sample size $M$ goes to infinity. Thus, for sufficiently large $M$ the evaluation of the scenario set can be regarded as the evaluation of probabilistic distributions.
From the theoretical point of view, strictly proper multivariate scoring rules are the first choice for evaluation, as they are able to identify the optimal forecast resp. the true distribution, see gneiting2007strictly, pinson2012evaluating. However, we want to remind us that even if forecast $A$ performs significantly better than forecast $B$ with respect to a strictly proper scoring rule, there is no guarantee that $A$ also performs better than $B$ in stochastic optimization problems (e.g. trading or storage optimization) where the forecasts are used as input. The optimal forecast would always yield optimal solutions in the stochastic optimization application. Thus, if $A$ is close to the optimal forecast with respect to a strictly proper multivariate scoring rule, the aforementioned risk that $B$ outperforms $A$ in the application is very limited if the stochastic optimization problem is continuous in the stochastic argument. Unfortunately, this only holds for strictly proper multivariate scoring rules. For proper scoring rules which identify only some characteristics of the full predictive distribution, this is certainly not true. The range of available strictly proper scoring rules is very limited, and reduces basically to the energy score for our ensemble forecasting problem, see e.g. lerch2020simulation. Therefore, we consider also proper scoring rules which might allow further insights as they focus on specific characteristics of the full distribution. To draw statistically significant conclusions on the outperformance of the forecasts of the considered models we utilize also the diebold1995comparing test.
As mentioned, the only available multivariate strictly proper scoring rule is the energy score (ES) gneiting2007strictly\footnote{Also the multivariate log-score is a known strictly proper scoring rule for multivariate distributions. However, it requires that the underlying multivariate distribution is continuous and has a density. Due to the non-trade events this is not satisfied for our forecasting problem.}. We compute the ES loss function in the following way
where
and
with $\mathbf{P}^{d,s} = \left( P_1^{d,s}, P_2^{d,s}, \dots, P_T^{d,s} \right)$ and $\widehat{\mathbf{P}}_j^{d,s} = \left(\widehat{P}_{1,j}^{d,s}, \widehat{P}_{2,j}^{d,s}, \dots, \widehat{P}_{T,j}^{d,s}\right)$. The $\text{ED}^{d,s}$ component measures the distance between the simulated trajectories and the observed prices. On the other hand, $\text{EI}^{d,s}$ measures the spread between the simulations. To calculate the overall energy score we use an average
We mentioned that the ES evaluates the full predictive distribution, which includes the path dependency in the generated scenarios. To illustrate the appropriateness of the energy score to evaluate correctly an ensemble forecasting study, we perform a short experiment in the results section. We take the best performing model and modify it with 3 different copulas which we refer as maximum dependency, minimum dependency and independence. For the maximum dependency copula we consider the $T$-dimensional co-monotonicity copula defined by $ M_{\max}(\boldsymbol{u}) = \min(u_1 , u_2 ,\ldots, u_T)$. The minimum dependency copula $M_{\min}$ is constructed using pairwise bivariate counter-monotonicity copulas defined by $W(u_1, u_2) = \max(u_1 + u_2 - 1, 0)$. So $M_{\min}$ is the copula that is associated with the $T$-dimensional uniform random variable $\boldsymbol{U} = (U_1,\ldots, U_T)$ that satisfies $(U_{t},U_{t+1}) \sim W$ for all $t=1,\ldots, T-1$. The $T$-dimensional independence copula is defined by $M_{\text{ind}}(\boldsymbol{u}) = \Pi_{t=1}^T u_t$. The 3 new models are evaluated using all considered measures and compared to the original one.
As pointed out by pinson2012evaluating, the ES does not evaluate the ability of the trajectories to mimic specific characteristics of the stochastic process. Therefore, we also focus our evaluation on specific characteristics of the underlying multivariate distribution. In this purpose, we consider additionally the subsequent proper scoring rules. We utilize the mean absolute error (MAE) and the root mean squared error (RMSE), pinball score (PB) to evaluate the median, mean and selected quantile trajectories, respectively, see e.g. gneiting2011making. For evaluation of the marginal density fit of our scenarios, we consider the continuous ranked probability score (CRPS) and additionally the empirical coverage of specific prediction intervals gneiting2007strictly. Moreover, we consider the variogram score (VS) and Dawid-Sebastiani score (DSS) which are regularly used to evaluate multivariate distributions lerch2020simulation. Note that both measures are only proper scoring rules and correct model identification fails in general, see e.g. ziel2019multivariate for empirical examples. Other scoring rules evaluating e.g. marginal distribution characteristics or specific events might also be added if it is relevant for the desired application.
The RMSE is the optimal least squares measure, i.e. it is the strictly proper scoring rule for mean evaluation while MAE is strictly proper for median evaluation. They are widely used both by researchers and practitioners. Their formulas are given by
and
where $\widehat{P}_{t,j}^{d,s}$ is the $j$-th simulation of $P_t^{d,s}$ and $\text{med}_{j = 1, \dots, M} \left(\widehat{P}_{t,j}^{d,s} \right)$ is the median of $M$ simulated $\widehat{P}_{t,j}^{d,s}$ prices.
We approximate the CRPS using the pinball loss
for a dense equidistant grid of probabilities $r$ between 0 and 1 of size $R$, see e.g. nowotarski2018recent. In this study, we consider $r = \{0.01, 0.02,\dots, 0.99\}$ of size $R = 99$. $\text{PB}_{t,\tau}^{d,s}$ is the pinball loss with respect to probability $\tau$. Its formula is given by
where $Q_{j = 1, \dots, M}^{\tau}\left(\widehat{P}_{t,j}^{d,s}\right)$ is the $\tau$-th quantile of $M$ simulated $\widehat{P}_{t,j}^{d,s}$ prices. To calculate the overall CRPS value we use a simple average
We can also use the pinball loss to compare the models' performance in particular quantiles. In this purpose the following formula is used
As mentioned, we use also the empirical coverage of prediction intervals, precisely the $(\tau/2,1-\tau/2)$-prediction interval. The $\tau\%$-coverage is calculated using the following formula
where $\tau \in \{0.5, 0.9, 0.99\}$.
The variogram score was introduced by scheuerer2015variogram in a probabilistic forecasting exercise for meteorological data. We compute it by
The Dawid-Sebastiani score evaluates the first and second moments gneiting2011making. It corresponds to the log-score of the multivariate normal distribution. We calculate it by
where $\widehat{\boldsymbol{\mu}}^{d,s}$ and $\widehat{\boldsymbol{\Sigma}}^{d,s}$ are the sample mean vector and the sample covariance matrix of the predicted price ensemble.
However, the aforementioned measures do not allow us to make conclusions regarding the statistical significance. To do so we utilize the diebold1995comparing test which tests forecasts of model $A$ against the ones of model $B$. The DM test is mostly used to evaluate point forecasts, but with correctly defined loss differential series it can be successfully applied in the evaluation of probability forecasts. We derive the series using ES and CRPS what was already applied by e.g. muniain2020probabilistic and lerch2020simulation. We utilize the multivariate version of the DM test as ziel2018day. The multivariate DM test results in one statistic for each model which is computed based on the $S$-dimensional vector of losses per day. Therefore, denote $L_A^d = (L_A^{d,1}, L_A^{d,2}, \dots, L_A^{d,S})'$ and $L_B^d = (L_B^{d,1}, L_B^{d,2}, \dots, L_B^{d,S})'$ the vectors of out-of-sample losses for day $d$ of models $A$ and $B$, respectively. By $L_Z^{d,s}$ we mean the $\text{ES}^{d,s}$ and $\text{CRPS}^{d,s}$ losses of model Z, formally we choose
The multivariate loss differential series
defines the difference of losses in $||\cdot||_1$ norm. For each model pair, we compute the p-value of two one-sided DM tests. The first one is with the null hypothesis $\mathcal{H}_0: \mathbb{E}(\Delta_{A,B}^d) \leq 0$, that is to say the outperformance of the forecasts of model $B$ by the ones of model $A$. The second test is the reverse null hypothesis $\mathcal{H}_0: \mathbb{E}(\Delta_{A,B}^d) \geq 0$. Let us note that these tests are complementary, and we assume that the loss differential series is covariance stationary.
We divided this section into two subsections: in the first one, we inspect the in-sample characteristics and in the second one, we present the out-of-sample simulation results.
We start our study with an analysis of the initial in-sample characteristics. Table (ref) shows the estimated coefficient values of model Mix.t.mu.sigma based on the initial in-sample data. The table reports the values for every hourly product, and it is split to 3 sub-tables, each regarding different parameter of the t-distribution. The first sub-table presents coefficients of the model described by equation (ref). Variable $\Delta P_{t-1}^{d,s}$ appears to be statistically significant for most of the hours. However, raising the lag decreases the significance. This behaviour goes in the direction of weak-form efficiency concluded by narajewski2019econometric.
The second sub-table shows coefficients of the model presented in equation (ref). Here, we see that all the variables using lagged absolute price differences are mostly significantly different from 0. Moreover, the coefficients of $|\Delta P_{t-1}^{d,s}|$ and $|\Delta P_{t-2}^{d,s}|$ are relatively high. Surprisingly, the day-ahead forecast of total load is mostly irrelevant. The day-ahead forecast of solar generation is significant mainly during the day-peak and in the evening. The day-ahead forecast of wind onshore generation appears to have a big positive impact on the volatility of price differences, in contrast to the wind offshore forecast. The behaviour of weekday dummies gives some light to our mixed expectancies -- they indicate a different behaviour of traders on Monday at night and on Saturday and Sunday during the day. On weekends the volatility is higher, likely due to higher bid-ask spreads on weekends. The lagged values of $\alpha_t^{d,s}$ have significant, negative impact on the volatility of price differences. It means that if there was trading at times $t-1$ and $t-2$, then the standard deviation would be lower. Let us note that the values of the coefficients are very similar among all hours except hours 14 to 16. For these hours the estimates of intercept and absolute price differences deviate heavily from the estimates of the remaining hours. A possible reason for this may be a few extreme outliers which were observed for these hours and for the others not. The table presents no values for the P-splines, because they are non-parametric functions. The last sub-table shows the estimate values for $g_3(\nu_t^{d,s})$ which we assumed to be constant. We show it anyway to gain an insight in the magnitude of the degrees of freedom. Let us recall that $g_3^{-1}(\nu) = \exp(\nu)+2$. Applying this to the estimate results in values of $\nu_t^{d,s}$ between around 3.7 and 6.6. Thus, the innovations are not extremely heavy tailed, and it is reasonable to apply asymptotic statistic for validation and interpretation.
Figure (ref) shows the initial in-sample P-splines $h_1(P_{t-1}^{d,s})$ and $h_2(t)$. We see that in case of both variables, the smoothing functions are non-linear. Extreme values of most recent price $P_{t-1}^{d,s}$ result in most cases in high rise of volatility. On the other hand, the values between 0 and 50 EUR/MWh have rather marginal impact on the variance of the price differences. An interesting effect can be seen in Figure (ref). We see that until 60 minutes before the delivery the impact on the volatility is on a similar, negative level among all products. Then, in the last 30 minutes of trading the volatility rises substantially above zero. This behaviour can be misinterpreted as a result of the closure of XBID as in Figure (ref). However, this plot is based on the initial in-sample data, i.e. the data between 16th July 2015 and 14th July 2016. Therefore, the effect of XBID could not be in the data as it was introduced on 18th June 2018. Figure (ref) is analogous to Figure (ref), but based on the first year of XBID, i.e. the data between 18th June 2018 and 17th June 2019. Comparing the two figures concludes that the introduction of XBID has an impact on the volatility of the price differences decreasing it even lower before the XBID closure and rising it even higher just after it. Interestingly, this is in contrary to the paper of kath2019modeling who concluded that there is no evidence for the influence of XBID on the price volatility.
Figure (ref) presents a price trajectory and a decomposition of fitted $g_2(\sigma_t^{d,s})$ of the product with delivery at 12:00 for the last 7 days of in-sample data. For the sake of readability we grouped the components of the model for standard deviation similarly as in equation (ref). Let us note that the absolute price differences and fundamental regressors have big, positive impact on the volatility of price differences. We also observe overall higher volatility on the weekend, i.e. the second and third day on the plot, than on the week. Note that in this specific example the impact of the non-linear price due to $h_1$ looks rather negligible. However, the price level in these seven trading sessions is always between 0 and 40 EUR/MWh where we expect minor impacts.
Now, we turn ourselves to the analysis of the simulated trajectories. Figure (ref) shows the first out-of-sample simulation exercise of prices of product with delivery at 12:00 on 15.07.2016. The trajectories are simulated from Mix.t.mu.sigma model and it can be easily compared to the simulations from Gaussian random walk presented in Figure (ref). It is clear that in this example the trajectories of the mixture model are less volatile than the random walk.
Table (ref) shows the values of utilized error measures. The Naive model performs very well overall. Moreover, it gives the best results in terms of 50% coverage. Very similar results to the naive model gives the MV.t which assumes the trajectories to follow a multivariate t-distribution. Having a look at the performance of MV.N we see that indeed the t-distribution is better in modelling of the trajectories. The LQR-based models are according to most of the considered measures worse than the Naive or MV.t. It is worth to emphasize the very bad ability to model the mean and median trajectories and on the other hand quite good 99% coverage and the values of VS and DSS. Let us note a very bad performance of the Gaussian random walk. Model RW.N is clearly the worst. Having a look at its coverage values, we conclude that its simulations are too volatile. The introduction of t-distribution to random walk yields already a big improvement. Another step in our modelling, the usage of simple mixture distribution of the Dirac distribution $\delta_0$ and the random walk with innovations from t-distribution do not improve the results. However, the next step, i.e. modelling of the probability $\pi_t^{d,s}$ with model (ref) improves the results, but still they are not better than the ones of model \textbf{RW.t}. Moreover, modelling of the expected value as in equation (ref) also worsens the performance substantially. All these models are clearly worse than the \textbf{Naive} considering almost every measure. The last change to the mixture model, i.e. modelling of the standard deviation according to the formula in equation (ref) lowers the errors significantly. Modelling of the expected value in addition to the standard deviation brings a little improvement. Model \textbf{Mix.t.mu.sigma} is marginally better than \textbf{Mix.t.sigma} and it turns out to be the best model in terms of ES, CRPS, VS, DSS, MAE and RMSE. A little disturbing are the values of the 50%- and 90%-coverage which are too high for the mixture models. This means that it is very likely that the results can be still improved. On the other hand, they capture better the behaviour in the tails than the \textbf{Naive} model.
The values of the error measures in Table (ref) may suggest that both ES and CRPS evaluate the same thing -- the marginal distribution. To emphasize that ES evaluates also the quality of the generated scenarios, we perform a small experiment on the model Mix.t.mu.sigma. In Table (ref), we compare the performance of the Mix.t.mu.sigma with its copies modified using 3 copulas: maximum dependency, minimum dependency and independent. This results in the same marginal distributions, mean and median trajectories, and coverage values, but in completely different ensembles. This is depicted by the values of the measures -- the CRPS remains unchanged while the ES, VS and DSS changed drastically. Let us note the enormous aggravation of the DSS which is particularly sensitive to changes of the dependency structure.
Figure (ref) shows the models' performance over all products in terms of energy score. A very interesting is the case of model RW.N. Usually it is not that much worse than the other random walks, but for hours 14-16 the error explodes. A look into the data explains the situation clearly -- there were a few in-sample observations of extreme price differences. The normal distribution assumes exponentially decaying tails, and thus the model overestimated the variance. This indicates clearly that the t-distribution is better in this purpose as it was unaffected by these events. Furthermore, we observe that models with modelled $\sigma_t^{d,s}$ are uniformly better than the others.
Figure (ref) presents the models' performance over time to delivery. The values rise as the time goes, but it is rather not surprising. An interesting behaviour can be observed from 150 to 100 minutes before the delivery. In this time range the errors of the random walk models and the mixture models that assume constant standard deviation rise significantly in comparison to the other models. It is also the time of decreasing volatility in Figure (ref).
Pinball Score values over quantiles $\tau$ are depicted in Figure (ref). Let us note that the gain from the forecasting of central quantiles is marginal, and it is inline with other studies regarding the ID$_3$-Price in the German intraday market narajewski2019econometric,janke2019forecasting. On the other hand, models Mix.t.sigma and Mix.t.mu.sigma gain a lot from the forecasting of quantiles outside the centre, performing especially well in the tails. In relation to the naive benchmark, the error is around 30% lower in the lower tail and around 25% lower in the upper tail. Let us also note that the LQR-based models give quite good results in the tails, but lose a lot in the centre, compared to the naive or to the models with non-constant $\sigma_t^{d,s}$.
Figure (ref) shows the results of the DM test using two types of losses: the $\text{ES}^{d,s}$ and the $\text{CRPS}^{d,s}$. Before applying the test we conducted on the multivariate loss differential series $\Delta_{A,B}^d$ three tests that evaluate the null hypothesis that a unit root is present in the series against the alternative that the data is stationary or trend-stationary. We used the Dickey-Fuller test dickey1979distribution, the Augmented Dickey-Fuller test said1984testing, and the Phillips-Perron test phillips1988testing. In 99% of cases the obtained p-values were smaller than 0.01, so we reject the null hypothesis. This indicates in our case that the loss differential series is covariance stationary. Only for the differences with RW.N the Dickey-Fuller test reported no significance for rejecting the null hypothesis. This may be caused by the bad capturing of the marginal distribution of the RW.N. The results of the DM test show that the difference between the forecasts of models Mix.t.mu.sigma and Mix.t.sigma is insignificant. Moreover, these models give better forecasts than all the other considered models. It is worth to emphasize a very good performance of the Naive model, but it is not surprising after taking a look at Table (ref).
We conducted an ensemble forecasting study in the German Intraday Continuous Market which is novel in two ways. The first way, this study is the first one that raises the issue of price trajectory simulation and ensemble forecasting in continuous intraday electricity markets. The second way, the study uses a very clever mixture of distributions that is fitted to the data. The results are very satisfying and showing that it is possible to successfully model the volatility in the German Intraday Continuous Market. The study was carried out using the data from the German market, but the generality of this method and the organization of the European electricity markets ensure a possible application to other markets, especially the markets participate in XBID which covers exchanges like EPEX, Nordpool and OMIE.
Obviously, the proposed method can be developed further. One of possible directions is using other external processes like the traded volume or price of nearby hours as regressors. Although, this is a non-trivial task and could easily lead to the accumulation of errors. Another possibility is utilization of other probability distribution. The not perfect coverage of the best performing model indicates that there is still some space for improvement. This issue could be addressed with some post-processing method as well.
This research article was partially supported by the German Research Foundation (DFG, Germany) and the National Science Center (NCN, Poland) through BEETHOVEN grant no. 2016/23/G/HS4/01005.