Multivariate Simulation-based Forecasting for Intraday Power Markets: Modelling Cross-Product Price Effects

Intraday electricity markets play an increasingly important role in balancing the intermittent generation of renewable energy resources, which creates a need for accurate probabilistic price forecasts. However, research to date has focused on univariate approaches, while in many European intraday electricity markets all delivery periods are traded in parallel. Thus, the dependency structure between different traded products and the corresponding cross-product effects cannot be ignored. We aim to fill this gap in the literature by using copulas to model the high-dimensional intraday price return vector. We model the marginal distribution as a zero-inflated Johnson's $S_U$ distribution with location, scale and shape parameters that depend on market and fundamental data. The dependence structure is modelled using latent beta regression to account for the particular market structure of the intraday electricity market, such as overlapping but independent trading sessions for different delivery days. We allow the dependence parameter to be time-varying. We validate our approach in a simulation study for the German intraday electricity market and find that modelling the dependence structure improves the forecasting performance. Additionally, we shed light on the impact of the single intraday coupling (SIDC) on the trading activity and price distribution and interpret our results in light of the market efficiency hypothesis. The approach is directly applicable to other European electricity markets.


Introduction
Intraday electricity markets are used to balance the short-term intermittency of renewable energy assets. The increasing penetration of wind and solar yields the need for accurate probabilistic price forecasts. This need is  underscored by the heightened volatility in light of the European energy crisis in 2022/23. To date, the literature on (probabilistic) electricity price forecasting on intraday markets remains scarce [33,20,48,32] and is exclusively focused on univariate approaches. However, in most European countries, the intraday market is a continuous forward market, where all delivery periods for a delivery day d are traded in parallel, thus univariate approaches neglect the complex dependency structure in these markets. As products close with the physical delivery of electricity, it is not possible to "glue" subsequent trading sessions together, as it is commonly done with equity markets. Additionally, while the spot market is driven by the absolute level of fundamentals such as wind and solar production, the intraday prices are influenced more by the changes in forecasts [50]. Our work focuses on these challenges in modelling the dependency structure and fills the according gap in the literature. Based on the work of [33,20] on the marginal distribution of the price process, we employ copulas to model the time-dependent correlation. We model the dependency parameter as latent variable using beta regression. We validate our results in a forecasting study for the German intraday electricity market. Our results indicate that modelling the dependence structure between the different trading session improves the forecasting performance. Additionally, we provide evidence on market efficiency in the German short-term market. Our fundamental and parametric approach allows us to shed further light on the impact of the cross-border shared order books of the single intraday coupling (SIDC) and the driving factors of the distribution parameters.
Let us give an illustrative example for the trading schedule of the German intraday electricity market. Trading for physical delivery on d starts on the previous day d − 1, 15:00 hours and lasts till few minutes before the start of physical delivery of electricity. For example, trading for the delivery on d, 18:00-19:00 started on d − 1 at 15:00 hours and closes at d, 17:55 hours. During the beginning of this trading window, traders will be able to trade power in many neighbouring delivery hours, from delivery at d 16:00-17:00 (which closes at d − 1 15:55), to all delivery periods for the next delivery day d + 1 (which start trading at d, 15:00 hours). Figure 1 shows all trades in the two trading sessions for delivery on 2022-12-01 and 2022-12-02 for all hourly products on the German intraday market.
The literature on intraday electricity price forecasting can be divided into three main groups: (1) paper which treat the intraday market in a similar fashion as the day-ahead market and predict (index) prices along the delivery time line [48,9,23] and (2) paper which predict prices along the trading time for single delivery periods [42,20,32,33,21,30]. The correlation between different trading windows has, to the best knowledge of the authors, not been investigated so far, while the correlation between day-ahead and intraday (index) prices has been the subject of studies such as [3,9]. Lastly, (3), empirical, in-sample studies on price formation in intraday electricity markets have been conducted by [24,25,22].
In light of the aforementioned challenges and the gap identified in the literature, our contributions are: • We develop a global model for the marginal distribution of intraday electricity prices in the German market and extend previous work from [33,20] by taking the whole trading window into account.
• As novelty, we analyse the correlation and dependence structure in the German intraday electricity market and develop a multivariate, probabilistic forecasting model for the German intraday electricity markets to take cross-product effects into account.
• We validate our approach in an extensive simulation study for the German intraday electricity market.
• Using a parametric approach, our methods and models shed light on the driving variables such as trading activity and renewable forecasts in the intraday market, but also on the impact of the market structure and SIDC.
Our main strategy is the canonical inference-for-margins [37] approach commonly used for copula modelling and can be summarized as follows: We use the probabilistic models developed by [33,20] as a starting point to gaussianize the intraday market observations. We use the pseudo-Gaussian observations to fit the dependency structure. For forecasting, we simulate (multivariate) Gaussian random variates and use the inverse probability integral transformation to receive samples in the desired marginal distribution. Within the energy markets literature, similar approaches have been used for simulating wind and load forecast deviations [47,7,8], day-ahead electricity prices [see e.g. 28,39,5] and the design of hedging strategies [40]. Our results show that modelling the dependence structure leads to significantly improved forecasting performance compared to univariate approaches. However, we find that time-dependent modelling of the dependence structure is of little added value compared to constant dependence. Additionally, we provide new insight on the effects of the opening and closing of the cross-border order books during SIDC. We interpret our results in light of the market efficiency hypothesis and discuss reflections on modelling already highly volatile prices during a period of increased uncertainty. Figure 2 gives an illustrative example of our 24-dimensional forecast Our results offer multiple avenues for further research. First, we restrict ourselves to Gaussian dependence structures. Future work might improve our methods by using copulas that reflect possible tail dependence effects or use Vine-copulas to better approximate the dependence structure. Secondly, we note that the zero-inflated Johnson's S U distribution offers a good, but not yet perfect fit for the marginal distribution in the intraday market and hence further research in the marginal distribution of intraday electricity prices is required. Third, in light of the continued development of SIDC and the planned introduction of intraday auctions within the SIDC system and the possible integration of interconnector cables in the SIDC system [14], our results on the impact of SIDC on the trading activity provide a fruitful starting point for further research [for an early work on the topic see 22].
What is more, our results are also relevant for researchers and practitioners working on stochastic optimization  of bidding strategies for the intraday markets. For storage assets such as batteries and pumped-hydro, modelling the dependency structure is important as charging (pumping) positions some periods depends on the ability to discharge (generate) later during the day and therefore depends on the dependency structure. Recent works as [6,35,17] use sampled paths for the intraday market, but model the dependency structure only implicitly, if at all, and thus might produce too optimistic results.
The remainder of this paper is structured as follows: the following Section 2 gives a detailed introduction to the German short-term electricity market. Section 3 introduces our data set, preparation and summary statistics. We present our modelling approach in Section 4. Section 5 describes the forecasting study design and scoring rules. Finally, Section 6 and 7 present and discuss our results and conclude this paper.

Market Description
Electricity markets are structured as forward markets. The German short-term electricity market consists of three major parts: (1) the daily spot auction, (2) the continuous intraday market and (3) the balancing market. The following discussion focuses on the spot and intraday markets. The spot market is the main electricity market in Germany. It is organized as a daily auction at noon on which electricity for all 24 delivery hours on the following day is traded. The following intraday market is used to balance deviations in forecasts after the day-ahead market. It is organized as continuous trading similar to equity markets. The continuous trading starts at d − 1, 15:00 hours and closes shortly before delivery. After the delivery period ends, remaining imbalances between the traded position and the actual production are settled in the balancing market with the TSO. However, strict market regulation in Germany prohibit explicit active position taking in the balancing market. Figure 3 depicts the daily procedure for a single delivery hour.
Let us introduce some nomenclature to ease the following discussion of the market structure of the spot and intraday market. We refer to the delivery time as the time of actual production power, while the trading time refers to the time at which a trade for a certain delivery period is conducted. As a general rule, we try to denote delivery time in superscript, while we denote trading time in subscript. We refer to a trading session on the intraday market as  The spot market in Germany is organized by EPEX Spot and Nordpool AS with shared order books. On the spot market, electricity for all 24 hours for the following day is traded. The market is organized as a pay-as-cleared auction and its order book closes at d − 1, 12:00 hours. Results are published at d − 1, 12:42 hours. The minimum price is currently set to -500 EUR/MWh and the maximum price is set to 3000 EUR/MWh.   [15]. Note that we have omitted trading windows during days to save space.
The intraday market is structured as continuous forward market. For all hourly delivery periods with delivery on d, trading starts at d − 1, 15:00 hours and ends 5 minutes before the actual start of delivery. Within the Single Intraday Coupling (SIDC), the order books of all major continuous intraday markets across Europe are coupled, as long as there is sufficient cross-border transmission capacity. The coupling proceeds in two waves: first, at d − 1, 18:00 hours, the order books of Germany, Denmark, Sweden, Poland and Norway and Netherlands 1 are coupled. At d − 1, 22:00, France, Netherlands, Belgium, Austria, the Czech Republic, Hungary, Romania follow [34]. All shared order books close 60 minutes before the start of physical delivery and trading resumes with Germany wide delivery. 30 minutes before the start of physical delivery, the Germany wide delivery and trading resumes on a TSO/grid zone-level. Finally, 5 minutes before delivery, the grid-zone trading closes as well. Note that all hourly (and also half-hourly and quarter-hourly) delivery periods for a delivery day d are traded in parallel, as it shown in Figure 4. A more detailed overview on intraday electricity markets can be found in [43] and [49].

Data
The following chapter gives a brief introduction of the data used in this paper and the required pre-processing. We use intraday transaction data from EPEX, the anonymous day-ahead spot auction bid curves from EPEX, and wind, solar and demand forecasts from SMARD respectively ENTSO-E. We also provide summary statistics. As a rule of thumb, superscript indices denote delivery periods while subscript indices denote trading time. We hope this makes the forward market structure of the short-term electricity markets more clear.
Intraday transactions are individual trades conducted on the continuous intraday market. We use only trades conducted with either the buy-or sell leg in one of the 4 German grid zones. As trading happens continuously in the intraday market, transactions are irregular spaced in time and need to be aggregated. In line with [20,33,32,42], we aggregate all trades for delivery period d, h on a 15-minute equidistant grid along the trading time (denoted with t), where t = 0 denotes the first 15 minutes after trading start. Let P d,h t denote the volume-weighted average price of all trades with delivery period d, h belonging to bucket t and α d,h t denote a Boolean indicator whether there was at least one trade. The spot price for each delivery period is denote as P d,h Spot . We drop all trades in the local trading phase in the last 30 minutes to the start of physical delivery. The full aggregation process can be seen in Figure  5. Summary statistics for the price differences are given in Table 1. We especially note an increasing volatility in year 2022 driven by the Russian invasion in Ukraine and the energy crisis in Europe. Additionally, we note that trading activity in general increases as the share of no-trade event decreases. Figure 6 gives the pairwise correlation between price changes for all 24 × 24 intraday delivery periods during the training set. We already note a cluster of high correlation for the night hours and for the afternoon peak hours. MAD denotes the median absolute deviation and IQR denotes the interquartile range.  We use wind on-and offshore, solar and demand forecasts from ENTSO-E. The data is aggregated to hourly resolution using a simple arithmetic average. Forecasts are generated by the transmission system operator for each delivery period d, h and available at the day-ahead stage, i.e. latest at d−1, 12:00 o'clock. We denote the forecasts as WindOn d,h , WindOff d,h , Solar d,h and Load d,h . For all data, we adjust the daylight saving times by (back-) filling the missing hour in spring and averaging the double hour in autumn as it is standard in the electricity price forecasting literature.     Additionally, we employ a metric for the merit-order regime. In the classical model, the merit-order is defined as the supply side of the electricity market, sorted by the marginal production costs. The intersection of the supply and demand curve gives the market price. Depending on the slope of the merit-order, changes in the supply or demand have different impacts on the market price. [25] and [20] have shown that the different merit-order regimes explain the size and volatility of price changes in the intraday markets. However, modelling the merit-order in short-term power markets is not straight forward and several approaches have proposed. We follow the approach of [20] by using anonymous bid and offer curves from the day-ahead market to model the intraday merit order. The curves are published by EPEX Spot around 14:30 and are therefore available at the time of forecasting.
The overview of the strategy for the calculation of the merit-order regime coefficient MO d,h is given in Figure 8. We take the anonymous supply and demand curves from the auction and transform these into an elastic supply curve and an inelastic demand curve using the same transformation as in [20,26]. This transformation is based on the idea that buying 50 MW up to a price of 100 EUR/MWh is the same as buying 50 MW at any price and placing a sell order at 100.01 EUR/MWh. Intuitively, we can therefore create a price independent buy curve and move the price-dependent bids to the sell side. The resulting demand and supply curves are depicted in the middle panel of Figure 8. Note that the equilibrium price at the intersection between supply and demand curve is unchanged. Lastly, we calculate the slope around the equilibrium price as finite difference quotient. For the exact calculation, we refer the reader to [20].

Models
Our approach follows the widely used inference for margins approach for copula models. We first use a univariate model to estimate the time-varying conditional marginal distribution, apply the probability integral transform and estimate the copula distribution respectively the time-varying dependence parameter in the second step. The following two subsections 4.1 and 4.2 describe our modelling in more detail. Section 4.3 describes the general simulation set-up. Subsection 4.4 describes the different (nested) models for our forecasting study and our benchmark models.
Let us remark that strictly, Sklar's theorem [44] is only valid for continuous marginal distributions, as F X (X) ∼ U [0,1] is not true for discrete distributions F. This leads to the issue that the copula distribution is not necessarily identifiable with discrete or mixed discrete-continuous marginals. A practical approach to alleviate this issue is to fill the gaps induced through the discreteness of the marginal distribution by some uniform 'jitters', which results in the so-called checkerboard copula, a strategy we employ in this paper as well [18,16].

Marginal Model for Electricity Prices
We model the distribution of the price changes ∆P d,h t as a mixture distribution D d,h t to account for the zero-inflation in price changes.
is a mixture distribution from a continuous distribution and the Dirac distribution with an atom at 0, denoted as δ 0 .
where α d,h t = 1 indicates a trade event and is modelled as . A similar approach is taken by [33,20], who estimate the mixture distribution in a two-stage procedure. Owing to stylized facts on intraday electricity prices, [33,20] assume F to follow the (skewed) Student-t or Johnson's S U distribution. With regards to [20] remarks on the issues with estimation stability using the skew t-distribution, we generally use Johnson's S U distribution in this work. We use a two-step estimation procedure for the estimation. First, we estimate the conditional mixing probability π d,h 5. Trading variables. We use first three lagged prices, absolute lagged prices and the first three lagged values of α t to account for the auto-regressive nature in the distribution parameters. We also use a spline for the level of the spot price, denoted as ReLU(P d,h Spot ). ReLU splines are piecewise linear approximations on the domain of an explanatory variable. For an arbitrary explanatory variable x we define a set of thresholds τ ∈ T x and we define the ReLU splines as thereby the clipped values on the domain can contribute with individual slope coefficients. This allows an efficient approximation of non-linear functional relationships while preserving linearity in the coefficients [31]. ReLU splines are especially useful for variables with known and bounded domain such as the trading time t in our application, but can in theory be used for any continuous variable. Even though the thresholds T x can be chosen arbitrarily, we generally use an equidistant grid.
The logistic regression model used to estimate the conditional mixing probabilities π d,h t is defined as Dummies for SIDC. See Table 2 for the exact definitions Day-of-the-week dummies.
ReLU splines for trading time per delivery hour.
where β 0 denotes the intercept, β 1 to β 4 are the coefficients for the fundamental wind, solar and load forecasts, β 5 to β 9 are the coefficients for the single intraday coupling dummies SIDC {O,1,2,C,L} (d, h, t). The dummies denote whether the cross-border order books are open, the start of the 1st and 2nd wave, the closing of the cross-border order books 1 hour before the gate closure and the last periods of (local) trading.
Closing of the cross-border order books one hour before delivery.
All periods after the closing of cross-border order books. of the SIDC related dummies. We include day-of-the-week dummies for the delivery day d with the coefficients β 10 to β 15 . The last term denotes a ReLU spline for the trading time t on an equidistant grid T t , with ⊗ denoting the Kronecker product. We conduct a grid search on the validation data set to select the step size of the grid T t and the regularization parameter based on the validation accuracy. The step size is constant across all h. An exemplary fit for the ReLU spline along the trading time t can be seen in Figure 9. The estimation is regularized using the L 2 norm, which we choose from an exponential grid on 0.001 to 1000. We also evaluate the Bayesian Information Criterion (BIC) for each combination of step size and L 2 regularization and find that the results align.
The following paragraph describes the probabilistic model for the marginal distribution F with time-varying location, scale and shape parameters. Our general framework follows the generalized additive models for location, scale and shape introduced by [41]. Let Y = (Y 1 , Y 2 , ..., Y n ) be a vector of n independent observations Y i and Y i have the probability (density) function where each distribution parameter can be a smooth function of explanatory variables. Denote θ k i = (µ i , σ i , ν i , τ i ) = (θ 1 , θ 2 , θ 3 , θ k ) as the n × k parameter vector with k location, scale and shape parameters. We have Let g k (·) be a known monotonic link function for each distribution parameter relating θ k to explanatory variables through an additive model where X k is a n × J k known design matrix of J k exogenous regressors. Additionally, the model can consists of non-linear effects as in classical generalized additive models. Note that each distribution parameter can have an individual design matrix X k .
As noted, we assume Johnson's S U distribution for the marginal distribution of ∆P d,h t . The probability density function of Johnson's S U distribution is defined as follows: where −∞ ≤ µ ≤ ∞, 0 < σ ≤ ∞, 0 < ν ≤ ∞, −∞ ≤ τ ≤ ∞ represent the location, scale, tail and skewness parameters and their respective domains. To ensure all distribution parameters are within their domain, we employ the link functions g σ (x) = ϵ + log(1 + exp(γ · x)), g ν (x) = ϵ + log(1 + exp(γ · x)), where the link functions for µ and τ are the identity and the link functions for σ and ν are known as Softplus link function with constants ϵ = 1 −3 and γ = 0.1 to improve numerical stability [45]. Formally, we define the model as follows: Dummies for SIDC.
Trading time t ReLU spline.
Time to delivery T − t ReLU spline.
Dummies for SIDC.
Day-of-the-week dummies.
We model the conditional location, scale and shape parameters based on fundamental variables. We guide our selection from the literature [20,33,21,32] and regularize the estimation to avoid overfitting.
• The location parameter µ is modelled through three autoregressive lags of ∆P d,h t . • The scale parameter σ is modelled through the the fundamental wind, solar and demand forecasts and the coefficient for the merit-order regime. Additionally, we include day-of-the week and hourly dummies to account for different baseline volatility. We include a ReLU spline for the trading time and the time to delivery to capture the effects of increasing volatility towards gate closure and a ReLU spline for the spot price level to account for the impact of different price regimes.
• The tail parameter ν is explained using a similar, but slightly reduced set of features as the scale parameter σ.
• The skewness parameter τ is explained only using day-of-the-week and hourly dummies. This characterization is backed by the recent findings of [20].
As the modelling of the distributions higher moments is dependent on the quality of the lower moments' models [on this issue, see. e.g. 51], we refrain from making the models for ν and τ as complex as the model for the scale parameter σ and decrease complexity accordingly. For the expected value of intraday price changes, many studies have shown indications of weak-form market efficiency [20,32,33,48,21,27], hence we keep the model simple as possible without comprising the quality of the estimation for σ, ν and τ .
Our estimations are implemented using the scikit-learn [38] and the tensorflow-probability packages [13,1] in python. We automatically tune the hyperparameters of our model using the well-known optuna library, a Bayesian framework specifically tailored towards the optimization of hyperparameters of machine learning models [2]. The sampling space for the hyperparameter tuning framework can be found in Table 3. We initialize the coefficients for the estimation of the conditional location parameter µ as 0, while the remaining coefficient vectors are initialized uniformly sampled. We tune the L 1 -regularization for all distribution parameters, the learning rate and introduce a dropout layer during the training process to reduce the risk of overfitting [4,46,51]. During model fitting, we reserve 25% of our training data as validation set to employ early stopping if the training-validation loss does not improve after 25 epochs [51,29]. We run 250 iterations of the optuna algorithm and observe the best trial at iteration 48. Diagnostic plots can be found in the Appendix (see Figure 15).

Dependence
Remember that we simulate the T ×H price difference vector ∆P d t = (∆P d,1 t , ∆P d,2 t , ..., ∆P d,23 t ), where the different delivery periods can be correlated. Our strategy follows the general setup of [7,8]. We estimate the dependence as the covariance after we gaussianize our data in the spirit of Sklar's theorem.
In our copula-based modelling approach, we employ three different structures for the dependence with increasing complexity. The simplest model assumes that all delivery periods are independent and is denoted as Mix.Ind. The next complex model assumes a constant cross-product dependence across the trading window and is denoted as Mix.CD. Lastly, we estimate a time-varying dependence parameter across the trading time t. The resulting most complex model is denoted as Mix.TD.
We denote the dependence parameter between two delivery hours A, B as ρ A,B . For the constant dependence model, we estimate the pairwise dependence parameter as: The pair-wise estimation is necessary as not all delivery periods have the same trading length. For our time-varying dependence model Mix.TD, we fit ρ as a latent variable using the beta-regression. , (0, 12, ..., 128)).
To estimate the beta-regression, we use the correlation coefficient on the pseudo-Gaussian observations, which we map to the (0, 1) space required by the beta-regression model. We use the BIC to decide on the width of the grid of the spline on t and transform the estimated values back to the covariance using the inverse mapping and the empirical standard deviation of the pseudo-Gaussian observations.

Simulation Approach
We where it is important to note that the end of trading, T , depends on h. Therefore, the vector is not square but has the asymmetric trapezoid shape already visible in Figure 1 and

Benchmark Models
We employ three benchmark models from the literature. First, we employ the well performing naive model introduced by [33,20]. The model re-uses past price trajectories by sampling. We employ the model in two versions, first assuming independence between different delivery hours (i.e sampling each delivery hour h = 0, ..., H individually) and second by sampling the full T × H path vector. The third benchmark model is a arithmetic random walk model introduced by [27] drawing from the empirical price difference distribution. The following paragraphs introduce the benchmark models formally.
The Naive.Ind and Naive.Dep model is defined as where d ′ is a random day sampled from the training data set d ′ ∼ U({0, ..., d − 1}). Note that for the Naive.Ind we sample d ′ independent for each delivery hour h = 0, ..., H and for the Naive.Dep, we use the same d ′ for all h. [20,33] have shown that this type of benchmark models provides very good point and probabilistic forecasting performance in ensemble forecasting settings.
The RW.Emp has been used by [27] to generate paths for the intraday market in an application study for grid-scale storage optimization. The price process is defined as random walk, where the innovations are drawn from a discrete distribution of the centered empirical price changes. It is defined as: where again d ′ is a random day sampled from the training data set d ′ ∼ U({0, ..., d − 1}) and the second term centers all residuals to mean zero, ensuring the price process to be a martingale.
For all three benchmark models, the size of the training data d = 0, ..., d − 1 is a tuning parameter, which we optimize through a grid search.

Forecast Study and Evaluation Metrics
We employ the well-known rolling window forecasting study design. Our data comprises of the three years between 2020-01-01 and 2023-01-01. We split our data set in a training, validation and test data set. Our initial training set is 2020-01-01 to 2021-07-01, the validation set contains 2021-07-01 to 2022-01-01 and our test set contains the final year 2022-01-01 to 2023-01-01. The split is also depicted in Figure 7. We use the initial training set to develop and estimate our models and the validation data set to calibrate hyperparameters. We use the test set to run a rolling window forecasting study with monthly re-training of all models using the most recent two years of data. This is due to the high computational burden of training and tuning the probabilistic models in tensorflow-probability.
We evaluate our simulations from a point and probabilistic forecasting perspective using strictly proper scoring rules [19]. We evaluate the mean and median simulation paths using the well known root mean squared error (RMSE) and median absolute error (MAE). The marginal fit is evaluated using the continuous ranked probability score (CRPS), which we approximate on a dense grid of quantiles using the pinball score (PB). We evaluate the scenario paths using the energy score (ES). The energy score is the multivariate generalization of the continuous ranked probability score, taking into account the correlation structure. We establish significance using the Diebold-Mariano test [12,11] for comparing predictive accuracy. The following few paragraphs introduce our metrics formally.
The root mean squared error is defined as: The RMSE is a strictly proper scoring rule for the expected value. We also report the RMSE averaged additionally across the delivery hours h and delivery days d.

Denote the median trajectory over all M simulated trajectories as med(P d,h,[m] t
). The mean absolute error is defined as: We evaluate the probabilistic forecasting performance using the continuous ranked probability score (CRPS) and the energy score (ES). Both are strictly proper scoring rules for the marginal distribution respectively the multidimensional predictive distribution [19,52] and routinely used in probabilistic energy forecasting [20,33,5,36].
The CRPS is approximated using the pinball score on a dense grid of quantiles of our scenarios. Let Q d,h τ,t (P . The τ -PB can be defined as: The CRPS is the average across all quantile levels T = {0.01, .., 0.99} of length 99: For the energy score, we implement the K-band estimator as given in [52]: for an integer 1 ≤ K ≤ M and where we set P . ∥·∥ 2 denotes the L 2 or Euclidean norm. We use K = 10 as trade-off between computational complexity and estimation accuracy. We evaluate the energy score for the full scenario trajectory and for the last three hours of before the start of physical delivery for each model. The later evaluation acknowledges the importance of the last hours of trading and appeals to practitioners in the field. We use the Diebold-Mariano (DM) test to compare the predictive accuracy of the forecasts [36,12,11]. Intuitively, the DM-test evaluates the null hypothesis (H 0 ) that the difference in means between the loss series of two models is statistically significantly different from zero. Formally, for two models A and B, let L A and L B the loss series. The loss differential ∆ A,B is defined as for the i-norm. For each model pair, we test two one-sided tests for the null hypothesis (1) E ∆ i A,B > 0 and (2) E ∆ i A,B < 0, i.e. (1) the forecasts of model B outperform the forecasts of model A and (2) the forecasts of model A outperform the forecasts of model B. These tests are complimentary. Note that the Diebold-Mariano test assumes the loss differential series to be stationary. We test this assumption using the augmented Dickey-Fuller test [10].

Results and Discussion
Our parametric approach to modelling the distribution parameters allows us to derive some fundamental insight in the driving factors of the intraday price process. The following section presents first presents some in-sample results from our modelling and subsequently presents the results from our forecasting study.

Fundamental Analysis
Our parametric modelling set-up allows us to analyse the influence of the driving factors for the location, shape and scale parameters of the price distribution. Our focus here is on the impact of the SIDC on the trading activity and the impact of fundamental variables on the volatility.
The evolution of the trading activity respectively the share of no-trade events can be seen in Figure 9. We note that in the first hours of trading and that trading activity rises non-linearly towards gate closure. The opening of SIDC induces to spikes in trading activity when the cross-border order books are coupled. At this point, orders in markets with previously different price levels are instantly matched, leading high trading activity for short periods of time. After the SIDC coupling, trading activity increases towards the gate closure. The last trading periods have (almost) no no-trade events.
We show the contribution of individual groups of regressors on the scale parameter in Figure 10. Recall that we model the volatility by four main groups of regressors: fundamental forecasts, time-derived variables, SIDC-related variables, and variables related to the trading activity. Generally, we note a pattern of higher volatility in the beginning of the trading session, followed by decrease and an increase closer to delivery. We note that SIDC has a distinct impact on the scale parameter: The opening of the cross-border order books at 18:00 and 22:00 for the leads to clearly visible spikes in the volatility. Subsequently, the phase during which the order-books are coupled is characterized by lower volatility. The closing of the cross-border order books shortly before delivery leads to a spike in volatility. The dampening effect of the open SIDC order books on volatility is likely due to the increased liquidity available to market participants, while the volatility spikes during the opening at 18:00 and 22:00 hours can be explained by matching the order books at different price levels. Overall, the effects contradict the results of [22], who finds no effects of SIDC, but align with [33,20] on the effect of the SIDC closing period. Comparing the different delivery hours h = 0, 6, 12 and 18, we note the different impact of the time to delivery and trading time splines. The time-to-delivery spline kicks in early, but keeps rather constant shortly before the delivery. On the other side, the trading time t spline rises with trading time.
The influence of fundamental variables like wind, solar and demand forecasts is small. One reason for this might be, that these forecasts are generated at the day-ahead stage and are not updated throughout the trading window, as previous works [25,50,20] have shown that the intraday market is more impacted by forecast changes compared, while the level of forecasts is less important. Additionally, the conclusion that fundamental variables seem not to improve forecasts supports the notion of market efficiency as already indicated by [32,33,20].

Forecasting Performance
This section presents the results of the out-of-sample forecasting study. Aggregate error metrics are given in Table  4. Figure 11 presents the error metrics by the delivery hour h. Figure 13 gives the results of our pairwise Diebold-Mariano tests.
Remember that we have both, the naive model and the mixture model in in at least two versions: one version assuming independence between different delivery hours and at least one version that considers the correlation structure. Additionally, the RW.Emp considers the correlation structure implicitly. The results for the energy score show that considering the correlation structure leads to (significantly, see Figure 13) better forecasts than assuming independence if the remaining model structure is unchanged. This holds for moving from the Naive.Ind to the Naive.Dep and for moving from the Mix.Ind to Mix.CD. Interestingly, modelling the dependence structure in a potentially time-varying fashion does not improve the forecasting performance. For the energy score, the difference in ordering for the models between the full path and the last three hours of trading. Note however, that the scale of both is not directly comparable. This is an important result for modelling the intraday market in    applications such as battery/storage optimisation, where the neglecting the correlation structure can therefore lead to too optimistic results [35,27].
On an aggregate level, we see that the RW.Emp yields the best point forecasting performance (MAE and RMSE) and the Naive.Dep yields the best probabilistic forecasting performance. The Diebold-Mariano test confirms the statistical significance of superior probabilistic forecasting forecasts for the Naive.Dep. We note that the forecasting performance of the mixture models is not as good as the rather simple benchmark models. Additionally, we note that the Mix.Ind exhibits a lower CRPS than the Mix.CD and Mix.TD, which suggests some some cross-propagation of errors. On the other hand, the Mix.Ind yields worse scores for the energy score than its sister models including a dependence structure.
The hourly shape of forecasts errors throughout the day is depicted in Figure 11. It follows the typical shape of prices in electricity prices, we see low forecast errors in the morning and higher errors through the day. Let us note that the benchmark models exhibit a slightly higher MAE during the afternoon peak, but show lower RMSE and ES during the same periods. Figure 12 presents the PB across the quantile range and the delivery hour. We can see that the highest errors occur during the afternoon peak hours and in the higher distribution quantiles. Note that the pinball loss only measures the marginal fit to the distribution. In an analysis of the in-sample, transformed observations we also note that throughout the rolling window study, the calibration decreases and we experience an underdispersed forecast (see Figure 14 in the Appendix).
Our results also emphasize the importance of robust approaches to modelling and forecasting in periods of high  volatility and black swan events such as the Russian invasion of Ukraine and the subsequent energy crisis in Europe 2022/23. The Naive.Dep and Naive.Ind model already have shown very good probabilistic forecasting performance in the respective studies of [20,33] and are robust to extreme events, as out-of-support situations with extreme prices cannot happen for these models. Additionally, our results can be viewed in the light of the market efficiency hypothesis. Our result indicate that including more data, especially data that does not change throughout the trading period (such as day-ahead forecasts) does not improve the modelling or the price. The superior performance of the benchmark models, especially of the RW.Emp, which ensures the Martingale assumption, underscores this notion. Similar results with respect to market efficiency have been found by [32].

Conclusion
This paper presents a a forecasting study for multivariate, simulation-based forecasting for intraday electricity markets. We provide insight in the dependence structure of short-term electricity markets and extend previous works of [33,20] to include cross-product price effects.
We develop a probabilistic model for the marginal distribution of the intraday price path, accounting for the impact of fundamental driving variables on the location, scale and shape parameters. As novelty, we employ copulas to model the (time-dependent) dependency between different delivery periods and the according parallel trading sessions and allow the dependency parameter to be time-dependent. We validate our results in a forecasting study for the German intraday electricity market. Our results indicate that modelling the dependence structure between the different trading session improves the forecasting performance. Additionally, we provide evidence on market efficiency in the German short-term market. Our fundamental and parametric approach allows us to shed further light on the impact of the cross-border shared order books of the single intraday coupling (SIDC) and the driving factors of the distribution parameters. Our case study employs data from the German intraday electricity market, but our method is directly transferable to other European electricity markets.
While we are able to show that modelling the dependence structure improves the forecasting performance, our methods can surely be improved and our results offer a multitude of further research areas: First, a further investigation of the dependency structure seems worthwhile. However, the improved modelling of the correlation structure is dependent on having a suitable probabilistic marginal model at hand. We note that the proposed mixture of Johnson's S U distribution still does not provide a perfect fit and struggles to cope with periods of extreme volatility. Hence, the modelling of the price distribution is an important field for further research. Third, we provide new evidence on the impact of the single intraday coupling (SIDC) on the price distribution and trading activity. As the SIDC system is dynamically changing and little researched, this might be an interesting direction for further research into the micro-structure of intraday electricity markets. The literature on intraday markets is still scarce compared to the fast growth of renewable energy sources and intraday electricity markets.    : Diagnostic plots for our hyperparameter tuning using the optuna framework. We show the bivariate distribution of the explored hyperparameters and the scatter plot of each hyperparameter towards the negative log-likelihood on the diagonal. Color represents increasing the trial number from dark to yellow. The red dot represents the best trial.