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.
306,961 characters · 20 sections · 102 citation commands
Simulation-based Forecasting for Intraday Power Markets: Modelling Fundamental Drivers for Location, Shape and Scale of the Price Distribution
\setcounter{tocdepth}{2}
\deffootnote[1.5em]{0em}{1.5em}{\makebox[1.5em][l]{\thefootnotemark )\ }} Keywords: electricity price forecasting, volatility forecasting, intraday energy market, auction curves, gamlss
Intraday power markets gained significant importance throughout the past few years, most visible in sharply increased trading volumes. In the European power market structure, they provide traders, asset owners and marketers of intermittent \gls{RES} the opportunity to balance forecast errors arising after the day-ahead auction until five minutes before the beginning of the delivery period hirth2019. Increasingly, this balancing action is taken over by algorithmic trading strategies, for which reliable short-term price and volatility forecasts are necessary. \textcolor{Black}{At the same time, our results shed light on the influence of idiosyncratic features of intraday markets such as \gls{SIDC} on the price process and are thus valuable for policy makers concerned with short-term markets.}
The recent and still sparse literature on probabilistic forecasting in intraday markets janke2019, ziel2020b, uniejewski2019 and the markets' driving fundamentals has, so far, focussed on modelling the impact of renewable forecast ziel2017, kath2019, pape2016, gurtler2018, balardy2018 and forecast errors ziel2017, kulakov2020, kuppelwieser2021. A different strand of literature emerged around modelling of the merit-order effect for price changes and price elasticity kiesel2017, kiesel2020a, kiesel2020b, kulakov2019, balardy2018. To the best knowledge of the authors, only ziel2020a and baule2021 focus on modelling the volatility in intraday markets. This paper aims to generalize the above research by investigating the fundamental drivers of the location, scale and shape parameters of the intraday price return distribution. We use the \gls{GAMLSS} framework to model the distribution moments in a parametric and explainable fashion. Our contribution is thus two-fold: We are able to significantly improve forecasting performance compared to benchmarks and qualitatively analyse the impact of fundamental drivers for the distribution moments. While this paper focuses on the German intraday market, our methodology is transferable to any continuous intraday market such as France, Great Britain, Spain or Turkey.
This paper builds on the work of ziel2020a to develop a simulation-based probabilistic forecasting model for the path of the five-minute-volume-weighted price \gls{symb:P} between 185 and 30 minutes before delivery. Instead of directly modelling the price $P$ as it is common in day-ahead forecasting and is done in other forecasting studies on the intraday market uniejewski2019, janke2019, the first differences $\Delta P$ will be modelled. The path of $P$ is thus the cumulative sum of the initial price $P_0$ and all price differences in the forecasting period. Following the suggestions of ziel2020a, we assume the first differences \gls{symb:delta_P} follow a mixture distribution of the Dirac distribution $\delta_0$ \textcolor{Black}{with an atom at 0} and a continuous distribution \gls{symb:F}. In ziel2020a, the latter is assumed to be $t$-distributed without any detailed justification, except the observation that $\Delta P$ tends to be heavy tailed. This manuscript extends the approach in four dimensions:
The extended models are tested in a forecasting study and compared to the benchmark models. We evaluate the probabilistic forecasting performance by utilizing established probabilistic scoring rules and calibration measures. Statistical significance is evaluated by the widely used \gls{DM}-test ziel2020a, ziel2019b, weron2018, diebold2002. \textcolor{Black}{In our forecasting study, the GAMLSS-based model assuming Johnson's $S_U$ significantly outperforms all proposed benchmark models as well as the GAMLSS-based model assuming the popular skew-$t$ distribution. The GAMLSS-based model assuming the skew-$t$ distribution exhibits stark sensitivity towards outliers.} Qualitatively, our results indicate that price changes $\Delta P$ in the intraday market are influenced by the first lag, while other explanatory variables have little predictive power for the expected value of $\Delta P$. This result supports the notion of weak-form efficient markets already indicated by ziel2020b, kuppelwieser2021. We find evidence for a merit-order effect in the volatility and kurtosis of the distribution $\Delta P$. A steeper merit-order implies higher volatility and heavier tails. Additionally, the volatility rises with decreasing time to delivery and with the gate closure of XBID/SIDC, while kurtosis is more driven by trading-related variables such as lagged absolute price differences. We find that none of the included explanatory variables has predictive power for the skewness of the distribution.
The presented models and methodology are also of interest to practitioners in intraday markets. Path-based forecasts allow to price short-term asset optionality using Asian option valuation. Additionally, the explicit modelling of the volatility provides a starting point to introduce time-varying volatility to mathematical finance models for market making and position solving luckner2017, glas2020, aid2016, kath2020.
The remainder of this paper is structured as follows: Section (ref) gives a short introduction to the structure of the German short-term power markets. Section (ref) presents the data preparation of the intraday trade data, the forecast and outage data sets and the transformation of the day-ahead auction curves and related assumptions. Also, some exploratory data analysis is carried out in this Section. Section (ref) introduces the used models. The forecasting study design and scoring rules are discussed in Section (ref). Finally, Sections (ref) and (ref) present the results and conclude this paper.
\FloatBarrier
This section briefly introduces the relevant structure of the German power market. As we work with data from the day-ahead auction and the intraday market, the description focuses on these markets. Generally, denote the delivery day as \gls{symb:d} and the delivery hour as \gls{symb:s} for $s = 0, ..., S$ and $S=23$. Times are usually expressed in local time unless otherwise noted.\footnote{\gls{CET} respectively \gls{CEST} for Germany.} Electric power markets generally follow the structure of a forward market, where different delivery periods in the future can be traded almost up to the start of actual physical delivery. Figure (ref) shows the time line of the German short-term markets, for more details see viehmann2017. \textcolor{Black}{Let us generally note here that we place indices referring to the delivery periods $d,s$ as superscript, while placing indices relating to the time where the price is determined as subscript. The same holds for other variables as e.g. production forecasts.}
The spot market is organized as a pay-as-cleared auction. The order book closes on $d-1$ at 12:00 and first auction results are published around 12:42 on $d-1$. Official results shall be published at latest at 14:00 on $d-1$. From the bids submitted by the market participants, EPEX Spot calculates aggregated supply and demand curves for each delivery period. The intersection between supply and demand curves is the market clearing price $\ensuremath{P^{d,s}_\text{DA}}$. Additionally to normal bids, market participants can submit special bids such as block bids spanning more than one delivery period and linked bids, where execution is linked to neighbouring bids. The day-ahead spot price also serves as reference price for cascading financial futures. epex2020c publishes aggregate curves together with the official market results. \textcolor{Black}{The lower price level is set to -500 EUR/MWh, the upper level is set to 3000 EUR/MWh}.
\textcolor{Black}{The intraday market is structured as continuous pay-as-bid auction similar to financial markets. However, contrary to equity or currency markets, the individual trading sessions of the intraday electricity markets are not part of a larger process, as the intraday trading session ends with the physical delivery of power. Hence, trading sessions for the same delivery period on different delivery days might be driven by fundamentally different circumstances and need to be viewed separately.} Trading starts at 15:00, 15:30, 16:00 on $d-1$ for hourly, half-hourly and quarter-hourly products with delivery on day $d$. At 18:00 on $d-1$, cross-border trading within the \gls{SIDC} system, formerly known as \gls{XBID}, starts in Germany, Denmark, Netherlands, Norway and Poland.\footnote{For the sake of consistency, both the \gls{XBID} and \gls{SIDC} are referred to as \gls{SIDC} throughout this paper.} At 22:00 on $d-1$, the remaining countries of the core market area follow. Here, the intraday order books of all participating countries are shared and orders can be matched internationally as long as sufficient transmission capacity is available. For each product, the cross-border shared order books close one hour before delivery. \gls{SIDC} went live on June 18, 2018 nordpool2018. 30 minutes before delivery, the Germany-wide order book closes and trading resumes in local (control zone) products up to five minutes before delivery. \textcolor{Black}{Note that all open delivery periods are traded in parallel.} The market price limits are at $\pm9999$ EUR/MWh. The smallest possible price tick changed multiple times throughout the last few years and is currently set to 0.01 EUR/MWh. The smallest possible volume tick is 0.1 MW epex2020b, epex2018, viehmann2017.
\FloatBarrier
\textcolor{Black}{On the intraday market, trading happens continuously. Hence, the transactions are irregularly spaced and need to be aggregated. The following paragraphs and Figure (ref) give a brief overview of the aggregation. A detailed description can be found in Appendix (ref).} Trade data is obtained from epex2020a. The data consists of all hourly trades on the continuous intraday market between January 1st, 2016 and July 31st, 2020.
For each delivery period $d, s$ we aggregate all trades on an equidistant 5-minute grid by taking the volume-weighted average price within each bucket, denoted by $\ensuremath{P^{d,s}_{\text{ID},t}}$, where $t$ denotes the 5-minute interval (see panel 2 in Figure (ref)). We then take first differences $\ensuremath{\Delta P^{d,s}_{\text{ID},t}} = \ensuremath{P^{d,s}_{\text{ID},t}} - \ensuremath{P^{d,s}_{\text{ID},t-1}}$ (see panel 3). Lastly, we define a boolean variable $\ensuremath{\alpha^{d,s}_t}$, which takes the value 1 if there has been at least one trade within the 5-minute interval (see panel 4). As the trading sessions in the intraday market are of varying length for the different delivery periods and our simulation concerns the last 185 minutes of trading for each product, we define $t$ relative to the start of the physical delivery. $t=1$ denotes the first 5-minute interval in the simulation window, thus 185 to 180 minutes before the start of physical delivery and $t = 31 = T$ denotes the last 5-minute interval in the simulation window, 35 to 30 minutes before the start of physical delivery. Similar aggregation methods have been used by ziel2020a, ziel2020b and serafin2022.
\textcolor{Black}{Figure (ref) shows the relationship between the share of no-trade events, i.e. 5-minute intervals where $\alpha_t^{d,s} = 0$, relative to the time to delivery on the initial training set. With decreasing time to delivery, the probability of no-trade events decreases in a non-linear fashion. For periods close to 30 minutes to delivery, the share of no-trade events in the initial training data set is close to 0, while further away from delivery, there are more periods without trades. For products with delivery in the peak hours, there are less no-trade events at the beginning of the $\text{ID}_3$ period already.} Additionally, Table (ref) presents summary statistics for $\alpha_t$ and $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ for all 5-minute intervals with at least one trade, grouped by year. The share of 5-minute intervals where $\alpha_t = 1$, i.e. at least one trade happens happens, increases throughout the years. It is almost 1 from 2018 onwards, implying that there are barely any periods without trades. Accordingly, the number of observations for $\alpha_t$ and $\ensuremath{\Delta P^{d,s}_{\text{ID},t}} \mid \alpha_t = 1$ converge. \textcolor{Black}{We can thus identify two levels of time-varying behaviour of $\alpha^{d,s}_t$, first across the multiple years of the data set, but also second within each trading session. While we explicitly model the latter, the first will be coped with due to the set-up of a rolling window forecasting study.}
\textcolor{Black}{The mean and median values of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ are close to 0 across all years in the dataset. However the standard deviation is rather high and the extreme minima and maxima already hint at a leptokurtic distribution of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$.} The minima and maxima increase throughout the data set, \textcolor{Black}{while the 5% respectively 95% and the 10% respectively 90% quantiles are roughly constant.} Especially for 2020, the minima and maxima of more than 2000 respectively less than -2000 EUR/MWh are noteworthy. Driven by these larger outliers in 2020, the standard deviation of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ rises fourfold between 2019 and 2020, while staying roughly constant before. \textcolor{Black}{The more robust dispersion measures median absolute deviation (MAD) and the interquartile range (IQR) support this notion.}
Figure (ref) plots histograms for $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ in the initial training set in the hours $s \in \{4, 12, 20\}$ exemplary. These delivery hours represent the typical night, noon and afternoon peak load hours. The first plot focuses on the general shape of the distribution as well as the relation between intervals with and without trades. \textcolor{Black}{The center bar shows the relative weight of 5-minute intervals without trades (i.e. $\ensuremath{\alpha^{d,s}_t} = 0$) and 5-minute intervals with at least one trade (i.e. $\ensuremath{\alpha^{d,s}_t} = 1$), but small or no price changes. As visible already in Figure (ref), the share of 5-minute intervals with $\ensuremath{\alpha^{d,s}_t} = 0$ decreases strongly for delivery hours after 8.} In the second Figure, the tails of the distribution are shown together with fitted normal, student-$t$, and Johnson's $S_U$ distributions. \textcolor{Black}{Additionally, Figure (ref) (a) plots the pearson autocorrelation coefficients $r$ for $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ for each trading session for the first lag. Colour intensity corresponds to the coefficient size. For the first lag, slight positive autocorrelation is present for the morning hours, while some negative correlation is visible for noon and evening hours. Figure (ref) (b) shows the according $p$-values for the test statistic $r \cdot \sqrt{n-2} / \sqrt{1 - r^2}$ where $n$ is the number of 5-minute intervals in the trading session. We find that for around one-fifth of all trading sessions, the lag 1 autocorrelation coefficient is significant at the 10%-level and for only 15% of all trading sesions, the lag 1 autocorrelation is significant at the 5%-level. For lags 2 and 3, we find even less significant autocorrelation (see Figures (ref) and (ref) in Appendix (ref)).}
\FloatBarrier
The relationship between the realized variance and the time to delivery in the initial training data is shown in Figure (ref). The volatility increases slightly until 60 minutes before the start of the delivery and rises sharply between 60 and 30 minutes before the start the delivery. \textcolor{Black}{Here, we note three levels of time-varying behaviour: first, across the full data set, volatility is increasing. Second, within each day, the volatility moves with the peak/off-peak hours. Third, within each trading session, volatility increases towards the end of the trading session.}
\textcolor{Black}{To analyse the stationarity properties of the differenced and un-differenced price series, we apply the augmented Dickey-Fuller test to each simulation window individually and report aggregate results in Table (ref). For the majority of the trading windows, we find stationarity of the price differences and unit-root behaviour in the prices. We note though, that due to the heteroskedasticity present in the individual trading windows, the underlying assumptions of the ADF-test might be violated. Together with the results of lohndorf2022value, who aggregate trades in the intraday market on a 1-hour grid and report similar results for the ADF-test at the 10%-level, we conclude that the price changes in the intraday market are stationary.}
\textcolor{Black}{Intra-daily updated renewable production forecasts used in this paper are provided by Statkraft Markets and generated by statkraft2020a. Day-ahead demand / system load forecasts are obtained from entsoe2020. The forecasts have a 15-minute delivery period resolution and denote the expected produced power by all assets of the respective technology in Germany in MW. Forecasts are sampled to hourly frequency using a simple arithmetic average. A new update is available every hour. The first forecast version is issued several days before the delivery day, the latest version usually after the end of the delivery period due to ex-post updates. Let $\ensuremath{\hat{W}_{v}^{d, s}}, \ensuremath{\hat{S}_{v}^{d, s}}$ denote forecasts for wind and solar production for delivery period $d, s$ available at time $v$. Forecasts for demand are not updated as regularly, hence intraday-updates are not considered in this paper. We denote demand forecasts as $\ensuremath{\hat{D}_{\text{DA}}^{d, s}}$. Note that the issuance time $v$ of a new forecast does not necessarily correspond to the timing of trades on the continuous market or the 5-minute intervals used to aggregate these trades. For any forecasting study, it is important to keep in mind the information set at the point of forecasting. The start of the simulation is set to 185 minutes before the start of the delivery period. Hence, forecast versions and updates can only be considered if they are available earlier than 185 minutes before the start of the delivery period.\footnote{For example, for a product with delivery on September 1st, 12:00 to 13:00, all forecast versions available until September 1st, 8:55 can be used. For the product with delivery 13:00 to 14:00, all versions up to 9:55 can be used.} For each delivery period, two forecast versions deserve special attention: First, the latest forecast available before 12:00 on $d-1$, the deadline for submission of bids to the spot auction, is referred to as the day-ahead forecast $\ensuremath{\hat{W}_{\text{DA}}^{d, s}}$ and $\ensuremath{\hat{S}_{\text{DA}}^{d, s}}$. Secondly, the newest forecast available before the start of the simulation, i.e. at which $v \geq b(d,s) - 185$ holds, is denoted as the intraday forecast $\ensuremath{\hat{W}_{\text{ID}}^{d, s}}$ and $\ensuremath{\hat{S}_{\text{ID}}^{d, s}}$.}
\textcolor{Black}{ An initial analysis showed that individual forecast updates immediately before the start of the simulation carries little predictive power for the whole simulation period of three hours. Therefore, the forecast updates are aggregated. We consider two aggregated measures for forecast changes: first, the aggregated change between the production forecasts at the day-ahead stage and the production forecasts at the start of the simulation. Second, we employ the volatility of all forecast changes between the day-ahead stage and the production forecasts at the start of the simulation. Let us generally define the change between two forecast versions $v_1, v_2$ as $\ensuremath{\Delta\hat{W}_{v1, v2}^{d, s}} = \ensuremath{\hat{W}_{v_2}^{d, s}} - \ensuremath{\hat{W}_{v_1}^{d, s}}$ with $v_2$ being the newer forecast. }
\textcolor{Black}{
Analogously, $\ensuremath{\Delta\hat{W}_{\text{DA}, \text{ID}}^{d, s}}, \ensuremath{\Delta\hat{S}_{\text{DA}, \text{ID}}^{d, s, +}}, \ensuremath{\Delta\hat{S}_{\text{DA}, \text{ID}}^{d, s, -}}$ and $\ensuremath{\sigma^{d,s}_{\text{DA},\text{ID}}(\Delta \hat{S})}$ are defined for the solar production forecasts. Panels (a) - (c) in Figure (ref) show the day-ahead versions of wind, solar and demand forecasts. For wind and solar, the change between day-ahead and intraday versions and the standard deviation are plotted as well. }
Under the \gls{REMIT}, market participants are required to report non-availabilities of their assets and make this information available to all other market participants in order to avoid insider trading. In practice, this obligation is fulfilled by market participants by submitting non-availability messages to an inside information platform eu2011, lazarczyk2018, acer2020. We retrieve unavailability messages from the eex2020 market transparency platform for all non-availabilities regarding the delivery periods between January 1st, 2016 and July 31st, 2020. A non-availability message is defined by the date of publication, beginning and end of the non-availability, the type of non-availability, i.e. whether it has been planned or unplanned, the fuel type of the unavailable asset as well as the unavailable capacity in MW. The outage messages are aggregated to the total non-available generation capacity for the delivery period $d, s$ known at the time of the spot auction. Additionally, the outages are aggregated to the total non-available generation capacity known at the start of the simulation for a delivery period $d,s$. Sub-hourly outages are taken into account with the respective share of the full hour. The differences between the level of outages day-ahead and at the start of the simulation is calculated similar as the difference in the forecasts: $\ensuremath{\Delta O_{\text{DA}, \text{ID}}^{d, s}} = \ensuremath{O_{\text{ID}}^{d, s}} - \ensuremath{O_{\text{DA}}^{d, s}}$. Afterwards, the difference $\Delta O^{d,s}_{\text{DA},0}$ is split into planned and unplanned outages denoted as $\ensuremath{\Delta O_{\text{DA}, \text{ID}}^{d, s, \text{planned}}}$ and $\ensuremath{\Delta O_{\text{DA}, \text{ID}}^{d, s, \text{unplanned}}}$. Figure (ref) (d) plots the aggregated outage data; Table (ref) gives summary statistics.
\FloatBarrier
The impact of forecast errors on $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ depends on the steepness of the merit-order kiesel2020a, kiesel2020b, kulakov2020. Following this thought, the volatility of the intraday price should also be influenced by the slope of the merit-order. If the market is in a steep merit-order regime, even small volume changes might have a high price impact. Thus, the expected price impact of changes in (\gls{RES}) supply is stronger. Under uncertainty of future \gls{RES} forecast updates, the expected volatility should increase with the steepness of the merit-order. If the market price corresponds to a rather flat region of the merit-order, the price impact of changes in \gls{RES} supply should be smaller and hence the volatility of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ should be lower. This thought will be the main intuition for the addition of a merit-order slope to the model for the volatility of the distribution of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$.
There are different methods to model the merit-order used in practice and academia. Fundamental methods as developed by pape2016, gurtler2018, beran2019 are complex, data-intensive and rely heavily on assumptions. For this reason, kiesel2020a, kiesel2020b develop an econometric model based on he2013 by fitting the relationship between demand forecasts and day-ahead prices to an exponential function. This yields an analytically traceable function, whose slope can easily be calculated as the first derivative. This paper develops a further method to derive the slope of the merit-order by using the day-ahead auction curves as a proxy for the supply stack. This approach is based on balardy2018 and kulakov2020 and has three advantages compared to the approach of he2013: First, the auction curves combine all market and availability information available on $d-1$ and do not depend on a longer time frame for the estimation of the function coefficients. Second, by using the auction curves, there is no need to assume an explicit functional form for the merit-order. Lastly, the auction curves also represent negative prices, while the exponential function is only defined on the positive real line.
However, it is also important to discuss the drawbacks attached to modelling the intraday merit-order based on day-ahead information in general and attached to the auction curves especially. First, the available generation capacity can (and due to \gls{RES} will) change between $d-1$ and $d$, leading to shifts in the merit-order. Second, power plants might not be as flexible intraday as in a day-ahead planning horizon due to ramping behaviour, start-up costs or constraints due to grid service delivery. On the contrary, some power plants might be optimised predominantly intraday and not on the day-ahead auction if they are at-the-money \textcolor{Black}{pape2016}. Especially for the auction curves, two further problems arise: First, the demand and supply curve are both elastic curves, contrary to the common assumption of largely inflexible demand in energy markets. This problem is addressed by applying the transformation introduced by kulakov2020, kulakov2019 in the following paragraph. Thereby, all elasticity from the demand curve is shifted to the supply curve, which yields a perfectly inelastic demand and elastic supply curve. Second, the day-ahead auction curves as provided by EPEX Spot only contain standard bids. Thus, linked, block and other complex bids are excluded from the curves, which removes information about the available generation capacity. This problem cannot be addressed simply and needs to be kept in mind for the further interpretation of the results.
The intuition behind the transformation of the auction curves is outlined in detail in kulakov2019 and coulon2014, so here only a brief introduction is given. Figure (ref) shows that the demand curve at the day-ahead auction is elastic, which is at odds with the common assumption of few price elastic consumers of electricity, especially at short notice knaut2016, coulon2014. However, producers and consumers have the chance to sell/purchase their energy not only on the spot auction, but also in the OTC and derivative markets. In addition, there might be market participants that own both assets on the supply and demand side. Thus, arbitrage opportunities between the two markets arise that can be used by the trader. coulon2014 and kulakov2019 consider this effect by flipping the elasticity from the demand curve to the supply curve, hence obtaining a perfectly inelastic (vertical) demand curve and an elastic supply curve to incorporate those effects. \textcolor{Black}{The core idea here is that, at the day-ahead auction, placing a buy order for a volume $x$ for a price $y$ is the same placing a buy order for the volume $x$ at the maximum price and placing a sell order with volume $x$ for the price $y$ + the smallest tick. kulakov2019 elaborate in detail on the econometric framework, which is adopted in this paper and the implications for the different market participants.}
Figure (ref) shows the supply and demand curves from the spot auction for the delivery day June 1st, 2017 for hour $s=9$. The intersection of supply and demand yields the spot price $P^{d,s}_\text{Spot}$. The notation follows largely kulakov2019, kulakov2020. Define the supply and demand curves as a mapping of volumes to prices by $ \text{\textit{SUP}}_{WS} : \; (0, \infty) \rightarrow \left[P_\text{min}, P_\text{max}\right] $ and $ \text{\textit{DEM}}_{WS} : \; (0, \infty) \rightarrow \left[P_\text{min}, P_\text{max}\right]. $ Due to strict monotonicity the inverse ${\text{\textit{SUP}}_{WS}}^{-1}$ and ${\text{\textit{DEM}}_{WS}}^{-1} $ always exist. Hence, $\text{\textit{SUP}}_{WS}^{d,s}(q) = P$ is the supply or sell curve and $\text{\textit{DEM}}_{WS}^{d,s}(q) = P$ is the demand curve at the spot auction for delivery day $d$ and hour $s$ relating the volume $q$ ought to be sold/bought to the according price $P$. The inelastic demand in the wholesale market can be calculated by $DEM^{d,s}_\text{inelastic} = {DEM_{WS}^{d,s}}^{-1}(P_\text{min})$ where $P_\text{min} = -500$ EUR/MWh is the minimum price at the day-ahead auction epex2018. The transformed inverse supply curve can be written as:
As the curves are monotonic, ${SUP^{d,s}}^{-1}(z)$ also defines $SUP^{d,s}(q)$. As it is clearly visible in Figure (ref), the original equilibrium is reached at the point $P_\text{DA}^{d,s} = SUP^{d,s}(DEM^{d,s}_\text{inelastic})$. For the transformed curves it now holds that for the resulting clearing price, shifting $DEM^{d,s}_\text{inelastic}$ by some quantity $x$ equals shifting ${SUP^{d,s}}^{-1}(z)$ by $-x$, as the demand is perfectly inelastic.
Under the assumptions that the merit-order does not change significantly between day-ahead and intraday and that the transformed supply curve is a reasonable proxy for the merit-order, the the implied intraday demand and the slope coefficient for the merit-order can be derived. The first assumption is implicitly already made by kiesel2020a, kiesel2020b. The second assumption is discussed above. The last known 5-minute-interval-\gls{VWAP} before the start of the simulation is $\ensuremath{P^{d,s}_{\text{ID},0}}$. Under the \gls{MEH}, this price should reflect all changes to demand and supply. As all flexibility is already included in the supply curve, the implied intraday inelastic demand at $t=0$ can be calculated as $\text{\textit{DEM}}^{d,s}_{\text{implied}} = {\text{\textit{SUP}}^{d,s}}^{-1} (\ensuremath{P^{d,s}_{\text{ID},0}}).$ As already mentioned in Chapter (ref), the lower and upper price limits at the day-ahead auction are $\left[-500, 3000\right]$ EUR/MWh, while in the intraday market these are set to $\left[-9999, 9999\right]$ EUR/MWh. Hence, it might be possible that $\ensuremath{P^{d,s}_{\text{ID},0}}$ is outside the domain of ${\text{\textit{SUP}}^{d,s}}^{-1} (z)$. This case, however, does not occur in the dataset used in this paper. In the spirit of balardy2018 and kulakov2019, the measure for the elasticity $\textit{MO}^{d,s}_q$ is calculated as a finite central difference quotient of the transformed supply curve around $\text{\textit{DEM}}^{d,s}_{\text{implied}}$:
where $q = \{500, 1000, 2000\}$ MWh. It is defined in $\text{EUR}/\text{MWh}^2$ and is the steepness of the auction curves around the price level at $t=0$. Intuitively, it can be interpreted as the expected price change for a 1 MWh change in supply. In this paper, three values for $q$ are tested as there is some arbitrariness in choosing this value. In the literature, kulakov2019 choose $q = 100$ MWh, balardy2018 chooses $q = 500$ MWh. In this paper, slightly higher $q$ are selected to accommodate the fact that the standard deviation of $\ensuremath{\Delta\hat{W}_{\text{DA}, \text{ID}}^{d, s}}$ and $\ensuremath{\Delta\hat{S}_{\text{DA}, \text{ID}}^{d, s}}$ is roughly between 500 MW and 1300 MW (see Table (ref)). Values of $q < 500$ MWh thus might not catch the full range of volume changes occurring during the intraday trading. These volume changes in turn lead to movements along the merit-order. Figure (ref) (c) shows the intuition of the slope coefficient for $q = 2000$. \textcolor{Black}{Figure (ref) shows boxplots of the slope coefficients for $q = 2000$ MWh by the price level. Clearly, the slope increases with increasing price level, which is in line with the classic merit order model. However, we see the slope rising as well for small and negative prices. This observation is at odds with the common assumption of zero-marginal cost renewable production, which would imply a flat lower end of the merit order. However, many renewable assets are part of subsidy schemes, making their effective marginal costs negative. These assets are sold to the market even for negative prices, as long as the subsidy paid per produced MWh offsets negative selling prices.}
\FloatBarrier
Recall from the introduction that the price differences $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ follow certain distribution that we denote $G^{d,s}$ which is a mixture distribution
with the Dirac distribution $\delta_0$, the continuous distribution $F^{d,s}$ and the Bernoulli variable $\alpha_t^{d,s}$ with probability $\pi^{d,s}_t$. The two-stage approach is:
This introduces a dependence structure between the parameters of $\delta_0$ and $F$, as the probability $\pi^{d,s}$ is explained by past realisations of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ and $\ensuremath{\alpha^{d,s}_t}$. In the following three sections, the logistic model, the \gls{GAMLSS} model and the used benchmarks are presented.
\textcolor{Black}{The binary variable $\ensuremath{\alpha^{d,s}_t}$ will be modelled by a regularized logistic regression model tibshirani1996, meier2008 in the implementation of rglmnet using coordinate descent}. Generally, for a logistic model
for the Bernoulli variable $\alpha$ with probability $P(\alpha = 1) = \pi$, the \gls{LASSO} estimator $\hat{\mathbf{\beta}}^\text{LASSO}$ is given by
where $\lambda$ is a tunable shrinkage parameter. The corresponding log-likelihood $l$ is given by
where $\tilde{\mathbf{X}}$ is a standardisation of $\mathbf{X}$. The parameter $\lambda$ is optimised from an exponential grid of 100 values by choosing the minimum \gls{BIC}, i.e. $\lambda^\text{opt} = \text{arg} \; \text{min} \; \text{BIC} \left( \lambda \right)$ using the glmnet package by rglmnet.
Here, the logit function for $\ensuremath{\alpha^{d,s}_t}$ is explained by four components: the impact of past price differences, the time to maturity and weekday effects, fundamental variables such as \gls{RES} forecasts, outages and the slope of the merit-order, and a regression on averaged past $\ensuremath{\alpha^{d,s}_t}$. Intuitively, the probability of trades should rise with higher $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$, closer to delivery, with increasing wind and solar forecasts and with increased recent trading activity measured by past $\ensuremath{\alpha^{d,s}_t}$, but decrease on the weekends and the transition day Monday.
where $\bar{\alpha}^{d,s}_{t-j} = 1/j \cdot \sum_{i=1}^j \alpha_{t-i}^{d,s}$, the average of the last $j$ observed values of $\alpha^{d,s}_t$. \textcolor{Black}{This approach to transforming lagged values is similar to HAR-type models found in the field of financial econometrics}. A thorough description of the model is omitted here and can be found in ziel2020a. SAT($d$), SUN($d$), and MON($d$) are dummies for the weekday of $d$. TTD($t$) is a set of dummies for $t$. Accordingly, the model has more than 70 coefficients of which some tend to be highly correlated. To avoid problems with over fitting and multicollinearity, the model is estimated using the \gls{LASSO} of tibshirani1996. \textcolor{Black}{Note however, that for the one year training set used in this paper, we have ${365 \cdot T = 365 \cdot 31 = 11315}$ observations and are still in a setting where $n \gg{} p$. The number of observations is sufficiently larger than the number of parameters, hence identification is not an issue here.}
\textcolor{Black}{ This chapter briefly introduces the \gls{GAMLSS}-framework used to model
The \gls{GAMLSS} is an extension of the \gls{GAM} introduced by hastie1987, hastie1990. It allows to model not only the expected value of the a variable $Y \sim F$ but also the higher moments under a wide range of continuous and discrete distributions $F$. For an in-depth treatment we refer the reader to rigby2005, rigby2007, stasinopoulos2018gamlss and the manual of the R-package gamlss rigby2017. The notation in the following paragraphs follows the notation of aforementioned sources. }
\textcolor{Black}{ We first introduce the framework in an abstract notation. Following the mathematical formulation we will relate the abstract notation to the notation of the price differences. Let be $\boldsymbol{Y} = (Y_1, Y_2, ..., Y_n)$ be a vector of $i = 1, ..., n$ independent observations $Y_i$. The \gls{GAMLSS}-framework assumes that $Y_i$ have the probability density function
where each of the distribution parameters can be a smooth function of the explanatory variables. We denote as $\boldsymbol{\theta}_i = (\theta_{i,1} \theta_{i,2}, \theta_{i,3}, \theta_{i,4}) = (\mu_i, \sigma_i, \nu_i, \tau_i)$ the vector of $k=1,..., 4$ distribution parameters which are usually known as the location, scale and shape parameters $\theta_{i,k}$. For the distributions used in this paper, $\nu_i$ denotes the skewness and $\tau_i$ denotes the kurtosis. $\boldsymbol{\theta}$ is a matrix whose individual components have the indices $i$ and $k$. The vectors $\boldsymbol{\theta}_i$ and $\boldsymbol{\theta}_k$ are defined along the 2 axis of $\boldsymbol{\theta}$. Formally, we have
For each $k$, let $g_k(\cdot)$ be a known and monotonic link function that relates the distribution parameters $\boldsymbol{\theta}_k$ to the predictor $\boldsymbol{\eta_k}$. We consider the \gls{GAMLSS} model equation
where $\boldsymbol{X}_k$ is a $n \times J_k$ fixed design matrix and ${{\boldsymbol{\beta}}'_k = (\beta_{1,k}, \beta_{2,k},..., \beta_{J_k,k})}$ is a parameter vector of length $J_k$. The link functions $g_k$ ensure that the estimated distribution parameters fulfil the necessary assumptions concerning their support. To improve the robustness of the estimation, the following link functions are used:
We use $g_{\text{ident}}$ for the location parameters $\mu$ in both distribution assumptions, and additionally for the skewness parameter $\nu$ of the skewed t-distribution which also has support $(-\infty, \infty)$. $g_{\text{logident}}$ is introduced to avoid the exponential inverse for large estimates, thus improving the robustness of the estimation ziel2021m5, ziel2020a. We utilize it for all scale parameters $\sigma$. In addition, $g_{\text{logshift2}}$ is simply the natural logarithm shifted to 2 to preserve the condition $\nu > 2$ for the scale of the skewed $t$-distribution. For the remaining parameters we consider $g_{\text{log}}$. }
\textcolor{Black}{ rGamlssLasso extend the \gls{GAMLSS} framework to allow for regularized \gls{LASSO} estimation. As with the logistic model, we employ the \gls{BIC} to select the optimal shrinkage parameter $\lambda$. The adaptive \gls{LASSO} estimator $\boldsymbol{\beta}\ast_k$ is used. It is defined as
with the weights vector $\boldsymbol{\hat{w}}_k = 1 / \mid \hat{\beta} \mid^\gamma$. $\hat{\beta}$ denotes a root-$n$ consistent estimator such as ordinary least squares zou2006. }
\textcolor{Black}{ Let us now relate the abstract notation $Y_i$ to $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$. The distribution $F(\boldsymbol{\theta}_t^{d,s})$ is fitted to all ${\ensuremath{\Delta P^{d,s}_{\text{ID},t}} \mid \ensuremath{\alpha^{d,s}_t} = 1}$. The abstract index $i = 1, ..., N$ is replaced by the combination of the superscript index $d = 1, ..., D$ and the subscript index $t = 1, ..., T$. We fit 24 models each day, one for each delivery period $s$. The delivery periods are treated as independent. Thereby, we yield an estimated vector of four distribution parameters $\hat{\boldsymbol{\theta}}_{t}^{d,s} = (\widehat{\mu}_t^{d,s}, \widehat{\sigma}_t^{d,s}, \widehat{\tau}_t^{d,s} \widehat{\nu}_t^{d,s})$ and accordingly parameter estimates $\boldsymbol{\beta}^{d,s}_{k}$ that condition $\hat{\boldsymbol{\theta}}_t^{d,s}$ on our explanatory variables. Analogously, we can also define the the vector $\boldsymbol{\theta}_{k}^{d,s} = (\theta_{1, k}^{d,s}, ..., \theta_{T, k}^{d,s})$ along the time-axis $t$. We explain all moments of the distribution by the same set of explanatory variables. ziel2020a choose $F$ as Student's $t$-distribution. Here, we extend their choice to the skewed Student's $t$-distribution and Johnson's $S_U$ distribution. Both distributions have four parameters. A short description of the distributions used can be found in the Appendix (ref). }
\textcolor{Black}{For each distribution parameter $k$, the model for $\hat{\theta}_{t, k}^{d,s}$ reads:}
for $k = 1, 2, 3, 4$. The intercept is only included for $k \geq 2$ as the price differences are assumed to be centred around 0, as indicated by the summary statistics in Table (ref). The volatility is expected to rise with higher absolute past price differences, on weekends, and with higher \gls{RES} generation. The strong changes between day-ahead and intraday \gls{RES} and demand forecasts should also imply higher volatility. We expect the volatility to decrease with more recent trading activity measured by lagged $\alpha^{d,s}_t$. $\text{SIDC}(d, t)$ is a dummy variable taking the value 1 for $d \geq$ June 18, 2018 and $26 \leq t \leq 31$, indicating that the cross-country order books are closed. $f_\text{TTD}(t)$ models the non-linear impact of the time to delivery and takes the form $f_\text{TTD}(t) = 1 / \sqrt{T-t+1}$. It is a deterministic transformation of the variable $t$ and can thus be calculated ex-ante. As argued already in Section (ref), we expect a steeper merit-order regime to lead to higher price volatility. Similar expectations hold for the kurtosis, i.e. we expect a steep merit-order regime to lead to heavier tails of the distribution.
\FloatBarrier
Lastly, some simple benchmark models are introduced. Even though the main focus of this paper is on modelling the volatility and its influencing factor, simple benchmarks can serve as a valuable benchmark to identify potential areas for model improvement. \textcolor{Black}{As our study shares the conceptual set-up with ziel2020a, it is natural to employ similar benchmark models. Additionally, we compare our approach to classical time series methods such as \gls{ARIMA} models. The following section introduces these models in more detail.}
ziel2020a introduce six simple benchmark models to evaluate the value-added by more complex models, briefly described in the following. For the exact specification we refer the reader to their work.
\textcolor{Black}{The closeness of intraday electricity markets to traditional equity markets invites the use of classical time series models as benchmark. However, some attention to the unique time-structure of the intraday market is necessary: As already noted in Section (ref), for all delivery periods $s$ on day $d$, trading starts at 15:00 on $d-1$. The same delivery hour on two following delivery days can have overlapping intraday trading sessions. Thus, we cannot simply combine all trading sessions of a product, as it is done in equity markets. We therefore can only estimate our time series models on the price differences between the start of the trading period and the start of the simulation window. The GAMLSS-based approach does not suffer from this limitation as we learn the coefficients from past data of the simulation windows directly. The following paragraphs introduce the time series benchmark models formally.}
\textcolor{Black}{ The classic \gls{ARIMA}($p, k, q$) model is defined as follows:
is an ARIMA($p, k, q$) process with drift $\delta / (1 - \sum{} \varphi_i$). $L$ denotes the lag operator. A full treatment of ARIMA models can be found in e.g. shumway2017. We estimate the ARIMA($p, k, q$) models using the auto.arima() function in the forecast package rForecast. The function uses a stepwise approach to fit the lag order for $p$ and $q$ based on the \gls{BIC} and performs the KPSS unit-root tests to evaluate the integration order $k$. For each delivery period $d,s$, we fit the model on all 5-minute intervals between the start of trading on $d-1$, 15:00 and the start of the simulation period. The models are denoted as Auto.ARIMA. } \FloatBarrier
\FloatBarrier
\textcolor{Black}{We employ the well-known rolling window forecasting study design, which is common in energy price forecasting ziel2015, bunn2018, weron2018,janke2019, uniejewski2019, ziel2020a. This setting reduces the impact of structural breaks within the data and ensures a robust setting for the comparison of predictive performance using the \gls{DM}-test diebold2002, diebold2015. The scheme is visualized in Figure (ref).} We train one model for each delivery hour on 365 days of in-sample data and issue forecasts for the next delivery. Subsequently, the training data set is shifted forward by one day, the models are re-trained for each delivery hour and forecasts are issued for the next day and henceforth. Keeping the length of the training set constant we thus move through the test set. Our full data set ranges from January 2016 to August 2020, holding in total $N = 1618$ days. The training set length is fixed to $D = 365$ days. The test set holds $L = 1256$ days.
\FloatBarrier
Forecasts are issued for all delivery hours $s = 0, ..., 23$. For each delivery hour, the forecast consists of \textcolor{Black}{$j = 1, ..., M$ paths with $M = 1000$ paths of $t = 1, ..., 31$ steps.} \textcolor{Black}{Generally, let variables with superscript $[j]$ denote simulated values on path $j$ and hence $\ensuremath{{P}_{\text{ID},t}^{d,s,[j]}}$ denotes the simulation for step $t$ for delivery period $d,s$ in the path $j$.} The vector notation $\ensuremath{{\boldsymbol{P}}_{\text{ID}}^{d,s, [j]}} = (\ensuremath{{P}_{\text{ID},1}^{d,s,[j]}}, ..., \ensuremath{{P}_{\text{ID},31}^{d,s,[j]}})$ is frequently used in the chapter on error metrics. For the simulation of the paths, an algorithm similar to the recursive Euler-Maruyama-Scheme is used asmussen2007, ziel2020a. Each simulation starts 185 minutes and ends 30 minutes before the start of physical delivery. For each simulation step $t$ and path $j$, the boolean variable $\ensuremath{{\alpha}_{t}^{d,s,[j]}}$ is simulated $M$ times from the Bernoulli distribution $B(\ensuremath{\widehat{\pi}_{t}^{d,s,[j]}})$ and the price difference $\ensuremath{\Delta{}P_{\text{ID},t}^{d,s,[j]}}$ is sampled $M$ times from the distribution $F(\ensuremath{\widehat{\boldsymbol{\theta}}_{t}^{d,s,[j]}}) = F(\ensuremath{\widehat{\mu}_{t}^{d,s,[j]}}, \ensuremath{\widehat{\sigma}_{t}^{d,s,[j]}}, \ensuremath{\widehat{\nu}_{t}^{d,s,[j]}}$ and $\ensuremath{\widehat{\tau}_{t}^{d,s,[j]}})$. The price $\ensuremath{{P}_{\text{ID},t}^{d,s,[j]}}$ is then calculated as:
The algorithm is visualized in Figure (ref). The estimates for $\ensuremath{\widehat{\pi}_{t}^{d,s,[j]}}, \ensuremath{\widehat{\mu}_{t}^{d,s,[j]}}, \ensuremath{\widehat{\sigma}_{t}^{d,s,[j]}}, \ensuremath{\widehat{\nu}_{t}^{d,s,[j]}}$ and $\ensuremath{\widehat{\tau}_{t}^{d,s,[j]}}$ for the first step $t = 1$ are all equal, but begin to differ from $t \geq 2$ onwards as the paths develop individually. Therefore, the prediction matrix needs to be updated dynamically for each path and each step.
\FloatBarrier
The mean and median trajectory are evaluated using the \gls{RMSE} and \gls{MAE} respectively. For the probabilistic evaluation of the generated scenarios the \gls{ES}, \gls{CRPS} and the empirical coverage ratio are used. Additionally, the \gls{WS} is used to evaluate the coverage of an $(1-\alpha) \cdot 100 \%$-\gls{PI}. The \gls{ES}, \gls{CRPS} and the \gls{WS} are strictly proper scoring rules weron2018, gneiting2007, ziel2019b. To draw conclusions about the statistical significance of the difference in forecasting performance for each model, the \gls{DM}-test is used. All measures are widely employed in academia and practice.
Formally, the \gls{RMSE} and \gls{MAE} are defined as:
and:
For an $(1-\alpha) \cdot 100 \%$-\gls{PI} with the lower and upper bounds $L_t, U_t$ and prediction interval width $\delta_t = \hat{U}_t - \hat{L}_t$, the empirical \gls{CR} is defined as:
The $\text{WS}_t^{d,s}$ is defined as:
and aggregated as:
For both, \gls{CR} and \gls{WS}, the upper and lower bounds of the $(1- \alpha) \cdot 100 \%$ \gls{PI} are defined by the respective quantiles $\hat{L}_t^{d,s} = Q^{\alpha/2}_{j=1,...,M}(\ensuremath{{P}_{\text{ID},t}^{d,s,[j]}})$ and $\hat{U}_t^{d,s} = Q^{1-\alpha/2}_{j=1,...,M}(\ensuremath{{P}_{\text{ID},t}^{d,s,[j]}})$, where $Q^\tau_{j=1,...,M}(\ensuremath{{P}_{\text{ID},t}^{d,s,[j]}})$ denotes the $\tau$-th quantile of $M$ simulated $\ensuremath{{P}_{\text{ID},t}^{d,s,[j]}}$ prices. Comparing both, \gls{WS} and \gls{CR}, one can see how the \gls{WS} penalizes for an observation outside the interval and rewards the forecaster at the same time for a more narrow \gls{PI}. Contrary to the \gls{CR}, the \gls{WS} is a strictly proper evaluation measure weron2018.
The \gls{CRPS} gneiting2007, weron2018 is approximated by the \gls{PB}
for a dense equidistant grid of probabilities $\mathcal{T} = \{0.01, ... 0.99\}$ of size $R = 99$. $\text{PB}^{d,s}_{t,\tau}$ denotes the pinball loss for probability $\tau$. The formula is given by:
The overall \gls{CRPS} is calculated by the average:
The pinball loss is also used to evaluate the performance of different models in specific quantile levels. For this reason, the \gls{PB} is aggregated as follows:
To measure the quality of the generated paths, ziel2020a propose the \gls{ES}. It is a generalisation of the \gls{CRPS} for two dimensions. Thereby, not only the approximation of the marginal distribution is evaluated, but the generated multivariate distribution gneiting2007, ziel2019b:
The average yields the overall energy score for each model:
\textcolor{Black}{The aforementioned measures provide insight in the accuracy of different forecasting models. To evaluate the statistical significance of the difference in forecast accuracy of two models $A$ and $B$, the \gls{DM}-test diebold2002, diebold2015 is routinely employed in the field of energy price forecasting weron2018, ziel2018day, janke2019. It originally stems from the field of point forecasting, however diebold2015 notes that the test is agnostic to the scoring rule used to evaluate forecasts. Hence, using strictly proper probabilistic scoring rules, such as the \gls{CRPS} and \gls{ES} loss, the \gls{DM} test can be applied to probabilistic forecasts as well diebold2015, weron2018. Following ziel2020a and ziel2018day, the \gls{DM}-test is employed in a multivariate fashion. Hence, let $\ensuremath{\boldsymbol{L}_{A}^{d}} = (\ensuremath{L_{A}^{d, 1}}, ..., \ensuremath{L_{A}^{d, S}})$ and $\ensuremath{\boldsymbol{L}_{A}^{d}} = (\ensuremath{L_{B}^{d, 1}}, ..., \ensuremath{L_{B}^{d, S}})$ denote the out-of-sample loss vectors for model $A$ and $B$ for day $d$ and delivery period $s$ of length $N$. For models $A, B$, the $N \times S$ vector of losses are reduced to an $N \times 1$ vector by taking the 1-norm. The difference between both is the loss differential used in the DM-test.
For example, for the \gls{ES} and the Naive model, the loss vector $\boldsymbol{L}_{\text{Naive}}^d = (\text{ES}_{\text{Naive}}^{d, 1}, ..., \text{ES}_{\text{Naive}}^{d, S})$.}
\textcolor{Black}{We test the loss differential series for stationarity using the augmented Dickey-Fuller (ADF) test dickey1979, dickey1981 and reject the $H_0$ of unit root at the 5% significance level for all loss differential series.} harvey1997 propose the usage of the $t$-distribution with $\nu = N-1$ degrees of freedom rather than the normal distribution, as well as the introduction of a bias correction. Formally, the corrected test statistic is defined as $t_\text{DM}^{\text{HLN}, h=1} = \sqrt{\frac{N+3}{N}} \cdot t_\text{DM} \sim t(0, 1, N-1) $ under the $H_0$, where $N$ is the length of the loss differential series $\Delta L^d_{A,B}$ and $h$ denotes the forecast horizon. \textcolor{Black}{The standard deviation is computed using an autocorrelation-consistent estimator.} For each model pair, two one-sided \gls{DM}-tests are computed. The first test has the $H_0$ that the forecasts of model $A$ are significantly better than the forecasts of model $B$. For the second test, the $H_0$ is that the forecasts of model $B$ are significantly better than the forecasts of model $A$. These tests are complimentary. We use the implementation in the R-package forecast hyndman2008, rForecast.
\FloatBarrier
The following chapter presents the results of the forecasting study. It is split into two parts: First, we show the error metrics for the out-of-sample analysis. Additionally, we show the in-sample coefficients for the model using Johnson's $S_U$ distribution.
\FloatBarrier
First, the aggregate error statistics will be presented, followed by the scoring rules considering the marginal fit relative to the time to delivery and the quantile range $\mathcal{T}$. Statistical significance is evaluated using the Diebold-Mariano test.
The Naive performs best in terms of RMSE and MAE, while Mix.JSU performs best across the probabilistic evaluation using the CRPS and ES scoring rules. Its superior performance in terms of ES is statistically significant according to the \gls{DM}-test. \textcolor{Black}{The GAMLSS-based model assuming the skew-$t$ distribution however shows a very poor performance for hour 6, which yields an overall poor performance. For this delivery period, we can trace the high error back to outliers and extreme $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ larger than 2000 EUR/MWh on March 11th, 2020. This indicates that Johnson's $S_U$ is more robust towards outliers. An investigation of the loss time series for Mix.JSU and Mix.SST shows the deteriorating forecasting performance of the Mix.SST after March 11th, 2020 clearly (see Figure (ref) in Appendix (ref)). The Auto.ARIMA performs surprisingly bad in terms of the RMSE and MAE and somewhat better in terms of the CRPS and ES. With respect to the other benchmark models, we see an overall mixed performance. We note a worse performance for the benchmark models for the probabilistic measures CRPS and ES compared to the \textbf{Naive} and \textbf{Mix.JSU}. }
Figure (ref) shows the ES relative to the delivery hours $s$. Relative to the Naive, the \gls{GAMLSS}-based models show an improved forecasting performance in the peak hours (Plot (ref) b). The error of the RW.N and Mix.SST models explodes for hour $s=6$.
The \gls{PB} over $\mathcal{T} = 0.01, ..., 0.99$ is shown in Figure (ref). Again, subfigure (a) represents absolute values and (b) depicts all models relative to the Naive. All models show similar performance in the central quantiles, as already indicated by the very close values for the \gls{MAE}. Relative to the Naive, most other benchmark models show worse performance in the tails of the distribution. The Mix.JSU shows an improved modelling of the tails compared to the Naive. The Mix.SST again shows a weak performance given by its sensitivity to outliers. The development of the \gls{CRPS} throughout the simulation window is shown in Figure (ref). Again, (a) shows absolute values while (b) shows the error relative to Naive. The CRPS is rising through the simulation window, especially for the last 60 to 30 minutes of trading. \textcolor{Black}{The relative error of most models towards the \textbf{Naive} decreases throughout the simulation window, however, it increases for the \textbf{Auto.ARIMA}. This might indicate that learning the model parameters of past trading sessions can be beneficial compared to learning the parameters only from the trading session of interest, before the start of the simulation, as market behaviour changes throughout the session.}
The results are largely confirmed as statistically significant by the Diebold-Mariano-Test. Figure (ref) shows the $p$-values for the pairwise \gls{DM}-tests for the \gls{ES} and \gls{CRPS}. The lower the $p$-value, the more significant is the difference in the forecasting performance, which implies that the model on the column (or $x$-axis) outputs superior forecasts than the model on the row (or $y$-axis). Generally, the $p$-values for the \gls{CRPS} and \gls{ES} are rather close. This makes sense, as a good coupling to the path's distribution should be closely related to a good fit on the marginal distribution. The other way, however, is not necessarily true. Inside the group of the benchmark models, the Naive model is confirmed as the superior model as it yields significantly better forecasting performance than all other benchmark models. \textcolor{Black}{The Auto.ARIMA is significantly better as the RW.N only.} The Mix.JSU yields significantly better forecasting performance than all other models in terms of the ES. The Mix.SST yields significantly worse forecasting accuracy in terms of both CRPS and ES than all other models, which is expected given the results shown in Figures (ref) to (ref).
\FloatBarrier
Given the strong probabilistic forecasting performance of the Mix.JSU we turn to an in-sample analysis of the estimated coefficients. \textcolor{Black}{Compared to black-box deep learning algorithms, the parametric \gls{GAMLSS} framework used in this paper allows for explainable machine learning by quantitatively and qualitatively analysing the estimated coefficients. Hence, we can gain further insight in the driving variables for all distribution parameters.} Tables (ref) to (ref) present the estimated \textcolor{Black}{scaled} coefficients for $d$ = January 22nd, 2017, the first out-of-sample day. \textcolor{Black}{Scaled coefficient correspond to mean-variance scaled inputs. Hence, the coefficients are hence unit-free and can be compared in the magnitude.} The background colouring indicates the share of non-zero estimates for the whole out-of-sample data set. Green indicates that few estimates are set to zero by the sparsity property of the LASSO, the darker the red, the more estimates are set to zero.
For $\mu$, only the first lag of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ shows more than a couple non-zero values for the first day of the test set. This variable yields non-zero estimates as well across the test set for the late morning to afternoon peak hours. This result is similar to the findings of ziel2020b, who find the most recent price to be among the most important features for forecasting the $\text{ID}_3$ as well as with the results of kiesel2020a, kiesel2020b, who find that lagged prices are an important predictor. The fact that other fundamental and trading related information, especially intraday forecast changes, do not yield additional predictive power suggests that this information is contained in the price already. These results support the notion of weak-form market efficiency already indicated by ziel2020b and kuppelwieser2021.
For the volatility $\sigma$, we present coefficients in similar fashion in Table (ref). For the first day of the test set, we yield non-zero estimates for the coefficients for the merit-order slope, for the intercept and for the transformed time to delivery. For a few hours, the coefficient for lagged values of $\ensuremath{\alpha^{d,s}_t}$ has a negative non-zero estimate as well. The large and positive coefficients for the merit-order slope confirm our initial assumption that the shape of the merit-order is a driving factor for the volatility in intraday markets. Intuitively, this is derived from the observation that on a steep merit-order, a slight change in supply or demand has a higher impact on the price than in a flat regime. Moving this effect from a threshold variable for the size of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ to the volatility parameter of the distribution of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ thus generalizes the results of kiesel2020a, kiesel2020b. Contrary to baule2021, we find little predictive power for the spread between spot and intraday price as well as for the fundamental forecasts and their intraday changes for the volatility. Remember that the coefficients in Table (ref) correspond to January 22nd, 2017, well before the introduction of SIDC. Thus, the SIDC variable is zero for this training period. We show the evolution of the estimated coefficient across the rolling training set in Figure (ref). After the launch of SIDC on June 13, 2018, the dummy is first included in the rolling training set. A sizeable positive estimate is visible, i.e. the volatility rises after gate closure of the cross-border shared order books 60 minutes before delivery. The effect is the strongest in 2019 and 2020 for the morning and afternoon peak hours and less clear for the solar peak hours around noon. Our findings are consistent with ziel2020a and contradict kath2019, who finds no evidence of rising volatility due to SIDC.
For the skewness parameter $\nu$ none of the variables apart from the intercept yield non-zero coefficient estimates. The intercept is slightly negative in the night hours and positive in the morning and afternoon hours. We conclude that the intraday price returns do not exhibit any strong skewness within the individual trading sessions.
Lastly, we turn to Table (ref) giving the estimated coefficients for the kurtosis parameter $\tau$. We find a negative impact of lagged $\ensuremath{\alpha^{d,s}_t}$ and a positive impact of lagged $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$. Thus, we expect the distribution of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ to be lighter-tailed if there has been no trade in the preceding 15 minutes of trading. On the other hand, large absolute price changes in the previous 15 minutes of trading increase $\tau$ and thus the heaviness of the tails. The impact of lagged $\ensuremath{\alpha^{d,s}_t}$ and $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$ is more pronounced during the night hours. During the day hours, there are some none-zero estimates for wind and solar forecasts. This is consistent with kiesel2020b's finding that the behaviour of night contracts is more driven by trading-related variables than fundamentals. For $\tau$, we find a negative impact of the merit-order slope parameter. This implies that a steep-merit leads to heavier tails for the distribution of $\ensuremath{\Delta P^{d,s}_{\text{ID},t}}$. Thus, if the merit-order is steep, not only the volatility level is elevated, but also the likelihood of spikes is higher. Lastly, we find that with decreasing time to delivery the heaviness of the distribution's tails decreases.
\FloatBarrier
This paper develops a simulation-based forecasting model for the intraday price process in the last three hours of each product's trading window. We expand the key work of ziel2020a in four dimensions by (i) investigating distributions with potential skewness and modelling all moments explicitly, (ii) adding intra-daily forecast updates and (iii) a novel measure for the merit-order slope, derived from day-ahead auction curves, and (iv) employing a regularized estimation using the \gls{GAMLSS}-\gls{LASSO} for all distribution moments.
Our results are two-fold: \textcolor{Black}{First, we show that the proposed method is able to generate high quality ensembles for the intraday markets, whose predictive performance is significantly better than benchmark models such as random walk or \gls{ARIMA}-type processes on a wide range of probabilistic scoring rules. The improvement in accuracy is especially distinct in the tails of the predictive distribution. Thus, our results can be applied directly to trading problems as proposed by serafin2022 or plugged into any optimization method relying on accurate sampling methods.} Second, the GAMLSS framework's explicit traceability and the regularized estimation allows to draw conclusions on the impact of explanatory variables. Qualitatively, our results for the expected value of the intraday return distribution imply weak-form efficient markets, as the inclusion of additional variables does not improve the prediction of the expected value significantly. Additionally, we find evidence for a merit-order effect in the volatility and kurtosis of the return distribution. A steep merit-order regime leads to higher volatility and heavier tails. What is more, we find that the volatility rises with decreasing time to delivery and rises with the closure of the pan-European order book sharing (SIDC). On the other hand, the kurtosis is driven by trading-related variables such as trade events and lagged prices. We find however, that the skewness is close to zero for all hours, and that none of the analysed variables show predictive power.
This paper's result opens several new research strings: the models used can be improved by the inclusion of cross-product effects and neighbouring products as additional input variables. However, due to the structure of intraday markets with parallel and overlapping trading sessions, this task is non-trivial. A second interesting research avenue is the relationship between trading volume, liquidity and volatility in intraday markets. Further research is also needed to better understand the impact of fundamental variables for modelling the volatility, kurtosis and skewness of the distribution of intraday price returns. The influence of the merit-order shape as explanatory variable for the volatility warrants further research into its modelling for short-term markets.
\FloatBarrier
This paper is based on research conducted during a joint project of Simon Hirsch and Statkraft Trading GmbH. Simon Hirsch is grateful to Statkraft, especially Patrick Otto, Dr.\ Konstantin Wiegandt and Dr.\ Daniel Gruhlke for the support received while writing his thesis. The authors are grateful to energy & meteo systems GmbH for providing the forecasts used in the paper. The views and opinions expressed in this paper are the author's own and do not reflect the views of Statkraft Trading GmbH or energy & meteo systems GmbH. The authors are grateful to helpful discussions at the 30. GEE Doctoral Workshop, Essen, 2022.
Due to the commercial nature of production forecasts the dataset remains confidential and cannot be shared.
Simon Hirsch is employed by Statkraft Trading GmbH. The authors declare no conflict of interest.