Simulation-based Forecasting for Intraday Power Markets: Modelling Fundamental Drivers for Location,Shape and Scale of the Price Distribution

Abstract

During the last years, European intraday power markets have gained importance for balancing forecast errors due to the rising volumes of intermittent renewable generation. However, compared to day-ahead markets, the drivers for the intraday price process are still sparsely researched. In this paper, we propose a modelling strategy for the location, shape and scale parameters of the return distribution in intraday markets, based on fundamental variables. We consider wind and solar forecasts and their intraday updates, outages, price information and a novel measure for the shape of the merit-order, derived from spot auction curves as explanatory variables. We validate our modelling by simulating price paths and compare the probabilistic forecasting performance of our model to benchmark models in a forecasting study for the German market. The approach yields significant improvements in the forecasting performance, especially in the tails of the distribution. At the same time, we are able to derive the contribution of the driving variables. We find that, apart from the first lag of the price changes, none of our fundamental variables have explanatory power for the expected value of the intraday returns. This implies weak-form market efficiency as renewable forecast changes and outage information seems to be priced in by the market. We find that the volatility is driven by the merit-order regime, the time to delivery and the closure of cross-border order books. The tail of the distribution is mainly influenced by past price differences and trading activity. Our approach is directly transferable to other continuous intraday markets in Europe.

1. Introduction

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 Renewable Energy Sources (RES) the opportunity to balance forecast errors arising after the day-ahead auction until five minutes before the beginning of the delivery period (Koch and Hirth, 2019). Increasingly, this balancing action is taken over by algorithmic trading strategies, for which reliable short-term price and volatility forecasts are necessary. At the same time, our results shed light on the influence of idiosyncratic features of intraday markets such as Single IntraDay Coupling (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 (Janke and Steinke, 2019; Narajewski and Ziel, 2020a; Uniejewski et al., 2019; Cramer et al., 2022) and the markets’ driving fundamentals has, so far, focused on modelling the impact of renewable forecast (Ziel, 2017; Kath, 2019; Pape et al., 2016; Gürtler and Paulsen, 2018; Balardy, 2022) and forecast errors (Ziel, 2017; Kulakov and Ziel, 2020; Kuppelwieser and Wozabal, 2021). A different strand of literature emerged around modelling of the merit-order effect for price changes and price elasticity (Kiesel and Paraschiv, 2017; Kremer et al., 2021, 2020; Kulakov and Ziel, 2019; Balardy, 2022). To the best knowledge of the authors, only Narajewski and Ziel (2020b) and Baule and Naumann (2021) 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 Generalized Additive Models for Location, Shape and Scale (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 analyze 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 Narajewski and Ziel (2020b) to develop a simulation-based probabilistic forecasting model for the path of the five-minute-volume-weighted price $P$ between 185 and 30 minutes before delivery. Instead of directly modelling the price $P$ as it is common in dayahead forecasting and is done in other forecasting studies on the intraday market (Uniejewski et al., 2019; Janke and Steinke, 2019), the first differences $\mathrm{\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 Narajewski and Ziel (2020b), we assume the first differences $\mathrm{\Delta }P$ follow a mixture distribution of the Dirac distribution ${\delta }_0$ with an atom at 0 and a continuous distribution $F$ . In Narajewski and Ziel (2020b), the latter is assumed to be $t$ -distributed without any detailed justification, except the observation that $\mathrm{\Delta }P$ tends to be heavy tailed. This manuscript extends the approach in four dimensions:

1. We study in more detail the distribution assumption for $F$ . We investigate more distributions, including those with potential skewness. The distribution for $F$ is selected from the skew-t and the Johnson’s $S_U$ . These distributions have been used successfully to model asset returns and/or have been applied to forecasting in energy markets (Serinaldi, 2011; Gianfreda and Bunn, 2018; Bunn et al., 2018).

2. One natural starting point to improve the models is adding intra-daily updated forecasts for wind and solar generation to reflect the information set available to market participants better. High forecast updates indicate that market participants with RES assets need to solve larger positions in the intraday market. Capturing this effect aims at improving the modelling of the location and volatility of $\mathrm{\Delta }P$ . The findings will be discussed in the light of the market efficiency hypothesis indicated by Narajewski and Ziel (2020a); Kuppelwieser and Wozabal (2021). Additionally, the variance of individual forecast updates is used to model the volatility of $\mathrm{\Delta }P$ . This is motivated by the assumption that uncertainty regarding the weather situations is transmitted to the price formation.

3. A measure for the intraday price elasticity is derived from the day-ahead auction curves, which are used as proxy for the intraday merit-order and included in the modelling of the volatility (Kulakov and Ziel, 2019, 2020). This is motivated by Kremer et al. (2021, 2020), who show that the price impact of forecast errors depends on whether the market is in a flat or steep merit-order ’regime’ and the work of Balardy (2022), who shows that the intraday bid-ask spread can be explained by the price elasticity derived from the day-ahead auction curves. Additionally, the predictive power of the day-ahead level of outages and the change in planned and unplanned outages between day-ahead and intraday is tested.

4. As in Narajewski and Ziel (2020b), we utilize the GAMLSS framework for parameter estimation. However, we allow for a more flexible parameter training approach utilizing an automatic variable selection for all distribution moments by utilising the adaptive Least Absolute Shrinkage and Selection Operator (LASSO) estimation technique for the GAMLSS model as well. This procedure was similarly used in Ziel (2021) in the context of the M5 forecasting competition.

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 Diebold-Mariano (DM)-test (Narajewski and Ziel, 2020b; Ziel and Berk, 2019; Nowotarski and Weron, 2018; Diebold and Mariano, 2002). 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 $\mathrm{\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 $\mathrm{\Delta }P$ . This result supports the notion of weak-form efficient markets already indicated by Narajewski and Ziel (2020a); Kuppelwieser and Wozabal (2021). We find evidence for a merit-order effect in the volatility and kurtosis of the distribution $\mathrm{\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 flexibility using Monte-Carlo option valuation. Additionally, the explicit modelling of the volatility provides a starting point to introduce time-varying volatility to stochastic optimization models for market making, position solving and sequential bidding (von Luckner et al., 2017; Glas et al., 2020; Aïd et al., 2016; Kath and Ziel, 2020; Kraft et al., 2023).

The remainder of this paper is structured as follows: Section 2 gives a short introduction to the structure of the German short-term power markets. Section 3 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 4 introduces the used models. The forecasting study design and scoring rules are discussed in Section 5. Finally, Sections 6 and 7 present the results and conclude this paper.

2. Structure of the German Power Market

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 $d$ and the delivery hour as $s$ for $s=0,\dots ,S$ and $S=23$ . Times are usually expressed in local time unless otherwise noted. ¹ 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 1 shows the time line of the German shortterm markets, for more details see Viehmann (2017). 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.

Figure 1

Daily procedure in the German short-term power markets (based on EPEX SPOT SE, 2018; Narajewski and Ziel, 2020b; Nordpool SE, 2018).

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 $P_{\mathrm{DA}}^{d,s}$ . 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. EPEX SPOT SE (2020a) publishes aggregate curves together with the official market results. The lower price level is set to -500 EUR/MWh, the upper level is set to 3000 EUR/MWh .

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 SIDC system, formerly known as Cross Border Intraday (XBID), starts in Germany, Denmark, Netherlands, Norway and Poland. ² 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. SIDC went live on June 18, 2018 (Nordpool SE, 2018). Thirty minutes before delivery, the Germany-wide order book closes and trading resumes in local (control zone) products up to five minutes before delivery. Note that all open delivery periods are traded in parallel. The market price limits are at $\pm 9999$ 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 (EPEX SPOT SE, 2020b, 2018; Viehmann, 2017).

3. Data and Exploratory Analysis

3.1 Intraday Trade Data

On the intraday market, trading happens continuously. Hence, the transactions are irregularly spaced and need to be aggregated. The following paragraphs and Figure 2 give a brief overview of the aggregation. A detailed description can be found in Appendix A.1. Trade data is obtained from EPEX SPOT SE (2020c). The data consists of all hourly trades on the continuous intraday market between January 1st, 2016 and July 31st, 2020.

Figure 2

Overview of the data preparation. The first figure shows the raw trade data. The color refers to the traded volume. Below, the 5 -minute Volume-Weighted Average Price (VWAP) $P$ and their first differences $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ are shown. Lastly, ${\alpha }_t^{d,s}$ is shown. The plots show the data for $d=$ January 1st, 2016 for delivery period $s=12$ in the last 6 hours to delivery.

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 $P_{\mathrm{ID}.,t}^{d,s}$ , where $t$ denotes the 5-minute interval (see panel 2 in Figure 2). We then take first differences $\mathrm{\Delta }{P}_{\mathrm{ID},t,t}^{d,s}=P_{\mathrm{ID},t}^{d,s}-P_{\mathrm{ID},t-1}^{d,s}$ (see panel 3). Lastly, we define a boolean variable ${\alpha }_t^{d,s}$ , 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 Narajewski and Ziel (2020b,a) and Serafin et al. (2022).

Figure 3 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 $\mathrm{I}{\mathrm{D}}_3$ period already. Additionally, Table 1 presents summary statistics for ${\alpha }_t$ and $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d.s}\mid {\alpha }_t=1$ converge. We can thus identify two levels of time-varying behaviour of ${\alpha }_t^{d,s}$ , 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.

Figure 3

Share of no-trade events over time to delivery and delivery hour for the initial oneyear training set from January 1st to December 31st, 2016. Colour corresponds to the $z$ -axis. Low values correspond to few no-trade events, high values correspond to many no-trade events.

Table 1

Summary statistics for ${\alpha }_t^{d,s}$ and $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ for 185 to 30 minutes before delivery. $Q_{\tau }$ denotes the empirical $\tau ·100\mathrm{\%}$ -quantile. Note that all $\mathrm{\Delta }{P}_{\mathrm{ID}.t}^{d,s}\mid {\alpha }_t^{d,s}=0$ are 0 by definition.

Variable	Summary Statistics	$2016$	$2017$	$2018$	$2019$	$2020$
${\alpha }_t^{d,s}$	Count $(n)$	257424	252960	268584	267096	158472
	Mean $(\mu )$	0.82	0.87	0.95	0.98	0.98
	Std. $(\sigma )$	0.39	0.34	0.23	0.14	0.15
$\mathrm{\Delta }{P}_{\mathrm{II},,t}^{d,s}\mid {\alpha }_t^{d,s}=1$	Count $(n)$	209995	219258	253931	261721	154597
	Mean $(\mu )$	0.02	0.02	-0.00	-0.00	-0.01
	Std. $(\sigma )$	1.57	2.06	1.94	2.15	8.13
	MAD	0.94	1.11	1.03	0.91	1.10
	IQR	1.17	1.33	1.23	1.01	1.07
	Skewness	1.07	1.94	8.02	12.25	0.65
	Kurtosis	131.50	322.14	902.99	2544.57	64968.44
	Min	-76.38	-127.70	-108.22	-214.16	-2161.80
	$Q_{0.10}$	-1.40	-1.61	-1.52	-1.29	-1.45
	$Q_{0.25}$	-0.57	-0.65	-0.62	-0.51	-0.54
	$Q_{0.50}$	0.00	0.00	0.00	-0.01	0.00
	$Q_{0.75}$	0.60	0.68	0.61	0.50	0.53
	$Q_{0.90}$	1.43	1.64	1.49	1.26	1.40
	Max	82.16	119.88	183.94	226.26	2164.23

The mean and median values of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ . The minima and maxima increase throughout the data set, while the $5\mathrm{\%}$ respectively $95\mathrm{\%}$ and the $10\mathrm{\%}$ respectively $90\mathrm{\%}$ 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 $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ rises fourfold between 2019 and 2020, while staying roughly constant before. The more robust dispersion measures median absolute deviation (MAD) and the interquartile range (IQR) support this notion.

Figure 5 plots histograms for $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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. The center bar shows the relative weight of 5 -minute intervals without trades (i.e. ${\alpha }_t^{d,s}=0$ ) and 5 -minute intervals with at least one trade (i.e. ${\alpha }_t^{d,s}=1$ ), but small or no price changes. As visible already in Figure 3, the share of 5 -minute intervals with ${\alpha }_t^{d,s}=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. Additionally, Figure 6 (a) plots the Pearson autocorrelation coefficients $r$ for $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 6 (b) shows the according $p$ -values for the test statistic $r·\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\mathrm{\%}$ -level and for only $15\mathrm{\%}$ 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 19 and 20 in Appendix B).

The relationship between the realized variance and the time to delivery in the initial training data is shown in Figure 4. 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. 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.

Figure 4

Volatility development over time to delivery and delivery hour for the initial training one-year set from January 1st to December 31st, 2016. Colour corresponds to the $z$ -axis. Note the inverted $y$ -axis (Time to delivery) compared to Figure 3 to ease visualization.

Figure 5

Histograms of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ for $=4,12,20$ in the initial one-year training set.

Figure 6

Autocorrelation of $P_{\mathrm{ID},t}^{d,s}$ per trading window for lags 1 and according $p$ -values. The first heat maps show the size of the correlation coefficient by delivery day $d$ and hour $s$ , second shows according $p$ -values. Lags 2 and 3 can be found in Appendix 9 .

To analyze 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 2. 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 Löhndorf and Wozabal (2022), who aggregate trades in the intraday market on a 1-hour grid and report similar results for the ADF-test at the $10\mathrm{\%}$ -level, we conclude that the price changes in the intraday market are stationary.

Table 2

Aggregate results of the augmented Dickey-Fuller unit-root tests. Tests are applied to each simulation window individually. The table reports the share of simulation windows where the test result implies stationarity, i.e. rejection of the $H_0$ of a unit root.

	YearSignificance	$2016$	$2017$	$2018$	$2019$	$2020$
Variable	$\alpha =0.01$	0.675	0.658	0.633	0.596	0.555
$\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$	$\alpha =0.05$	0.750	0.744	0.717	0.684	0.641
$P_{\mathrm{ID},t}^{d,s}$	$\alpha =0.01$	0.059	0.055	0.055	0.059	0.088
	$\alpha =0.05$	0.128	0.125	0.122	0.117	0.157

3.2 Renewables Forecasts and Outages

Intra-daily updated renewable production forecasts used in this paper are provided by Statkraft Markets and generated by energy & meteo systems GmbH (2020). Day-ahead demand / system load forecasts are obtained from ENTSO-E (2021). 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 ${\hat{W}}_v^{d,s},{\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 ${\hat{D}}_{\mathrm{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. ³ 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 ${\hat{W}}_{\mathrm{DA}}^{d,s}$ and ${\hat{S}}_{\mathrm{DA}}^{d,s}$ . Secondly, the newest forecast available before the start of the simulation, i.e. at which $v\ge b(d,s)-185$ holds, is denoted as the intraday forecast ${\hat{W}}_{\mathrm{ID}}^{d,s}$ and ${\hat{S}}_{\mathrm{ID}}^{d,s}$ .

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 $\mathrm{\Delta }{\hat{W}}_{v1,v2}^{d,s}={\hat{W}}_{v_2}^{d,s}-{\hat{W}}_{v_1}^{d,s}$ with $v_2$ being the newer forecast.

• The day-ahead to simulation forecast update is defined as: $\mathrm{\Delta }{\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s}={\hat{W}}_{\mathrm{ID}}^{d,s}-{\hat{W}}_{\mathrm{DA}}^{d,s}$ . The symmetry of the impact of forecast errors on $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ is a disputed topic in the literature (Ziel, 2017; Kremer et al., 2021, 2020) and has, so far, not been explored for the volatility of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ . To address this issue and test for possible asymmetric effects, $\mathrm{\Delta }{\hat{W}}_{\mathrm{DA}}^{d,s}$ and $\mathrm{\Delta }{\hat{S}}_{\mathrm{DA},\mathrm{ID}}^{d,s}$ are split in positive and negative updates: $\mathrm{\Delta }{\hat{W}}_{\mathrm{DA}}^{d,s,+\mathrm{ID}}=\max (\mathrm{\Delta }{\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s},0)$ and $\mathrm{\Delta }{\hat{W}}_{\mathrm{DA}}^{d,s,\mathrm{ID}}=\left|\min (\mathrm{\Delta }{\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s},0)\right|$ .

• The standard deviation of the forecast updates should reflect the uncertainty about the weather situation. Highly volatile forecast updates are mirrored to starkly changing positions in renewable energy asset portfolios and should thus exercise an influence on the volatility of the price process due to quickly changing demand and supply. Let

$${\mathcal{V}}^{d,s}=[d-1,12:00\ge v\ge b(d,s)-185]$$

be the set of all forecast versions received between the day-ahead auction and the start of the simulation. The difference between two consecutive forecast versions $v,v-1$ is denoted as $\mathrm{\Delta }{\hat{W}}_{v-1,v}^{d,s}$ . Then ${\sigma }_{\mathrm{DA},\mathrm{ID}}^{d,s}(\mathrm{\Delta }\hat{W})={\sigma }_{v\in \mathcal{V}d,s}\left(\mathrm{\Delta }{\hat{W}}_{v-1,v}^{d,s}\right)$ denotes the standard deviation of all differences between two consecutive forecasts received between the day-ahead auction and the start of the simulation. It is worth noting that due to the schedule of the day-ahead and intraday markets, $\mathcal{V}$ is larger for later delivery hours, thus more forecast versions are considered for ${\sigma }_{\mathrm{DA},0}^{d,s}(\mathrm{\Delta }\hat{W})$ for later delivery hours.

Analogously, $\mathrm{\Delta }{\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s},\mathrm{\Delta }{\hat{S}}_{\mathrm{DA}A,\mathrm{ID}}^{d,s,+},\mathrm{\Delta }{\hat{S}}_{\mathrm{DA},\mathrm{ID}}^{d,s,-}$ and ${\sigma }_{\mathrm{DA},\mathrm{I}}^{d,s}(\mathrm{\Delta }\hat{S})$ are defined for the solar production forecasts. Panels (a)-(c) in Figure 7 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.

Figure 7

Day-ahead RES production forecasts, forecast updates and their standard deviation and the aggregated outage notifications.

Under the Regulation on wholesale Energy Market Integrity and Transparency (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 (European Commission, 2011; Lazarczyk and Le Coq, 2018; Agency for the Cooperation of Energy Regulators (ACER), 2020). We retrieve unavailability messages from the European Energy Exchange (EEX) AG (2020) 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: $\mathrm{\Delta }{O}_{\mathrm{DA}}^{d,s,\mathrm{ID}}=O_{\mathrm{ID}}^{d,s}-O_{\mathrm{DA}}^{d,s}$ . Afterwards, the difference $\mathrm{\Delta }{O}_{\mathrm{DA},0}^{d,s}$ is split into planned and unplanned outages denoted as $\mathrm{\Delta }{O}_{\mathrm{DA},\mathrm{ID}}^{d,s,\mathrm{panned}}$ and $\mathrm{\Delta }{O}_{\mathrm{DA},\mathrm{ID}}^{d,s,\mathrm{sunned}}$ . Figure 7 (d) plots the aggregated outage data; Table 3 gives summary statistics.

Table 3

Summary Statistics for wind, solar and demand forecasts, and the reported outages. Day-ahead to intraday forecast changes and the standard deviation of versioned updates. Indices $d,s$ and $t$ omitted for better readability. Values in MW.

	Wind			Solar			Demand	Outages
	$\hat{W}$	$\mathrm{\Delta }\hat{W}$	$\sigma (\mathrm{\Delta }\hat{W})$	$\hat{S}$	$\mathrm{\Delta }\hat{S}$	$\sigma (\mathrm{\Delta }\hat{S})$	$\hat{D}$	$O$	$\mathrm{\Delta }{O}_{\mathrm{planned}}$	$\mathrm{\Delta }{O}_{\mathrm{unplanned}}$
Count ( $n$ )	38856	38856	38856	38856	38856	38856	38856	38856	38856	38856
Mean $(\mu )$	12372.39	-59.60	182.78	4634.27	8.35	56.22	59422.74	15506.99	573.24	410.07
Std. $(\sigma )$	9410.80	1291.23	110.11	6994.23	521.30	84.29	10286.86	5498.09	711.15	486.55
Min	172.75	-9515.50	6.00	0.00	-4604.50	0.00	33927.90	5157.40	-3451.60	-1690.00
Max	47708.00	11383.25	1940.55	32474.50	5914.50	842.83	83494.50	33384.40	5569.00	3847.00

Table 4

Aggregate error statistics for the MAE, RMSE, CRPS, ES and the CR and WS for the $50\%,90\%$ and $99\%$ prediction interval. Colour indicates performance. The best value for each scoring rule is highlighted.

3.3 Spot Auction Curves and Elasticity

The impact of forecast errors on $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ depends on the steepness of the merit-order (Kremer et al., 2021, 2020; Kulakov and Ziel, 2020). 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 (RES) supply is stronger. Under uncertainty of future 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 RES supply should be smaller and hence the volatility of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ .

There are different methods to model the merit-order used in practice and academia. Fundamental methods as developed by Pape et al. (2016); Gürtler and Paulsen (2018); Beran et al. (2019) are complex, data-intensive and rely heavily on assumptions. For this reason, Kremer et al. (2021, 2020) develop an econometric model based on He et al. (2013) 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 Balardy (2022) and Kulakov and Ziel (2020) and has three advantages compared to the approach of He et al. (2013): 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 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 (on this issue see e.g. Pape et al., 2016). 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 Kulakov and Ziel (2020, 2019) 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, theday-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 Kulakov and Ziel (2019) and Coulon et al. (2014), so here only a brief introduction is given. Figure 8a 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 (Knaut and Paulus, 2016; Coulon et al., 2014). 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. Coulon et al. (2014) and Kulakov and Ziel (2019) 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. 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. Kulakov and Ziel (2019) elaborate in detail on the econometric framework, which is adopted in this paper and the implications for the different market participants.

Figure 8

Transformed auction curves for June 1, delivery hour $s=9$ and the calculation of the merit-order slope coefficient. Note how the intersection of $\mathrm{SU}{\mathrm{P}}^{d,s}\left({\mathrm{DEM}}_{\mathrm{inelastic}}^{d,s}\right)$ again yields $P_{\mathrm{DA}}^{d,s}$ for Panels (a) and (b). Figures truncated on the $y$ -axis to $[-20,100]$ .

Figure 8 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_{\mathrm{DA}}^{d,s}$ . The notation follows largely Kulakov and Ziel (2019, 2020). Define the supply and demand curves as a mapping of volumes to prices by $\mathrm{SU}{\mathrm{P}}_{\mathrm{WS}}:(0,\infty )\to [P_{\min },P_{\max }]$ and $\mathrm{DE}{\mathrm{M}}_{\mathrm{WS}}:(0,\infty )\to [P_{\min },P_{\max }]$ . Due to strict monotonicity the inverse ${{\mathrm{SUP}}_{\mathrm{WS}}}^{-1}$ and ${{\mathrm{DEM}}_{\mathrm{WS}}}^{-1}$ always exist. Hence, ${\mathrm{SUP}}_{\mathrm{WS}}^{d.s}(q)=P$ is the supply or sell curve and ${\mathrm{DEM}}_{\mathrm{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 ${\mathrm{DEM}}_{\mathrm{inelastic}}^{d,s}={\mathrm{DEM}}_{\mathrm{WS}}^{d,s-1}\left(P_{\min }\right)$ where $P_{\min }=-500\mathrm{EUR}/$ MWh is the minimum price at the day-ahead auction (EPEX SPOT SE, 2018,see Chapter 2 of this paper). The transformed inverse supply curve can be written as:

$$\begin{array}{c}\mathrm{SU}{P}^{d,s-1}(z)=\underset{\mathrm{inverted}\mathrm{wholesale}\mathrm{supply}\mathrm{curve}}{\underbrace{{\mathrm{SUP}}_{\mathrm{WS}}^{d,s-1}(z)}}+\underset{\mathrm{flipped}\mathrm{inverted}\mathrm{wholesale}\mathrm{demand}\mathrm{curve}}{\underbrace{{\mathrm{DEM}}_{\mathrm{inelastic}}^{d,s}-{\mathrm{DEM}}_{\mathrm{WS}}^{d,s-1}(z)}}\hfill \end{array}$$

(1)

As the curves are monotonic, $\mathrm{SUP}{P}^{d,s-1}(z)$ also defines $\mathrm{SU}{\mathrm{P}}^{d,s}(q)$ . As it is clearly visible in Figure 8, the original equilibrium is reached at the point $P_{\mathrm{DA}}^{d,s}=\mathrm{SU}{P}^{d,s}\left({\mathrm{DEM}}_{\mathrm{inelastic}}^{d,s}\right)$ . For the transformed curves it now holds that for the resulting clearing price, shifting ${\mathrm{DEM}}_{\mathrm{inelastic}}^{d,s}$ by some quantity $x$ equals shifting $\mathrm{SU}{P}^{d,s-1}(z)$ by $-x$ , as the demand is perfectly inelastic.

Under the assumptions that the merit-order does not change significantly between dayahead 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 Kremer et al. (2021, 2020). The second assumption is discussed above. The last known 5-minute-interval-VWAP before the start of the simulation is $P_{\mathrm{ID},0}^{d,s}$ . Under the Market Efficiency Hypothesis (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 ${\mathrm{DEM}}_{\mathrm{implied}}^{d,s}=\mathrm{SU}{\mathrm{P}}^{d,s-1}\left(P_{\mathrm{ID},0}^{d,s}\right)$ . As already mentioned in Chapter 2, the lower and upper price limits at the day-ahead auction are [-500,3000] EUR/MWh, while in the intraday market these are set to $[-9999,9999]$ EUR/MWh. Hence, it might be possible that $P_{\mathrm{ID}.0}^{d,s}$ is outside the domain of $\mathrm{SU}{\mathrm{P}}^{d,s-1}(z)$ . This case, however, does not occur in the dataset used in this paper. In the spirit of Balardy (2022) and Kulakov and Ziel (2019), the measure for the elasticity ${\mathrm{MO}}_q^{d,s}$ is calculated as a finite central difference quotient of the transformed supply curve around ${\mathrm{DEM}}_{\mathrm{implied}}^{d,s}$ :

$$\begin{array}{c}{\mathrm{MO}}_q^{d,s}=\frac{\mathrm{SU}{\mathrm{P}}^{d,s}({\mathrm{DEM}}_{\mathrm{implied}}^{d,s}+q)-\mathrm{SU}{\mathrm{P}}^{d,s}({\mathrm{DEM}}_{\mathrm{implied}}^{d,s}-q)}{2·q}\hfill \end{array}$$

(2)

where $q=\{500,1000,2000\}$ MWh. It is defined in EUR $/\mathrm{MW}{\mathrm{h}}^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\mathrm{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, Kulakov and Ziel (2019) choose $q=100\mathrm{MWh}$ , Balardy (2022) chooses $q=500$ MWh. In this paper, slightly higher $q$ are selected to accommodate the fact that the standard deviation of $\mathrm{\Delta }{\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s}$ and $\mathrm{\Delta }{\hat{S}}_{\mathrm{DA}}^{d,s}$ ,ID is roughly between $500\mathrm{MW}$ and 1300 MW (see Table 3). 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 8 (c) shows the intuition of the slope coefficient for $q=2000$ . Figure 9 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.

Figure 9

Boxplots for the merit-order slope $M_{2000}$ grouped by the price level. Colouring indicates the Number of observations, note the log-scale of the colorbar. Upper and lower end of the box correspond to $Q_{0.25}$ and $Q_{0.75}$ respectively.

4. Electricity Price Models

Recall from the introduction that the price differences $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ follow a certain distribution that we denote as $G^{d,s}$ which is a mixture distribution

$$\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}~G^{d,s}=(1-{\alpha }_t^{d,s}){\delta }_0+{\alpha }_t^{d,s}{F}^{d,s}$$

with the Dirac distribution ${\delta }_0$ , the continuous distribution $F^{d,s}$ and the Bernoulli variable ${\alpha }_t^{d,s}$ with probability ${\pi }_t^{d,s}$ . The two-stage approach is:

1. Model the probability ${\pi }_t^{d,s}$ for a trade-event ${\alpha }_t^{d,s}=1$ by a logistic regression model.

2. Model the distribution $F_t^{d,s}$ by skewed Student’s $t$ -distribution respectively the Johnson’s $S_U$ and estimate ${\mathit{\theta }}^{d,s}$ using the GAMLSS-framework (Stasinopoulos and Rigby, 2007; Stasinopoulos et al., 2017; Stasinopoulos and Rigby, 2005; Stasinopoulos et al., 2018).

This introduces a dependence structure between the parameters of ${\delta }_0$ and $F$ , as the probability ${\pi }^{d,s}$ is explained by past realizations of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ and ${\alpha }_t^{d,s}$ . In the following three sections, the logistic model, the GAMLSS model and the used benchmarks are presented.

4.1 Logistic regression model for $\alpha$

The binary variable ${\alpha }_t^{d,s}$ will be modelled by a regularized logistic regression model (Tibshirani, 1996; Meier et al., 2008) in the implementation of Friedman et al. (2010) using coordinate descent. Generally, for a logistic model

$$\begin{array}{c}\log \left(\frac{\pi }{1-\pi }\right)={\mathbf{X}}^T\mathit{\beta }\hfill \end{array}$$

(3)

for the Bernoulli variable $\alpha$ with probability $P(\alpha =1)=\pi$ , the LASSO estimator ${\hat{\beta }}^{\mathrm{LASSO}}$ is given by

$$\begin{array}{c}{\hat{\mathit{\beta }}}^{\mathrm{LASSO}}=\mathrm{argmin}(-l(\mathit{\beta },\tilde{\mathbf{x}})+\lambda \parallel \mathit{\beta }{\parallel }_1),\hfill \end{array}$$

(4)

where $\lambda$ is a tunable shrinkage parameter. The corresponding log-likelihood $l$ is given by

$$\begin{array}{c}l(\mathit{\beta },\tilde{\mathbf{x}})=\frac{1}{N}\sum _{i=1}^N{\alpha }_i{\tilde{\mathbf{X}}}_i^T\mathit{\beta }-\log (1+\exp \left({\tilde{\mathbf{X}}}_i^T\mathit{\beta }\right)),\hfill \end{array}$$

(5)

where $\tilde{\mathbf{X}}$ is a standardization of $\mathbf{X}$ . The parameter $\lambda$ is optimized from an exponential grid of 100 values by choosing the minimum Bayesian Information Criterion (BIC), i.e. ${\lambda }^{\mathrm{opt}}=\mathrm{argminBIC}(\lambda )$ using the glmnet package by Friedman et al. (2010).

Here, the logit function for ${\alpha }_t^{d,s}$ is explained by four components: the impact of past price differences, the time to maturity and weekday effects, fundamental variables such as RES forecasts, outages and the slope of the merit-order, and a regression on averaged past ${\alpha }_t^{d,s}$ . Intuitively, the probability of trades should rise with higher $\mathrm{\Delta }{P}_{\mathrm{ID},t,}^{d,s}$ , closer to delivery, with increasing wind and solar forecasts and with increased recent trading activity measured by past ${\alpha }_t^{d,s}$ , but decrease on the weekends and the transition day Monday.

$$\begin{array}{c}\multicolumn{1}{c}{\log \left(\frac{{\pi }_t^{d,s}}{1-{\pi }_t^{d,s}}\right)=\underset{\mathrm{Pricedifferences}}{\underbrace{{\beta }_0+\sum _{j=1}^3{\beta }_j\mathrm{\Delta }{P}_{\mathrm{ID},t-j}^{d,s}+\sum _{j=1}^6{\beta }_{3+j}\left|\mathrm{\Delta }{P}_{\mathrm{ID},t-j}^{d,s}\right|+{\beta }_{10}\sum _{j=7}^{12}\left|\mathrm{\Delta }{P}_{\mathrm{ID},t-j}^{d,s}\right|}}}\\ \multicolumn{1}{c}{\underset{\mathrm{Timedummies}}{\underbrace{+{\beta }_{11}\mathrm{MON}(d)+{\beta }_{12}\mathrm{SAT}(d)+{\beta }_{13}\mathrm{SUN}(d)+\sum _{j=1}^{31}{\beta }_{13+j}\mathrm{TTD}(t)}}}\\ \multicolumn{1}{c}{\underset{Day-aheadfundamental\mathrm{var}iables}{\underbrace{+{\beta }_{47}{\widehat{D}}_{DA}^{d,s}+{\beta }_{48}{\widehat{W}}_{DA}^{d,s}+{\beta }_{49}{\widehat{S}}_{DA}^{d,s}+{\beta }_{50}{O}_{DA}^{d,s}}}}\\ \multicolumn{1}{c}{\underset{\mathrm{Day}-\mathrm{aheadtointradayforecastupdates}}{\underbrace{+{\beta }_{51}\mathrm{\Delta }\widehat{W}{\mathrm{DA},\mathrm{ID}}^{d,s,+}+{\beta }_{52}\mathrm{\Delta }\widehat{W}{\mathrm{DA},\mathrm{ID}}^{d,s,-}+{\beta }_{53}\mathrm{\Delta }\widehat{S}{\mathrm{DA},\mathrm{ID}}^{d,s,}+{\beta }_{54}\mathrm{\Delta }\widehat{S}{\mathrm{DA}}^{d,s,-}}}}\\ \multicolumn{1}{c}{\underset{\mathrm{Standarddeviationofforecastupdates}}{\underbrace{+{\beta }_{55}{\sigma }_{\mathrm{DA}}^{d,s}(\mathrm{\Delta }\widehat{W})+{\beta }_{56}{\sigma }_{\mathrm{DA},\mathrm{ID}}^{d,s}(\mathrm{\Delta }\widehat{S})}}\underset{\mathrm{Intradaychangesinplannedandunplannedoutages}}{\underbrace{+{\beta }_{57}\mathrm{\Delta }{O}_{\mathrm{DA},\mathrm{ID}}^{d,s,\mathrm{planed}}+{\beta }_{58}\mathrm{\Delta }{O}_{\mathrm{DA},\mathrm{ID}}^{d,s,\mathrm{unned}}}}}\\ \multicolumn{1}{c}{\underset{\mathrm{Day}-\mathrm{aheadtot}-1\mathrm{pricespread}}{\underbrace{+{\beta }_{59}|P_{\mathrm{DA}}^{d,s}-P_{\mathrm{ID},\mathrm{t}-\mathrm{s}}^{d,s}|}}+\underset{\mathrm{Slopeofthemerit}-\mathrm{order}}{\underbrace{\sum _{j=1}^3{\beta }_{59+j}{\mathrm{MO}}_j^{d,s}}}+\underset{\mathrm{Regressionon}{\overline{\alpha }}_t^{d,s}}{\underbrace{\sum _{j=1}^{12}{\beta }_{62+j}{\overline{\alpha }}_{t-j}^{d,s}}},}\end{array}$$

(6)

where ${\overline{\alpha }}_{t-j}^{d,s}=1/j·{\sum }_{i=1}^j{\alpha }_{t-i}^{d,s}$ , the average of the last $j$ observed values of ${\alpha }_t^{d,s}$ . 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 Narajewski and Ziel (2020b). $\mathrm{SAT}(d),\mathrm{SUN}(d)$ , and MON $(d)$ are dummies for the weekday of $d\mathrm{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 LASSO of Tibshirani (1996). Note however, that for the one year training set used in this paper, we have $365·T=365·31=11315$ observations and are still in a setting where $n>>p$ . The number of observations is sufficiently larger than the number of parameters, hence identification is not an issue here.

4.2 GAMLSS Framework

This chapter briefly introduces the GAMLSS-framework used to model

$$\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}\mid {\alpha }_t^{d,s}=1~F.$$

The GAMLSS is an extension of the Generalized Additive Models (GAM) introduced by Hastie and Tibshirani (1987, 1990). It allows to model not only the expected value of the variable $Y~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 Stasinopoulos and Rigby (2005, 2007); Stasinopoulos et al. (2018) and the manual of the R-package GAMLSS (Stasinopoulos et al., 2017). The notation in the following paragraphs follows the notation of aforementioned sources. 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 $\mathbf{Y}=(Y_1,Y_2,\dots ,Y_n)$ be a vector of $i=1,\dots ,n$ independent observations $Y_i$ . The GAMLSS-framework assumes that $Y_i$ have the probability density function

$$f(y_i\mid {\mu }_i,{\sigma }_i,v_i,{\tau }_i)$$

where each of the distribution parameters can be a smooth function of the explanatory variables. We denote as ${\mathit{\theta }}_i=({\theta }_{i,1}{\theta }_{i,2},{\theta }_{i,3},{\theta }_{i,4})=({\mu }_i,{\sigma }_i,v_i,{\tau }_i)$ the vector of $k=1,\dots ,4$ distribution parameters which are usually known as the location, scale and shape parameters ${\theta }_{i,k}$ . For the distributions used in this paper, $v_i$ denotes the skewness and ${\tau }_i$ denotes the kurtosis. $\mathit{\theta }$ is a matrix whose individual components have the indices $i$ and $k$ . The vectors ${\mathit{\theta }}_i$ and ${\mathit{\theta }}_k$ are defined along the 2 axis of $\mathit{\theta }$ . Formally, we have

$$Y_i~F({\mu }_i,{\sigma }_i,v_i,{\tau }_i)\iff Y_i~F{(}_i)$$

For each $k$ , let $g_k(·)$ be a known and monotonic link function that relates the distribution parameters ${\mathit{\theta }}_k$ to the predictor ${\mathit{\eta }}_k$ . We consider the GAMLSS model equation

$$\begin{array}{c}{g}_k\left({\mathit{\theta }}_k\right)={\mathit{\eta }}_k={\mathbf{X}}_k{\mathit{\beta }}_k\hfill \end{array}$$

(7)

where ${\mathit{X}}_k$ is a $n\times J_k$ fixed design matrix and ${\mathit{\beta }}_k^{\prime }=({\beta }_{1,k},{\beta }_{2,k},\dots ,{\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:

$$g_{\mathrm{ident}}(z)=z$$

(8)

$$g_{\log }(z)=\log (z)$$

(9)

$$g_{\mathrm{logident}}(z)=\log (z)1(z\le 1)+(z-1)1(z>1)$$

(10)

$$g_{\mathrm{logshift}2}(z)=\log (z-2)$$

(11)

We use $g_{\mathrm{ident}}$ for the location parameters $\mu$ in both distribution assumptions and for the skewness parameter $v$ of the skewed t-distribution and Johnson’s $S_U$ distribution, which both have support $(-\infty ,\infty )$ for $v$ . For both distributions, we used the link function $g_{\mathrm{logident}}$ for the scale parameter $\sigma$ to avoid the exponential inverse for large estimates, thus improving the robustness of the estimation estimation (Ziel, 2021; Narajewski and Ziel, 2020b). The link function $g_{\mathrm{logshift}}$ is simply the natural logarithm shifted to 2 to preserve the condition $\tau >2$ for the scale of the skewed $t$ -distribution. For the remaining kurtosis parameter $\tau$ for Johnson’s $S_U$ we consider $g_{\log }$ .

Ziel and Muniain (2021) extend the GAMLSS framework to allow for regularized LASSO estimation. As with the logistic model, we employ the BIC to select the optimal shrinkage parameter $\lambda$ . The adaptive LASSO estimator $\mathit{\beta }{*}_k$ is used. It is defined as

$$\begin{array}{c}\mathit{\beta }{*}_k=\mathrm{argmi}{\mathrm{n}}_{\beta }{|y-\sum _{j=1}^J{x}_j{\beta }_{j,k}|}^2+{\lambda }_n\sum _{j=1}^J{\hat{w}}_{j,k}\left|{\beta }_{j,k}\right|\hfill \end{array}$$

(12)

with the weights vector ${\hat{\mathit{w}}}_k=1/|\hat{\beta }\mid \gamma$ . $\hat{\beta }$ denotes a root- $n$ consistent estimator such as ordinary least squares (Zou, 2006).

Let us now relate the abstract notation $Y_i$ to $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ . The distribution $F\left({\mathit{\theta }}_t^{d,s}\right)$ is fitted to all $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}\mid {\alpha }_t^{d,s}=1$ . The abstract index $i=1,\dots ,N$ is replaced by the combination of the superscript index $d=1,\dots ,D$ and the subscript index $t=1,\dots ,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{\mathit{\theta }}}_t^{d,s}=({\hat{\mu }}_t^{d,s},{\hat{\sigma }}_t^{d,s},{\hat{\tau }}_t^{d,s}{\hat{v}}_t^{d,s})$ and accordingly parameter estimates ${\mathit{\beta }}_k^{d,s}$ that condition ${\hat{\mathit{\theta }}}_t^{d,s}$ on our explanatory variables. Analogously, we can also define the vector ${\hat{\mathit{\theta }}}_k^{d,s}=({\theta }_{1,k}^{d,s},\dots ,{\theta }_{T,k}^{d,s})$ along the time-axis $t$ . We explain all moments of the distribution by the same set of explanatory variables. Narajewski and Ziel (2020b) 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 A.2.

For each distribution parameter $k$ , the model for ${\hat{\theta }}_{t,k}^{d,s}$ reads:

$$\begin{array}{c}\multicolumn{1}{c}{g_k\left({\theta }_{t,k}^{d,s}\right)={\beta }_{k1(k\ge 2),0}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Pricedifferences}}{\underbrace{\sum _{j=1}^3{\beta }_{k,j}\mathrm{\Delta }{P}_{\mathrm{ID},t-j}^{d,s}}}+\underset{\mathrm{Absolutepricedifferences}}{\underbrace{\sum _{j=1}^6{\beta }_{k,3+j}\left|\mathrm{\Delta }{P}_{\mathrm{ID},t-j}^{d,s}\right|+{\beta }_{10}\sum _{j=7}^{12}\left|\mathrm{\Delta }{P}_{\mathrm{ID},t-j}^{d,s}\right|}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Timedummies}}{\underbrace{{\beta }_{k,11}\mathrm{MON}(d)+{\beta }_{k,12}\mathrm{SAT}(d)+{\beta }_{k,13}\mathrm{SUN}(d)}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Day}-\mathrm{aheadfundamentalvariables}}{\underbrace{{\beta }_{k,14}{\hat{L}}_{\mathrm{DA}}^{d,s}+{\beta }_{k,15}{\hat{W}}_{\mathrm{DA}}^{d,s}+{\beta }_{k,17}{\hat{S}}_{\mathrm{DA}}^{d,s}+{\beta }_{k,18}{O}_{\mathrm{DA}}^{d,s}}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Day}-\mathrm{ahead}\mathrm{to}\mathrm{intraday}\mathrm{forecast}\mathrm{updates}}{\underbrace{{\beta }_{k,19}\Delta {\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s,+}+{\beta }_{k,20}\Delta {\hat{W}}_{\mathrm{DA},\mathrm{ID}}^{d,s,-}+{\beta }_{k,21}\Delta {\hat{S}}_{\mathrm{DA},\mathrm{ID}}^{d,s,+}+{\beta }_{k,22}\Delta {\hat{S}}_{\mathrm{DA},\mathrm{ID}}^{d,s,-}}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Standard}\mathrm{deviation}\mathrm{of}\mathrm{forecast}\mathrm{updates}}{\underbrace{{\beta }_{k,23}{\sigma }_{\mathrm{DA},\mathrm{ID}}^{d,s}\left(\Delta \hat{W}\right)+{\beta }_{k,24}{\sigma }_{\mathrm{DA},\mathrm{ID}}^{d,s}\left(\Delta \hat{S}\right)}}+\underset{\mathrm{Intraday}\mathrm{changes}\mathrm{in}\mathrm{planned}\mathrm{and}\mathrm{unplanned}\mathrm{outages}}{\underbrace{{\beta }_{k,25}\Delta O_{\mathrm{DA},\mathrm{ID}}^{d,s,\mathrm{planned}}+{\beta }_{k,26}\Delta O_{\mathrm{DA},\mathrm{ID}}^{d,s,\mathrm{unplanned}}}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Lagged}{\alpha }_{\mathrm{t}}^{\mathrm{d},\mathrm{s}}}{\underbrace{{\beta }_{k,27}{\alpha }_{t-1}^{d,s}+{\beta }_{k,28}{\alpha }_{t-2}^{d,s}}}+\underset{\mathrm{Day}-\mathrm{ahead}\mathrm{to}\mathrm{t}-1\mathrm{price}\mathrm{spread}}{\underbrace{{\beta }_{k,29}|P_{\mathrm{DA}}^{d,s}-P_{\mathrm{ID},\mathrm{t}-1}^{d,s}|}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{Slope}\mathrm{of}\mathrm{the}\mathrm{merit}-\mathrm{order}}{\underbrace{\sum _{j=1}^3{\beta }_{k,29+j}{\mathrm{MO}}_i^{d,s}}}+\underset{\mathrm{Time}\mathrm{to}\mathrm{delivery}\mathrm{and}\mathrm{SIDC}\mathrm{Closing}}{\underbrace{{\beta }_{k,32}f\mathrm{TTD}\left(t\right)+{\beta }_{k,33}\mathrm{SIDC}(d,t)}}}\end{array}$$

(13)

for $k=1,2,3,4$ . The intercept is only included for $k\ge 2$ as the price differences are assumed to be centred around 0, as indicated by the summary statistics in Table 1 . The volatility is expected to rise with higher absolute past price differences, on weekends, and with higher RES generation. The strong changes between day-ahead and intraday RES and demand forecasts should also imply higher volatility. We expect the volatility to decrease with more recent trading activity measured bylagged ${\alpha }_t^{d,s}$ . $\mathrm{SIDC}(d,t)$ is a dummy variable taking the value 1 for $d\ge$ June 18,2018 and $26\le t\le 31$ , indicating that the cross-country order books are closed. $f_{\mathrm{TTD}}(t)$ models the non-linear impact of the time to delivery and takes the form $f_{\mathrm{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 3.3, 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.

4.3 Benchmark Models

Lastly, some simple benchmark models are introduced. Even though the main focus of this paper is on modelling the volatility and its influencing factors, simple benchmarks can serve as a valuable comparison to identify potential areas for model improvement. As our study shares the conceptual set-up with Narajewski and Ziel (2020b), it is natural to employ similar benchmark models. Additionally, we compare our approach to classical time series methods such as AutoRegressive Integrated Moving Average (ARIMA) models. The following section introduces these models in more detail.

Narajewski and Ziel (2020b) 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.

• The Naive benchmark randomly draws past trajectories.

• For MV.N and MV.t, a multivariate normal respectively $t$ -distribution is fitted to the vector of price differences. Forecasts are randomly drawn from the distribution.

• The RW.N, RW.t and RW.t.mix.D are random-walk type of models, where the distribution parameters are estimated from the in-sample data. The RW.t.mix.D also includes a simple mixture term for ${\alpha }_t^{d,s}$ by estimating ${\pi }_t^{d,s}$ as the empirical mean of ${\alpha }_t^{d,s}$ .

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 2, 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.

The classic $\mathrm{ARIMA}(p,k,q)$ model is defined as follows:

$$\begin{array}{c}(1-\sum _{i=1}^p{q}_i{L}^i)(1-L)\Delta P_{\mathrm{ID},t}^{d,s}=\delta +(1+\sum _{i=1}^q{\theta }_i{L}^i){\epsilon }_t\hfill \end{array}$$

(14)

is an $\mathrm{ARIMA}(p,k,q)$ process with drift $\delta /(1-\sum ^{{\phi }_i}).L$ denotes the lag operator. A full treatment of ARIMA models can be found in e.g. Shumway and Stoffer (2017). We estimate the ARIMA( $p,k,q)$ models using the auto.arima() function in the forecast package (Hyndman et al., 2020). The function uses a stepwise approach to fit the lag order for $p$ and $q$ based on the BIC. We fix the integrationorder $k=0$ according to our initial analysis. ⁴ 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.

5. Forecasting Study and Evaluation

5.1 Study Design and Simulation Algorithm

We employ the well-known rolling window forecasting study design, which is common in energy price forecasting (see e.g. Ziel et al., 2015; Gianfreda and Bunn, 2018; Nowotarski and Weron, 2018; Janke and Steinke, 2019; Uniejewski et al., 2019; Narajewski and Ziel, 2020b). This setting reduces the impact of structural breaks within the data and ensures a robust setting for the comparison of predictive performance using the DM-test (Diebold and Mariano, 2002; Diebold, 2015). The scheme is visualized in Figure 10. 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.

Figure 10

Structure of the rolling window forecasting study. Blue denotes in-sample data, red denotes the out-of-sample forecast.

Forecasts are issued for all delivery hours $s=0,\dots ,23$ . For each delivery hour, the forecast consists of $j=1,\dots ,M$ paths with $M=1000$ paths of $t=1,\dots ,31$ steps. Generally, let variables with superscript $[j]$ denote simulated values on path $j$ and hence $P_{\mathrm{ID},t}^{d,s,j]}$ denotes the simulation for step $t$ for delivery period $d,s$ in the path $j$ . The vector notation ${\mathbf{P}}_{\mathrm{ID}}^{d,s,[j]}=(P_{\mathrm{ID},1}^{d,s,[]]},\dots ,P_{\mathrm{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 (Asmussen and Glynn, 2007; Narajewski and Ziel, 2020b). 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 ${\alpha }_t^{\mathrm{ds},[j]}$ is simulated $M$ times from the Bernoulli distribution $B\left({\hat{\pi }}_f^{d,s,[j]}\right)$ and the price difference $\mathrm{\Delta }{P}_{\mathrm{ID},t,[j]}^{d,t}$ is sampled $M$ times from the distribution $F\left({\hat{\theta }}_t^{d,s,[j]}\right)=F({\hat{\mu }}_{,{\sigma }_t}^{d,s,[j]}{\hat{v}}_t^{d,s,[j]}$ and ${\hat{\tau }}_t^{d,s,[j]})$ . The price $P_{\mathrm{ID},t}^{d,s,[j]}$ is then calculated as:

$$\begin{array}{c}{P}_{\mathrm{ID},t}^{d,s,j]}=P_{\mathrm{ID},t-1}^{d,s,[j]}\aleph P_{\mathrm{ID},t}^{d,s,[j]}{\alpha }_t^{d,s,[j]}\hfill \end{array}$$

(15)

The algorithm is visualized in Figure 11. The estimates for ${\hat{\pi }}_t^{d,s,[j]},{\hat{\mu }}_t^{d,s,[j]},{\hat{\sigma }}_t^{d,s,[j]},{\hat{v}}_t^{d,s,[j]}$ and ${\hat{\tau }}_t^{d,s,[j]}$ for the first step $t=1$ are all equal, but begin to differ from $t\ge 2$ onwards as the paths develop individually. Therefore, the prediction matrix needs to be updated dynamically for each path and each step.

Figure 11

Simulation Algorithm. After $t=3$ , the steps until $t=T$ are omitted. For each path $j$ , an individual regression matrix is created with the information of the path’s past development. Together, the results of the right, green column yield the path vector

5.2 Forecast Evaluation

The mean and median trajectory are evaluated using the Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) respectively. For the probabilistic evaluation of the generated scenarios the Energy Score (ES), Continuously Ranked Probability Score (CRPS) and the empirical coverage ratio are used. Additionally, the Winkler-Score (WS) is used to evaluate the coverage of an $(1-\alpha )·100\mathrm{\%}$ -Prediction Interval (PI). The ES, CRPS and the WS are strictly proper scoring rules (Nowotarski and Weron, 2018; Gneiting and Raftery, 2007; Ziel and Berk, 2019). To draw conclusions about the statistical significance of the difference in forecasting performance for each model, the DM-test is used. All measures are widely employed in academia and practice.

Formally, the RMSE and MAE are defined as:

$$\begin{array}{c}\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{NST}}\sum _{d=1s=1}^N\sum _{t=1}^S\sum _{\mathrm{ID},t}^T{(P^{d,s}-\frac{1}{M}\sum _{j=1}^M{P}_{\mathrm{ID},t}^{d,s,[j]})}^2}\hfill \end{array}$$

(16)

and:

$$\begin{array}{c}\mathrm{MAE}=\frac{1}{\mathrm{NST}}\sum _{d=1s=1}^N\sum _{t=1}^S\sum _{\mathrm{ID},t}^T|P_{\mathrm{ID}}^{d,s}-\mathrm{med}\left(P_{\mathrm{ID},t}^{d,[j]}\right)|\hfill \end{array}$$

(17)

For an $(1-\alpha )·100\mathrm{\%}$ -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 CR is defined as:

$$\begin{array}{c}\mathrm{CR}=\frac{1}{\mathrm{NST}}\sum _{d=1s=1t=1}^N\sum ^S\sum _t^T\{\begin{array}{c}1{\mathrm{forP}}_{\mathrm{ID},t}^{d,s}\in [{\hat{L}}_t^{d,s},{\hat{U}}_t^{d,s}]\hfill \\ \multicolumn{1}{c}{0\mathrm{else}.}\end{array}\hfill \end{array}$$

(18)

The ${\mathrm{WS}}_t^{d,s}$ is defined as:

$$\begin{array}{c}{\mathrm{WS}}_t^{d,s}=\{\begin{array}{c}{\delta }_t,{\mathrm{forP}}_{\mathrm{ID},t}^{d,s}\in [{\hat{L}}_t^{d,s},{\hat{U}}_t^{d,s}]\hfill \\ \multicolumn{1}{c}{{\delta }_t+\frac{2}{\alpha }({\hat{L}}_t^{d,s}-P_{\mathrm{ID},t}^{d,s}),{\mathrm{forP}}_{\mathrm{ID},t}^{d,s}<{\hat{L}}_t^{d,s}}\\ \multicolumn{1}{c}{{\delta }_t+\frac{2}{\alpha }(P_{\mathrm{ID},t}^{d,s}-{\hat{U}}_t^{d,s}),{\mathrm{forP}}_{\mathrm{ID},t}^{d,s}>{\hat{U}}_t^{d,s}}\end{array}\hfill \end{array}$$

(19)

and aggregated as:

$$\begin{array}{c}\mathrm{WS}=\frac{1}{\mathrm{NST}}\sum _{d=1}^N\sum _{s=1}^S\sum _{t=1}^T{\mathrm{WS}}_t^{d,s}\hfill \end{array}$$

(20)

For both, CR and WS, the upper and lower bounds of the $(1-\alpha )·100\mathrm{\%}\mathrm{PI}$ are defined by the respective quantiles ${\hat{L}}_t^{d,s}=Q_{j=1,\dots ,M}^{\alpha /2}\left(P_{\mathrm{ID},t}^{d,s,[j]}\right)$ and ${\hat{U}}_t^{d,s}=Q_{j=1,\dots ,M}^{1-\alpha /2}\left(P_{\mathrm{ID},t}^{d,s,[j]}\right)$ , where $Q_{j=1,\dots ,M}^{\tau }\left(P_{\mathrm{ID},t}^{d,s,j]}\right)$ denotes the $\tau$ -th quantile of $M$ simulated $P_{\mathrm{ID},t,t}^{P_{s,s}[j]}$ prices. Comparing both, WS and CR, one can see how the WS penalizes for an observation outside the interval and rewards the forecaster at the same time for a more narrow PI. Contrary to the CR, the WS is a strictly proper evaluation measure (Nowotarski and Weron, 2018).

The CRPS (see e.g. Gneiting and Raftery, 2007;Nowotarski andWeron, 2018) is approximated by thePinball-Score (PB)

$$\begin{array}{c}{\mathrm{CRPS}}_t^{d,s}=\frac{1}{R}\sum _{\tau \in \mathcal{T}}{\mathrm{PB}}_{t,\tau }^{d,s}\hfill \end{array}$$

(21)

for a dense equidistant grid of probabilities $\mathcal{T}=\{0.01,\dots 0.99\}$ of size $R=99$ . ${\mathrm{PB}}_{t,\tau }^{d,s}$ denotes the pinball loss for probability $\tau$ . The formula is given by:

$${\mathrm{PB}}_{t,\tau }^{d,s}=\{\begin{array}{c}(1-\tau )·(Q_{j=1,\dots ,M}^{\tau }\left(P_{\mathrm{ID},t}^{d,s,\left[j\right]}\right)-P_{t,j}^{d,s})\mathrm{for}{P}_{\mathrm{ID},t}^{d,s}\le Q_{j=1,\dots ,M}^{\tau }\left(P_{\mathrm{ID},t}^{d,s,\left[j\right]}\right)\hfill \\ \multicolumn{1}{c}{\tau ·(P_{\mathrm{ID},t}^{d,s}-Q_{j=1,\dots ,M}^{\tau }\left(P_{\mathrm{ID},t}^{d,s,\left[j\right]}\right))\mathrm{else}.}\end{array}$$

(22)

The overall CRPS is calculated by the average:

$$\begin{array}{c}\mathrm{CRPS}=\frac{1}{\mathrm{NST}}\sum _{d=1}^N\sum _{s=1}^S\sum _{t=1}^T{\mathrm{CRPS}}_t^{d,s}\hfill \end{array}$$

(23)

The pinball loss is also used to evaluate the performance of different models in specific quantile levels. For this reason, the PB is aggregated as follows:

$$\begin{array}{c}\mathrm{P}{\mathrm{B}}_{\tau }=\frac{1}{\mathrm{NST}}\sum _{d=1s=1}^N\sum _{t=1}^S\sum _{t,\tau }^T{\mathrm{PB}}_{t,s}^{d,s}\hfill \end{array}$$

(24)

To measure the quality of the generated paths, Narajewski and Ziel (2020b) propose the ES. It is a generalization of the CRPS for two dimensions. Thereby, not only the approximation of the marginal distribution is evaluated, but also the generated multivariate distribution (Gneiting and Raftery, 2007; Ziel and Berk, 2019):

$$\begin{array}{c}\mathrm{E}{\mathrm{S}}^{d,s}=\frac{1}{M}\sum _{j=1}^M{\Vert {\mathbf{P}}_{\mathrm{ID}}^{d,s}-{\mathbf{P}}_{\mathrm{ID}}^{d,s,[j]}\Vert }_2-\frac{1}{·M·(M-1)}\sum _{j=1}^M\sum _{i=j+1}^M{\Vert {\mathbf{P}}_{\mathrm{ID}}^{\mathrm{ds},[j]}-{\mathbf{P}}_{\mathrm{ID}}^{d,s,[i]}\Vert }_2\hfill \end{array}$$

(25)

The average yields the overall energy score for each model:

$$\begin{array}{c}\mathrm{ES}=\frac{1}{\mathrm{NS}}\sum _{d=1s=1}^N\sum _s^S\mathrm{E}{\mathrm{S}}^{d,s}\hfill \end{array}$$

(26)

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 DM-test (Diebold and Mariano, 2002; Diebold, 2015) is routinely employed in the field of energy price forecasting (Nowotarski and Weron, 2018; Ziel and Weron, 2018; Janke and Steinke, 2019). It originally stems from the field of point forecasting, however Diebold (2015) notes that the test is agnostic to the scoring rule used to evaluate forecasts. Hence, using strictly proper probabilistic scoring rules, such as the CRPS and ES loss, the DM test can be applied to probabilistic forecasts as well (see e.g. Diebold, 2015; Nowotarski and Weron, 2018). Following Narajewski and Ziel (2020b) and Ziel and Weron (2018), the DM-test is employed in a multivariate fashion. Hence, let ${\mathbf{L}}_A^d=(L_A^{d,1},\dots ,L_A^{d,S})$ and ${\mathbf{L}}_A^d=(L_B^{d,1},\dots ,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 each model $A$ and $B$ , the $N\times S$ vector of losses is 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.

$$\begin{array}{c}\Delta {\mathit{L}}_{A,B}^d=\Vert {\mathit{L}}_A^d{\Vert }_1-\Vert {\mathit{L}}_B^d{\Vert }_1.\hfill \end{array}$$

(27)

For example, for the ES and the Naive model, the loss vector ${\mathit{L}}_{\mathrm{Naive}}^d=({\mathrm{ES}}_{\mathrm{Naive}}^{d,1},\dots ,{\mathrm{ES}}_{\mathrm{Naive}}^{d,S})$ .

We test the loss differential series for stationarity using the augmented Dickey-Fuller (ADF) test (Dickey and Fuller, 1979, 1981) and reject the $H_0$ of an unit root at the $5\mathrm{\%}$ significance level for all loss differential series. Harvey et al. (1997) propose the usage of the $t$ -distribution with $v=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_{\mathrm{DM}}^{\mathrm{HLN}},h=1=\sqrt{\frac{N+3}{N}}·t_{\mathrm{DM}}~t(0,1,N-1)$ under the $H_0$ , where $N$ is the length of the loss differential series $\mathrm{\Delta }{L}_{A,B}^d$ and $h$ denotes the forecast horizon. The standard deviation is computed using an autocorrelation-consistent estimator. For each model pair, two one-sided 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 (Hyndman and Khandakar, 2008; Hyndman et al., 2020).

6. Results

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.

6.1 Out-of-sample Analysis: Forecasting Performance on Test Data

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 DM-test. 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 $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 21 in Appendix B). 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 Naive and Mix.JSU . The Auto.ARIMA model yields too low coverage ratio for the $90\mathrm{\%}$ and $95\mathrm{\%}$ PIs, which on the other side result in comparably very good Winkler scores for these PIs.

Figure 12 shows the ES relative to the delivery hours $s$ . Relative to the Naive, the GAMLSS-based models show an improved forecasting performance in the peak hours (Plot $12\mathrm{b}$ ). The error of the RW.N and Mix.SST models explodes for hour $s=6$ .

Figure 12

Plot (a) shows the ES and (b) its ratio to Naive over the delivery hours $s$ .

The PB over $\mathcal{T}=0.01,\dots ,0.99$ is shown in Figure 13. 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 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 driven by its sensitivity to outliers. The development of the CRPS throughout the simulation window is shown in Figure 14. 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. The relative error of most models towards the Naive decreases throughout the simulation window, however, it increases for the Auto.ARIMA. This might indicate that learning the model parameters of past trading sessions can be beneficial compared to learning the parameters

Figure 13

Plot (a) shows the $\mathrm{PB}$ and (b) its ratio to Naive over the quantile range $\mathcal{T}$ .

Figure 14

Plot (a) shows the CRPS and (b) its ratio to Naive over the time to delivery.

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 15 shows the $p$ -values for the pairwise DM-tests for the ES and CRPS. The lower the $p$ -value, the more significant is the difference in the forecasting performance, which implies that the modelon the column (or $x$ -axis) outputs superior forecasts than the model on the row (or $y$ -axis). Generally, the $p$ -values for the CRPS and 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. 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 12 to 14 .

Figure 15

Pairwise $p$ -values for the DM-test for the CRPS and ES loss. The closer the $p$ -value to 0, the more significant is the difference between the model on the column (better) and the model on the row (worse).

6.2 In-sample Analysis: Estimated Coefficients and their Development

Given the strong probabilistic forecasting performance of the Mix.JSU we turn to an in-sample analysis of the estimated coefficients. Compared to black-box deep learning algorithms, the parametric 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 5 to 8 present the estimated scaled coefficients for $d=$ January 22nd, 2017, the first out-of-sample day. Scaled coefficients correspond to mean-variance scaled inputs. Hence, the coefficients are 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.

Table 5

Estimated scaled coefficients for $\mu$ (expected value) on the first day of the test set.

Table 6

Estimated scaled coefficients for $\sigma$ (scale / volatility) on the first day of the test set.

Table 7

Estimated scaled coefficients for $v$ (skewness) on the first day of the test set.

Table 8

Estimated scaled coefficients for $\tau$ (kurtosis) on the first day of the test set.

For $\mu$ , only the first lag of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 consistent with the findings of Narajewski and Ziel (2020a), who find the most recent price to be among the most important features for forecasting the $\mathrm{I}{\mathrm{D}}_3$ as well as with the results of Kremer et al. (2021, 2020), 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 inthe price already. These results support the notion of weak-form market efficiency already indicated by Narajewski and Ziel (2020a) and Kuppelwieser and Wozabal (2021).

For the volatility $\sigma$ , we present coefficients in similar fashion in Table 6. 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 ${\alpha }_t^{d,s}$ 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 $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ to the volatility parameter of the distribution of $\mathrm{\Delta }{P}_{\mathrm{ID}.t}^{d,s}$ thus generalizes the results of Kremer et al. (2021, 2020). Contrary to Baule and Naumann (2021), 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 6 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 16. 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 Narajewski and Ziel (2020b) and contradict Kath (2019), who finds no evidence of rising volatility due to SIDC.

Figure 16

Estimated coefficients for the SIDC dummy across the test set.

For the skewness parameter $v$ 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 8 giving the estimated coefficients for the kurtosis parameter $\tau$ . We find a negative impact of lagged ${\alpha }_t^{d,s}$ and a positive impact of lagged $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ . Thus, we expect the distribution of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 ${\alpha }_t^{d,s}$ and $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ 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 Kremer et al. (2020)’s finding that the behaviour of night contracts is more driven by trading-related variables than by fundamentals. For $\tau$ , we find a negative impact of the merit-order slope parameter. This implies that a steep merit-order leads to heavier tails for the distribution of $\mathrm{\Delta }{P}_{\mathrm{ID},t}^{d,s}$ . 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.

7. Discussion and Conclusion

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 Narajewski and Ziel (2020b) 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 GAMLSS-LASSO for all distribution moments.

Our results are two-fold: 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 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 Serafin et al. (2022) 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 jumps 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 analyzed 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.