Next Article in Journal
Location-Based Optimized Service Selection for Data Management with Cloud Computing in Smart Grids
Next Article in Special Issue
Neural Network Based Model Comparison for Intraday Electricity Price Forecasting
Previous Article in Journal
Stationary Energy Storage System for Fast EV Charging Stations: Simultaneous Sizing of Battery and Converter
Previous Article in Special Issue
Modeling Intraday Markets under the New Advances of the Cross-Border Intraday Project (XBID): Evidence from the German Intraday Market
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Estimation and Simulation of the Transaction Arrival Process in Intraday Electricity Markets

House of Energy Markets and Finance, University of Duisburg-Essen, 45141 Essen, Germany
*
Author to whom correspondence should be addressed.
Submission received: 31 October 2019 / Revised: 20 November 2019 / Accepted: 26 November 2019 / Published: 27 November 2019
(This article belongs to the Special Issue Modeling and Forecasting Intraday Electricity Markets)

Abstract

:
We examine the novel problem of the estimation of transaction arrival processes in the intraday electricity markets. We model the inter-arrivals using multiple time-varying parametric densities based on the generalized F distribution estimated by maximum likelihood. We analyse both the in-sample characteristics and the probabilistic forecasting performance. In a rolling window forecasting study, we simulate many trajectories to evaluate the forecasts and gain significant insights into the model fit. The prediction accuracy is evaluated by a functional version of the MAE (mean absolute error), RMSE (root mean squared error) and CRPS (continuous ranked probability score) for the simulated count processes. This paper fills the gap in the literature regarding the intensity estimation of transaction arrivals and is a major contribution to the topic, yet leaves much of the field for further development. The study presented in this paper is conducted based on the German Intraday Continuous electricity market data, but this method can be easily applied to any other continuous intraday electricity market. For the German market, a specific generalized gamma distribution setup explains the overall behaviour significantly best, especially as the tail behaviour of the process is well covered.

1. Introduction

A consecutive growth of the number and volume of transactions in the intraday electricity markets has been observed since the introduction of these markets. This results in a higher concern of researchers regarding the intraday electricity markets. Both Uniejewski et al. [1] and Narajewski and Ziel [2] consider a very short-term point electricity price forecasting (EPF) of the ID 3 -Price index in the German Intraday Continuous market. In Andrade et al. [3] and Monteiro et al. [4], the authors conducted research regarding the electricity price forecasting in the Iberian intraday electricity market. They performed a probabilistic electricity price forecasting and a point EPF using artificial neural networks, respectively. Both Ziel [5] and Kulakov and Ziel [6] examine the impact of renewable energy forecasts, i.e., wind and solar energy, on the intraday electricity prices. The relationship between the fundamental regressors and the price formation in the intraday markets is studied by many other scientists, e.g., González-Aparicio and Zucker [7], Pape et al. [8], Kath and Ziel [9].
Due to the continuity of the intraday markets, an important aspect is the bidding behaviour of the market participants. This problem has been already examined by, among others, Kiesel and Paraschiv [10] and Aïd et al. [11]. In the following paper, we take a closer look at the transaction arrivals in the intraday electricity market. Figure 1 shows the trajectories of the counting processes that correspond to the transaction time arrivals. In the exercise, we assume that the transactions arrive in accordance with some time-dependent intensity function. Moreover, assuming the parametric probability distribution of the inter-arrival times, we can perform a maximum likelihood estimation of the parameters and then a meaningful forecasting study. In the exercise, we use a rolling window study, following the recommendations of Diebold [12].
Our approach to the transactions in the intraday electricity markets is, to the best of our knowledge, an innovative one despite its simplicity. In Graf von Luckner et al. [13], the authors modelled the intensities of the buy/sell orders in a more complex manner, but they do not perform any forecasting study.
The outcome of the study is very satisfying despite the simplicity of utilized methods. The paper contributes to the literature by
(1)
filling the gap in the literature regarding the intensity estimation of transaction arrivals and is a major contribution to the topic,
(2)
proposing a novel modelling approach for inter-arrival times using a time-varying generalized F distributions to capture the underlying uncertainties,
(3)
presenting a procedure to simulate trajectories from the estimated processes,
(4)
discussing functional evaluation criterion in forecasting studies for stochastic processes, especially counting processes,
(5)
presenting a rigorous forecasting study for the German Intraday Continuous electricity market data.
Note that the methodology can be easily applied to any other continuous intraday electricity market.
The paper is structured as follows. In the next sections, we present the setting and the modelling details. Then, we explain the estimation and simulation methods. In the next section, we briefly describe the German Intraday Continuous market. The sixth section introduces the forecasting study design and the evaluation measures are described in detail. Next, the empirical forecasting results are presented. The paper is concluded with a discussion of the results and further development possibilities.

2. Setting

In the majority of all continuous intraday markets there are S products traded for each day of delivery, e.g., S = 24 in a market with hourly products. For a certain day of delivery d a product s { 1 , 2 , , S } is traded in the trading period [ b ( d , s ) , e ( d , s ) ) prior to the beginning of the delivery. Here b ( d , s ) denotes the beginning of the trading and e ( d , s ) the end of the trading. Both b ( d , s ) and e ( d , s ) potentially depend on the day of trading d and the considered product s. However, in the majority of European intraday markets the end of trading does not depend on the product s but the beginning of trading does. Furthermore, if not mentioned otherwise all times are measured in hours.
During the trading period [ b ( d , s ) , e ( d , s ) ) we observe a series of n ( d , s ) intraday transaction times T d , s = ( T 1 d , s , T 2 d , s , T n ( d , s ) d , s ) satisfying b ( d , s ) < T 1 d , s , T i 1 d , s < T i d , s for i { 2 , , n ( d , s ) } and T n ( d , s ) d , s < e ( d , s ) . An example of trajectories of corresponding counting processes is presented in Figure 1 for S = 24 products in the German Intraday Continuous market. As mentioned, the beginning of trading time differs for each product. For instance, the trading period in the German Intraday Continuous market for the hourly product s = 1 with delivery starting at 00:00 is [ b ( d , 1 ) , e ( d , 1 ) ) = [ 9 , 0.5 ) and for the hourly product s = 24 with delivery starting at 23:00 is [ b ( d , 24 ) , e ( d , 24 ) ) = [ 32 , 0.5 ) .
Let us note that most of the transactions take place in the last hours of the trading period. The reason for this behaviour is the design of the intraday electricity market, i.e., its main purpose is to let the market participants react to the changes in production prediction. In the first hours of trading in the intraday market usually there is not much more information, when comparing to the day ahead market, but in the last hours before the delivery the difference is significant, and thus it is the most traded time period in this market. This pattern justifies the decision to parametrize the time in such a manner that the last hours of trading are indexed in the same way, disregarding the delivery time, as in Figure 1.

3. Modelling and Estimation

In the purpose of estimating transaction arrivals, we consider the series of inter-arrival times ( X i d , s ) = ( T i d , s T i 1 d , s ) , where i { 1 , 2 , , n ( d , s ) } , n ( d , s ) is the number of transactions on day d and product s and T 0 d , s = b ( d , s ) is the beginning of trading. As pointed out in Figure 1, only the latter hours of trading are of major interest for modelling. Hence, we focus on modelling only the part [ a ( d , s ) , e ( d , s ) ) by choosing a ( d , s ) such that b ( d , s ) < a ( d , s ) < e ( d , s ) . In the example Figure 1, a reasonable choice for a ( d , s ) could be, for example, 8 , 5 or 3 . Denote X d , s = ( X i d , s ) all inter-arrival times after a ( d , s ) , so they satisfy a ( d , s ) < T i d , s . Further, let l ( d , s ) be the smallest index such that a ( d , s ) < T l ( d , s ) d , s holds.
Now, we assume that the series of inter-arrival times X d , s is independent and follows a probability distribution with a parametric density function f X d , s ( x ; θ ) . Therefore, knowing that the inter-arrival times are independent, we can perform the maximum likelihood estimation of the unknown vector of parameters θ
θ ^ = arg max θ f X d , s ( x ; θ ) = arg max θ i = l ( d , s ) n ( d , s ) f X d , s ( x i d , s ; θ ) .
Naturally, to make the estimation less biased, we can estimate the parameters using more than one day of history of the transaction arrival times. Assuming the independence between them and that we estimate based on D days of history, we get the following maximum likelihood estimator
θ ^ = arg max θ f ( X 1 , s , X 2 , s , , X D , s ) ( x ; θ ) = arg max θ d = 1 D f X d , s ( x ; θ ) = arg max θ d = 1 D i = l ( d , s ) n ( d , s ) f X d , s ( x i d , s ; θ ) .
The maximum likelihood problem stated in Equation (2) is solved using Rsolnp package in R, which was implemented by Ghalanos and Theussl [14], based on the algorithm of Ye [15], which is the general non-linear augmented Lagrange multiplier method. Since the likelihood function may contain local maxima, it is very important to set correctly the lower and upper bounds and the starting parameters. The algorithm should handle with no big problem up to 10-parametric optimization, so in purpose of our study it is satisfactory. Nevertheless, the choice of the maximum likelihood optimization tool is not crucial as we have a low dimensional problem.
In the case of German Intraday Continuous market, we assume four distributions of the inter-arrival times: exponential, gamma, generalized gamma and generalized F. Each of the consecutive distributions is an extension of the previous one. The distributions are parametrized as follows:
  • exponential distribution Exp( λ ) with rate parameter λ > 0 ,
  • gamma distribution Gamma( α , β ) with shape and rate parameters α > 0 and β > 0 ,
  • generalized gamma distribution GenGam( μ , σ , Q ) with location, scale and shape parameters μ R , σ > 0 and Q R ,
  • generalized F distribution GenF( μ , σ , Q , P ) with location and scale parameters μ R and σ > 0 , and shape parameters Q R and P 0 .
The exponential and gamma distributions are well-know and thus do not need any special introduction. The exponential distribution has the following density function
f ( x ; λ ) = λ exp { λ x }
and the gamma distribution has the density function defined by
f ( x ; α , β ) = β α Γ ( α ) x α 1 exp { β x } .
Let us remind that the exponential distribution is a special case of the gamma distribution, i.e., if X Exp ( λ ) , then X Gamma ( α = 1 , β = λ ) .
The generalized gamma distribution is an extension of the gamma distribution by Stacy et al. [16], but in the study we use the parametrisation of Prentice [17], which is stated above. If γ Gamma ( Q 2 , 1 ) and w = log ( Q 2 γ ) / Q , then x = exp ( μ + σ w ) follows the generalized gamma distribution with probability density function
f ( x ; μ , σ , Q ) = | Q | ( Q 2 ) Q 2 σ x Γ ( Q 2 ) exp { Q 2 ( Q w exp { Q w } ) } .
The relationship between the gamma and generalized gamma distributions parametrized in such a manner is as follows. If X Gamma ( α , β ) , then X GenGam ( μ = log ( β / α ) , σ = 1 / α , Q = 1 / α ) .
The last and the most general distribution that we assume is the generalized F distribution described by Prentice [18]. Define s 1 = 2 ( Q 2 + 2 P + Q δ ) 1 and s 2 = 2 ( Q 2 + 2 P Q δ ) 1 , where δ = ( Q 2 + 2 P ) 1 / 2 . If w = ( log ( x ) μ ) δ / σ , then the probability density function of x is given by
f ( x ; μ , σ , Q , P ) = δ ( s 1 / s 2 ) s 1 exp { s 1 w } σ x ( 1 + s 1 exp ( w ) / s 2 ) ( s 1 + s 2 ) B ( s 1 , s 2 ) ,
where B ( s 1 , s 2 ) is the beta function. Let us note that if we possess a random variable X GenGam ( μ , σ , Q ) , then X GenF ( μ , σ , Q , P = 0 ) . We see clearly that all the consecutive distributions are superior to the previous ones. Figure 2 presents densities of exemplary GenF distributions. In the exercise, we use the flexsurv package in R by Jackson [19], which contains the implementation of the generalized gamma and F distributions.
We assume that the rate parameter λ , the rate and shape parameters β and α , and the location and scale parameters μ and σ are some deterministic functions dependent on time t and unknown vector of parameters θ . Thus, we model them using the following functions:
  • constant— f ( t ; θ ) = c , where θ = c ,
  • linear— f ( t ; θ ) = c + β 1 t , where θ = ( c , β 1 ) ,
  • quadratic— f ( t ; θ ) = c + β 1 t + β 2 t 2 , where θ = ( c , β 1 , β 2 ) ,
  • exponential— f ( t ; θ ) = c + e α 1 + α 2 t , where θ = ( c , α 1 , α 2 ) .
In the estimation of the parameters of generalized gamma and F distributions, we make use of the relationship between gamma and generalized gamma distributions, i.e., μ ( θ , t ) = log ( β ( θ 2 , t ) / α ( θ 1 , t ) ) and σ ( θ , t ) = 1 / α ( θ 1 , t ) . Moreover, we assume that α ( θ 1 , t ) cannot be more complex than β ( θ 2 , t ) and that the Q and P parameters are constant over time to delivery. By complexity of a function, we mean the number of parameters. Using this criteria, quadratic and exponential functions are equally complex.
The study consists in total of 37 models of the inter-arrival times process X d , s : four models with the assumption of the exponential distribution and 11 models for each other considered distribution. We abbreviate them by X.Y.Z, where X stands for the distribution and Y and Z stand for the types of the β ( θ 2 , t ) and α ( θ 1 , t ) functions, respectively. For instance, Gamma.Lin.Const stands for a model with assumed gamma distribution, linear β ( θ 2 , t ) and constant α ( θ 1 , t ) . Let us note that for the model with exponential distribution and λ ( θ , t ) = θ the corresponding counting process is a homogeneous Poisson process and this model is our basic benchmark. Making the intensity function λ ( θ , t ) non-constant is the first extension of the benchmark and changing the distribution to the more general one is the further extension. Moreover, in this study for non-constant rate and shape functions we actually assume that the functions are constant in short intervals, e.g., for exponential distribution, λ ( θ , T i d , s ) is assumed to be the intensity on the time interval [ T i d , s , T i + 1 d , s ) . This means that the corresponding counting process of a model with exponential distribution is a mixture of homogeneous Poisson processes.
Let us note that the models contain not more than eight parameters, so their estimation should not be a problem. Figure 3 shows an example of fitting the aforementioned models assuming the exponential distribution to one day of data for one product. The figure presents the observed trajectory and estimated cumulative intensity functions Λ ( B ) = B λ ( t ) d t . The time to delivery range is [ a ( d , s ) , e ( d , s ) ) = [ 3.25 , 0.5 ) , because in this particular exercise we aim to forecast the transaction arrivals during the ID 3 -Price period. This approach to the German Intraday Continuous was taken also by other researchers ([1,2]). Based on Figure 3, we may expect that considering the exponential distribution models, the models with quadratic and exponential intensity functions have similar performance in our problem.

4. Simulation

The trajectory simulation is relatively straightforward in our setting. Since we assume the distribution of the inter-arrival times X d , s = ( X i d , s ) , we can simply generate the inter-arrival times from the estimated distribution and calculate the next arrival time by adding the simulated inter-arrival time to the time of forecasting. The only part that is not obvious is the simulation of the first transaction arrival. That is to say, if we simulate at the time of the last observed transaction, there is no change, but if we fix the time of forecasting (as done in the study), then we have to truncate the distribution for the first observation at our starting point a ( d , s ) . This truncated distribution has the following density function
f X d , s θ , x i d , s | X i d , s > a ( d , s ) = g ( x i d , s ) 1 F X d , s ( a ( d , s ) ) ,
where g ( x ) = f ( x ) for all x > a ( d , s ) and g ( x ) = 0 otherwise, and F X d , s ( a ( d , s ) ) is the cumulative distribution function of X d , s . All the next transaction arrivals are simulated using the standard distribution with the density function f X d , s at the simulated times T ^ j , where j = i + 1 , i + 2 , until the desired end of forecasting.
An example of simulation of 20 trajectories in the German Intraday Continuous market is shown in Figure 4. The trajectories are simulated as described above, based on the generalized gamma distribution with quadratic α ( θ 1 , t ) and exponential β ( θ 2 , t ) functions, which was estimated using D = 28 days of data.

5. German Intraday Continuous Market

As mentioned, in the exemplary study we consider the German Intraday Continuous market. Trading in this market starts every day at 15:00 for hourly and at 16:00 for quarter-hourly products of the following day. Market participants can trade electricity until 30 min before the delivery in the whole market and until 5 min before the delivery in respective control zones. Figure 5 presents briefly the German electricity spot market, for more details see Viehmann [20].
An important measure in the German Intraday Continuous market is the ID 3 -Price index. The index is a volume-weighted price of all transactions taking place in the time interval between 3 h and 30 min before the delivery and it is calculated separately for each intraday product. The importance of the ID 3 -Price has been already noticed by the researchers and is a subject to modelling and forecasting by Uniejewski et al. [1] and Narajewski and Ziel [2]. In the latter paper, one can find a broader description and analysis of the ID 3 -Price and the German Intraday Continuous market. In both papers, the authors performed a very short-term point EPF of the ID 3 -Price. The outcome of Narajewski and Ziel [2] is the efficiency of the market, i.e., the volume-weighted price of the transactions in the last 15 min before forecasting appears to be the best model for the ID 3 -Price. In our study, we aim to forecast the time arrivals during the ID 3 -Price time interval. Naturally, we leave ourselves some time for calculation and decision-making and therefore the considered time frame is [ 3.25 , 0.5 ) which is exactly the same as in Narajewski and Ziel [2].
In Europe the majority of intraday electricity markets features a similar structure to the German intraday market. This holds especially for all markets that participate in the Cross-Border Intraday Project (XBID), see e.g., Kath [21]. It allows various participating electricity markets (e.g., Germany, France, Spain) to bid across borders, if inter-connector capacity allows doing so. From the modelling perspective it might be relevant to note that all electricity markets which participate in XBID close their markets for each product the same amount of time before delivery. Even though this is not always half an hour before delivery as in the German case an application to XBID markets should be no problem.

6. Data, Forecasting Study and Evaluation

In the following paper, as an example we perform a rolling window forecasting study based on the data from the German Intraday Continuous market. We consider a D = 28 -day window size with the initial in-sample data from 03.09.2017 to 30.09.2017 and forecast the next day arrivals, starting on 01.10.2017. We forecast the arrivals between 3 h 15 min and 30 min before the delivery. Our out-of-sample study is of size N = 365 , thus it spans the data range from 01.10.2017 to 30.09.2018. During each out-of-sample iteration M = 1000 trajectories are simulated. In the study, a multivariate approach is taken, which means that we create 24 separate models, each for every hourly product.
The data that we utilize in the study was obtained from EEX Transparency, and it consists of information regarding: the date of the delivery, the product type, market area, volume of the traded energy, price in EUR/MWh, the transaction ID and the time of the transaction. In our case, the only relevant informations are the date of the delivery, the product type and of course the time of the transaction. A small inconvenience regarding the data is the fact that the transaction times have minute grid. This makes many transactions have identical time arrival even if they weren’t made at the same time, e.g., four transactions with timestamp of 30.09.2017 16:01:00. We deal with the problem by distributing the transactions with the same timestamp T uniformly in the time range [ T , T + 1 min ) . Using the aforementioned timestamp as an example, the new timestamps are: 30.09.2017 16:01:00, 16:01:15, 16:01:30, 16:01:45.
Due to the lack of literature regarding the intensity estimation in the intraday markets, we cannot use any literature benchmark models or literature evaluation measures. Thus, in the study we simply compare the results of simulation of all the considered models and as evaluation measures we use the functional: bias (Bias), mean absolute error (MAE), root mean squared error (RMSE) and continuous ranked probability score (CRPS). We abbreviate the functional measures in a standard way, but the calculation is a little different. That is to say, let us denote by N d , s ( t ) the counting process of the true transactions on day d for product s and by N m d , s ( t ) the m-th simulation of the counting process N d , s ( t ) .
Let
ρ η , τ , p ( z ) = ( η z p + ( 1 η ) | z | p ) | τ 𝟙 ( z < 0 ) |
a loss function with η { 0 , 1 } , τ ( 0 , 1 ) , p 1 and 𝟙 as indicator function. Then we define
N ^ η , τ , p d , s ( t ) = arg   min z m = 1 M ρ η , τ , p ( N m d , s ( t ) z )
the sample ρ η , τ , p -estimate of the corresponding simulation sample N 1 d , s , , N M d , s . The special case ( η , τ , p ) = ( 0 , 0.5 , 2 ) in (8) corresponds to ordinary least squares (OLS) and ( η , τ , p ) = ( 0 , 0.5 , 1 ) to median regression. Thus, N ^ 0 , 0.5 , 2 d , s is the sample mean process, N ^ 0 , 0.5 , 1 d , s the sample median process and similarly N ^ 0 , τ , 1 d , s the sample τ -quantile process on which we focus especially. Now, we define the evaluation criteria
EVAL ( η 1 , τ 1 , p 1 ) , ( η 2 , τ 2 , p 2 ) d , s ( [ T 1 , T 2 ) ) = T 1 T 2 ρ η 1 , τ 1 , p 1 ( N d , s ( t ) N ^ η 2 , τ 2 , p 2 d , s ( t ) ) d t 1 p = j = 1 J ρ η 1 , τ 1 , p 1 ( N d , s ( t j ) N ^ η 2 , τ 2 , p 2 d , s ( t j ) ) Δ t j 1 p
on a time range [ T 1 , T 2 ) . J is the length of a grid of the time range [ T 1 , T 2 ) with t 0 = T 1 and t J = T 2 . The grid is defined by the jumps of both of the counting processes. Let us note that the transition from the integral to the sum is possible, because the difference of counting processes is a simple function. Moreover, we approximate the values of the evaluation measures by using a minute grid instead of the one defined by the jumps to reduce computational costs.
Now, we define the special cases of the evaluation measures which lead to the functional bias (Bias), the functional MAE, functional RMSE and functional pinball (PB) loss (or quantile loss) with respect to a probability τ :
Bias s = 1 N d = 1 N 2 EVAL ( 1 , 0.5 , 1 ) , ( 1 , 0.5 , 2 ) d , s
MAE s = 1 N d = 1 N 2 EVAL ( 0 , 0.5 , 1 ) , ( 0 , 0.5 , 1 ) d , s
RMSE s = 1 N d = 1 N 2 EVAL ( 0 , 0.5 , 2 ) , ( 0 , 0.5 , 2 ) d , s
PB τ , s = 1 N d = 1 N EVAL ( 0 , τ , 1 ) , ( 0 , τ , 1 ) d , s
We may observe that 2 PB 0.5 , s = MAE s . Additionally, we use the pinball loss to approximate the functional continuous ranked probability score (CRPS) by
CRPS s = 1 R τ r PB τ , s
for an equidistant grid of probabilities r between 0 and 1 of size R, see e.g., Nowotarski and Weron [22]. We choose r = ( 0.01 , 0.02 , , 0.99 ) of size R = 99 . In the purpose of comparing the models’ forecasting performance, we calculate the functional Bias, MAE, RMSE and CRPS based on M = 1000 trajectories in N = 365 out-of-sample iterations and for all s { 1 , , S } with S = 24 products. For creating summaries, it may be useful to average across all products and define
Crit = 1 S s = 1 S Crit s
for the considered criteria Crit { Bias , MAE , RMSE , PB τ , CRPS } .
To draw significant conclusions on the outperformance of the forecasts of the considered models, we also calculate the Diebold and Mariano [23] test, which tests forecasts of model A against forecasts of model B. In the following paper, we compute the multivariate version of the DM test as in Ziel and Weron [24]. The multivariate DM test results in only one statistic for each model that is computed based on the S-dimensional vector of losses for each day. Therefore, denote L d A = ( L d , 1 A , L d , 2 A , , L d , S A ) and L d B = ( L d , 1 B , L d , 2 B , , L d , S B ) the vectors of the out-of-sample losses for day d of the models A and B, respectively. By L d , s Z we mean the CRPS s loss for day d of model Z, formally we choose
L d , s Z = 1 R τ r EVAL ( 0 , τ , 1 ) , ( 0 , τ , 1 ) Z , d , s ( [ T 1 , T 2 ) )
with [ T 1 , T 2 ) = [ 3.25 , 0.5 ) . In the DM test we consider only the CRPS loss as it is the most important measure in our study. The multivariate loss differential series
Δ d A , B = | | L d A | | q | | L d B | | q
defines the difference of losses in | | · | | q norm, i.e., | | L d A | | q = s = 1 S | L d , s A | q 1 / q , where q { 1 , 2 } in our case. For each model pair, we compute the p-value of two one-sided DM tests. The first one is with the null hypothesis H 0 : E ( Δ d A , B ) 0 , i.e., the outperformance of the forecasts of model B by the forecasts of model A. The second test is with the reverse null hypothesis H 0 R : E ( Δ d A , B ) 0 , i.e., the outperformance of the forecasts of model A by those of model B. Let us note that these tests are complementary, and we perform them using two norms: | | · | | 1 and | | · | | 2 . Naturally, we assume that the loss differential series is covariance stationary.

7. Results

Table 1 presents the Bias, MAE, RMSE and CRPS measures (see Equation (15)) of the considered models on the interval [ T 1 , T 2 ) = [ 3.25 , 0.5 ) which was used for estimation. In the table, we observe that the lowest MAE and RMSE are obtained for model Gamma.Expon.Lin. At the same time we see that the values for the other models with gamma distribution and exponential rate function have similar performance, as well as model Exp.Expon, which appears to be the second best in terms of MAE and RMSE. These values indicate that most likely the difference in the performance of modelling the median and mean of models Gamma.Expon.Lin and Exp.Expon is not statistically significant. The models with generalized gamma and generalized F distributions clearly have difficulties in modelling the central parts of the distribution, but they handle the probabilistic forecasting well. The best performing model in terms of CRPS is GenGam.Quadr.Expon. The CRPS values of other models are mostly satisfying despite the models with constant parameters. These give the worst forecasts in terms of all considered measures.
On the other hand, the values of bias are interesting. We see there that the best models underestimate the true counting process. The reason for such an underestimation may be the fact that the German Intraday Continuous is constantly developing and the number of trades is growing every day. This suggests that there is still some space for improvement of the errors. In the remaining part of this section, we consider only selected best performing models: two best models per distribution and from each of the distributions, we choose also the best model with exponential rate function.
Figure 6 contains four graphs that present the performance of the selected models’ forecasts over products s { 1 , , S } . In Figure 6a, we observe that all the models despite the GenF.Lin.Const underestimate the true trajectory for all products. There is no clear pattern of the underestimation, but we see that the bias of best models in terms of MAE and RMSE is almost constant over products, whereas the bias of the best model in terms of CRPS varies significantly. Figure 6b–d show that all models handle the simulation of night hours’ transactions better than the simulation of the day hours’ transactions. An interesting spike in both measures can be observed for s = 6 . The reason for that might be some outliers as the spike is higher in terms of RMSE than in terms of MAE or CRPS. Furthermore, we see that the models with exponential and gamma distribution are clearly better than the others across all products in terms of MAE and RMSE. It is different for CRPS, where for most hours the best models are the GenGam models, but between hours 6 and 10 the gamma distributed models appear to be better.
Figure 7 is analogous to Figure 6, but the measures are calculated over the time to delivery. This means that we calculate the values (see Equation (15)) on many short time intervals, i.e., we apply a minute grid of the time-frame [ 3.25 , 0.5 ) . This way we can understand which part of the trajectory is forecasted the most and the least precisely. Figure 7a shows that the most underestimated is the last 1.5 h of the counting process, but the models with gamma and exponential distribution handle this period better than the GenGam and GenF models. Figure 7b,c show again the Gamma and Exp models with exponential rate function give mostly the lowest MAE and RMSE. In Figure 7d, we see that Gamma.Expon.Lin model has lowest CRPS in the first half of the forecasting horizon and, surprisingly, in the second half the lowest CRPS value is obtained for model GenF.Lin.Const. The performance of GenGam.Quadr.Expon is mostly stable over the time to delivery.
Figure 8 shows the pinball loss results of the models over probability values. Based on it, it is clear that the Gamma and Exp models forecast better the central quantiles, i.e., between 0.2 and 0.7, whereas the others are better forecasted by GenGam and GenF models. The figure indicates that the overall CRPS performance can be significantly improved.
To draw statistically significant conclusions, we perform a Diebold-Mariano test. The results are presented in Figure 9. Let us remind that we use the CRPS as the loss series as we believe that it is the most important measure in this study. Based on it, it is clear that the GenGam.Quadr.Expon model is significantly the best. Interestingly, the test in norm | | · | | 2 indicates a very good performance of the GenF models. Although, in this norm none of the other models is significantly better than the GenGam.Quadr.Expon.
Let us now take a look at the values of the time-varying coefficients of the model GenGam.Quadr.Expon. Let us recall that in this model we estimate the inter-arrival time X i d , s assuming generalized gamma distribution with the time-varying coefficients μ ( θ , t ) = log ( β ( θ 2 , t ) / α ( θ 1 , t ) ) and σ ( θ , t ) = 1 / α ( θ 1 , t ) and constant Q, where β ( θ 2 , t ) and α ( θ 1 , t ) are the rate and shape parameters of the standard gamma distribution. In this model β ( θ 2 , t ) is a quadratic function and α ( θ 1 , t ) is an exponential function. In Figure 10, we analyse the behaviour of μ and σ parameters over products at different time. Figure 10a shows the values of μ are similar over products, but differ significantly over time. That is to say, the closer we are to the gate closure, the lower the location coefficient, which means that the transactions appear more often. Figure 10b shows a very similar behaviour for the scale parameter σ . In this case the closer we are to the gate closure, the less vary the inter-arrival times.

8. Conclusions

We described the novel problem of estimation and simulation of the transaction arrival process in intraday electricity markets. The approach is not complicated and the presented methods are easy to implement. The paper fills the gap in the literature regarding the estimation and simulation of the transaction arrival process in the intraday electricity markets and thus is a major contribution to this field of research. The outcome of the conducted study is very satisfying.
Using the aforementioned approach, we utilized an exemplary rolling window forecasting study based on the German Intraday Continuous market. We assumed four probability distributions of the inter-arrival times: exponential, gamma, generalized gamma and generalized F distributions. We performed the maximum likelihood estimation of the distributions, assuming time-dependence of some of their coefficients. Then, using the estimated distributions we simulated new trajectories which we evaluated using the functional bias, MAE, RMSE and CRPS.
The results showed that the forecasting error can be significantly reduced, comparing to the most standard benchmark, which was the homogeneous Poisson process. The best in terms of forecasting of the central part of the distribution of the transaction arrivals, i.e., the mean or median are the exponential and gamma distributions with exponential rate function. In terms of the more meaningful and thus more important CRPS significantly the best forecasts were obtained from the generalized gamma model with quadratic rate function and exponential shape function.
This field of research can be easily developed further. A possible direction is considering other probability distributions of the inter-arrival times. Another possibility would be using more complex distributions’ parameter functions, such as Hawkes process-like, which is widely applied to modelling the transaction time arrivals in the financial markets, see, for example, Hewlett [25] or Bacry et al. [26]. The parameter functions could be also modelled using smoothing kernel or splines. To avoid overestimation, regularization methods should be considered.

Author Contributions

Conceptualization, M.N. and F.Z.; methodology, M.N. and F.Z.; software, M.N.; validation, M.N. and F.Z.; formal analysis, M.N. and F.Z.; investigation, M.N.; data curation, M.N.; writing—original draft preparation, M.N.; writing—review and editing, M.N. and F.Z.; visualization, M.N.; supervision, F.Z.

Funding

This research article was partially supported by the German Research Foundation (DFG, Germany) and the National Science Center (NCN, Poland) through BEETHOVEN grant no. 2016/23/G/HS4/01005.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Uniejewski, B.; Marcjasz, G.; Weron, R. Understanding intraday electricity markets: Variable selection and very short-term price forecasting using LASSO. Int. J. Forecast. 2019, 35, 1533–1547. [Google Scholar] [CrossRef]
  2. Narajewski, M.; Ziel, F. Econometric modelling and forecasting of intraday electricity prices. J. Commod. Mark. 2019, 100107. [Google Scholar] [CrossRef]
  3. Andrade, J.R.; Filipe, J.; Reis, M.; Bessa, R.J. Probabilistic Price Forecasting for Day-Ahead and Intraday Markets: Beyond the Statistical Model. Sustainability 2017, 9, 1990. [Google Scholar] [CrossRef]
  4. Monteiro, C.; Ramirez-Rosado, I.J.; Fernandez-Jimenez, L.A.; Conde, P. Short-Term Price Forecasting Models Based on Artificial Neural Networks for Intraday Sessions in the Iberian Electricity Market. Energies 2016, 9, 721. [Google Scholar] [CrossRef]
  5. Ziel, F. Modeling the impact of wind and solar power forecasting errors on intraday electricity prices. In Proceedings of the 2017 IEEE 14th International Conference on the European Energy Market (EEM), Dresden, Germany, 6–9 June 2017; pp. 1–5. [Google Scholar]
  6. Kulakov, S.; Ziel, F. The impact of renewable energy forecasts on intraday electricity prices. arXiv 2019, arXiv:1903.09641. [Google Scholar]
  7. González-Aparicio, I.; Zucker, A. Impact of wind power uncertainty forecasting on the market integration of wind energy in Spain. Appl. Energy 2015, 159, 334–349. [Google Scholar] [CrossRef]
  8. Pape, C.; Hagemann, S.; Weber, C. Are fundamentals enough? Explaining price variations in the German day-ahead and intraday power market. Energy Econ. 2016, 54, 376–387. [Google Scholar] [CrossRef]
  9. Kath, C.; Ziel, F. The value of forecasts: Quantifying the economic gains of accurate quarter-hourly electricity price forecasts. Energy Econ. 2018, 76, 411–423. [Google Scholar] [CrossRef]
  10. Kiesel, R.; Paraschiv, F. Econometric analysis of 15-minute intraday electricity prices. Energy Econ. 2017, 64, 77–90. [Google Scholar] [CrossRef]
  11. Aïd, R.; Gruet, P.; Pham, H. An optimal trading problem in intraday electricity markets. Math. Financ. Econ. 2016, 10, 49–85. [Google Scholar] [CrossRef]
  12. Diebold, F.X. Comparing predictive accuracy, twenty years later: A personal perspective on the use and abuse of Diebold–Mariano tests. J. Bus. Econ. Stat. 2015, 33, 1. [Google Scholar] [CrossRef]
  13. Graf von Luckner, N.; Cartea, Á.; Jaimungal, S.; Kiesel, R. Optimal Market Maker Pricing in the German Intraday Power Market; House of Energy Markets and Finance: Essen, Germany, 2017. [Google Scholar]
  14. Ghalanos, A.; Theussl, S. Rsolnp: General Non-Linear Optimization Using Augmented Lagrange Multiplier Method, R package version 1.16. 2015.
  15. Ye, Y. Interior Algorithms for Linear, Quadratic, and Linearly Constrained Non-Linear Programming. Ph.D. Thesis, Department of ESS, Stanford University, Stanford, CA, USA, 1987. [Google Scholar]
  16. Stacy, E.W. A generalization of the gamma distribution. Ann. Math. Stat. 1962, 33, 1187–1192. [Google Scholar] [CrossRef]
  17. Prentice, R.L. A log gamma model and its maximum likelihood estimation. Biometrika 1974, 61, 539–544. [Google Scholar] [CrossRef]
  18. Prentice, R.L. Discrimination among some parametric models. Biometrika 1975, 62, 607–614. [Google Scholar] [CrossRef]
  19. Jackson, C.H. flexsurv: A platform for parametric survival modeling in R. J. Stat. Softw. 2016, 70. [Google Scholar] [CrossRef] [PubMed]
  20. Viehmann, J. State of the German Short-Term Power Market. Zeitschrift für Energiewirtschaft 2017, 41, 87–103. [Google Scholar] [CrossRef]
  21. Kath, C. Modeling Intraday Markets under the New Advances of the Cross-Border Intraday Project (XBID): Evidence from the German Intraday Market. Energies 2019, 12, 4339. [Google Scholar] [CrossRef]
  22. Nowotarski, J.; Weron, R. Recent advances in electricity price forecasting: A review of probabilistic forecasting. Renew. Sustain. Energy Rev. 2018, 81, 1548–1568. [Google Scholar] [CrossRef]
  23. Diebold, F.; Mariano, R. Comparing Predictive Accuracy. J. Bus. Econ. Stat. 1995, 13, 253–263. [Google Scholar]
  24. Ziel, F.; Weron, R. Day-ahead electricity price forecasting with high-dimensional structures: Univariate vs. multivariate modeling frameworks. Energy Econ. 2018, 70, 396–420. [Google Scholar] [CrossRef]
  25. Hewlett, P. Clustering of order arrivals, price impact and trade path optimisation. In Proceedings of the Workshop on Financial Modeling with Jump Processes, Ecole Polytechnique, Palaiseau, France, 6–8 September 2006; pp. 6–8. [Google Scholar]
  26. Bacry, E.; Mastromatteo, I.; Muzy, J.F. Hawkes processes in finance. Mark. Microstruct. Liq. 2015, 1, 1550005. [Google Scholar] [CrossRef]
Figure 1. Trajectories of S = 24 counting processes corresponding to transaction times for day d = 03.09.2017 in the hourly German Intraday Continuous market. The dashed lines indicate the beginning of trading period b ( d , s ) = 8 s of product s { 1 , , S } with respect to the day of delivery d, the end of delivery e ( d , s ) = 0.5 .
Figure 1. Trajectories of S = 24 counting processes corresponding to transaction times for day d = 03.09.2017 in the hourly German Intraday Continuous market. The dashed lines indicate the beginning of trading period b ( d , s ) = 8 s of product s { 1 , , S } with respect to the day of delivery d, the end of delivery e ( d , s ) = 0.5 .
Energies 12 04518 g001
Figure 2. Illustration of the GenF distribution with different parameters.
Figure 2. Illustration of the GenF distribution with different parameters.
Energies 12 04518 g002
Figure 3. Observed trajectory of transaction counting process and fitted cumulative intensities of considered intensity functions of exponential distribution models for day d = 03.09.2017 and product s = 12 in the hourly German Intraday Continuous market.
Figure 3. Observed trajectory of transaction counting process and fitted cumulative intensities of considered intensity functions of exponential distribution models for day d = 03.09.2017 and product s = 12 in the hourly German Intraday Continuous market.
Energies 12 04518 g003
Figure 4. Observed (black) and simulated from the GenGam.Quadr.Expon(colourful) trajectories of transaction counting process for day d = 01.10.2017 and product s = 18 in the hourly German Intraday Continuous market.
Figure 4. Observed (black) and simulated from the GenGam.Quadr.Expon(colourful) trajectories of transaction counting process for day d = 01.10.2017 and product s = 18 in the hourly German Intraday Continuous market.
Energies 12 04518 g004
Figure 5. The daily routine of the German electricity market; d corresponds to the day of the delivery and s corresponds to the hour of the delivery.
Figure 5. The daily routine of the German electricity market; d corresponds to the day of the delivery and s corresponds to the hour of the delivery.
Energies 12 04518 g005
Figure 6. (a) Bias s , (b) MAE s , (c) RMSE s and (d) CRPS s of selected best performing models evaluated on [ T 1 , T 2 ) = [ 3.25 , 0.5 ) .
Figure 6. (a) Bias s , (b) MAE s , (c) RMSE s and (d) CRPS s of selected best performing models evaluated on [ T 1 , T 2 ) = [ 3.25 , 0.5 ) .
Energies 12 04518 g006
Figure 7. (a) Bias of the selected models, (b) MAE, (c) RMSE, (d) CRPS measures of selected best performing models in relation to the values of the model Exp.Expon and (e) errors of Exp.Expon over time to delivery.
Figure 7. (a) Bias of the selected models, (b) MAE, (c) RMSE, (d) CRPS measures of selected best performing models in relation to the values of the model Exp.Expon and (e) errors of Exp.Expon over time to delivery.
Energies 12 04518 g007
Figure 8. PB τ of selected best performing models over probabilities τ { 0.01 , 0.02 , , 0.99 } .
Figure 8. PB τ of selected best performing models over probabilities τ { 0.01 , 0.02 , , 0.99 } .
Energies 12 04518 g008
Figure 9. Results of the Diebold-Mariano test. (a) presents the p-values for the | | · | | 1 norm with the CRPS loss, (b) the values for the | | · | | 2 norm with the CRPS loss. The figures use a heat map to indicate the range of the p-values. The closer they are to zero (→ dark green), the more significant the difference is between forecasts of X-axis model (better) and forecasts of the Y-axis model (worse).
Figure 9. Results of the Diebold-Mariano test. (a) presents the p-values for the | | · | | 1 norm with the CRPS loss, (b) the values for the | | · | | 2 norm with the CRPS loss. The figures use a heat map to indicate the range of the p-values. The closer they are to zero (→ dark green), the more significant the difference is between forecasts of X-axis model (better) and forecasts of the Y-axis model (worse).
Energies 12 04518 g009
Figure 10. Average (a) location μ and (b) scale σ parameters of model GenGam.Quadr.Expon over products at time T i { 3.25 , 3 , , 0.5 } .
Figure 10. Average (a) location μ and (b) scale σ parameters of model GenGam.Quadr.Expon over products at time T i { 3.25 , 3 , , 0.5 } .
Energies 12 04518 g010
Table 1. Bias, MAE, RMSE and CRPS for the considered models on the interval [ T 1 , T 2 ) = [ 3.25 , 0.5 ) , see Equation (15). Bolded values indicate the lowest value in each column, instead of bias, where it indicates the value closest to 0.
Table 1. Bias, MAE, RMSE and CRPS for the considered models on the interval [ T 1 , T 2 ) = [ 3.25 , 0.5 ) , see Equation (15). Bolded values indicate the lowest value in each column, instead of bias, where it indicates the value closest to 0.
BiasMAERMSECRPS
Exp.Const144.113220.132155.464100.161
Exp.Lin1.703151.349114.71067.504
Exp.Quadr−28.784142.488108.48263.410
Exp.Expon−40.138138.866105.59261.761
Gamma.Const.Const156.004225.793158.639100.991
Gamma.Lin.Const−0.785150.973114.57366.631
Gamma.Lin.Lin24.151163.237123.18772.023
Gamma.Quadr.Const−28.934142.568108.52462.857
Gamma.Quadr.Lin−26.498143.301108.90862.413
Gamma.Quadr.Quadr−29.339146.516111.38564.006
Gamma.Quadr.Expon−22.970145.052110.12863.170
Gamma.Expon.Const−40.136138.952105.60761.218
Gamma.Expon.Lin−36.848138.829105.44860.427
Gamma.Expon.Quadr−40.102139.659105.95760.770
Gamma.Expon.Expon−41.907139.741106.09860.845
GenGam.Const.Const152.901252.987164.07290.424
GenGam.Lin.Const72.961189.489129.54568.444
GenGam.Lin.Lin−17.116160.708118.30860.203
GenGam.Quadr.Const67.173188.348129.67468.624
GenGam.Quadr.Lin−52.039161.245123.31262.396
GenGam.Quadr.Quadr−62.730165.170126.51564.582
GenGam.Quadr.Expon−63.908153.412118.43258.893
GenGam.Expon.Const69.276190.731130.24269.003
GenGam.Expon.Lin16.182174.413124.55664.196
GenGam.Expon.Quadr−59.219156.238120.49959.710
GenGam.Expon.Expon−0.735173.996125.89764.519
GenF.Const.Const77.112216.268146.78975.662
GenF.Lin.Const8.387170.659124.04361.562
GenF.Lin.Lin−54.312163.464126.61161.717
GenF.Quadr.Const4.748170.497124.34561.797
GenF.Quadr.Lin−101.139167.424134.29164.764
GenF.Quadr.Quadr−102.141172.586137.63068.198
GenF.Quadr.Expon−109.600164.070132.86463.636
GenF.Expon.Const8.630173.305125.66862.511
GenF.Expon.Lin−39.876168.606128.45562.641
GenF.Expon.Quadr−106.988164.796133.92463.842
GenF.Expon.Expon−41.044172.291130.45663.797

Share and Cite

MDPI and ACS Style

Narajewski, M.; Ziel, F. Estimation and Simulation of the Transaction Arrival Process in Intraday Electricity Markets. Energies 2019, 12, 4518. https://0-doi-org.brum.beds.ac.uk/10.3390/en12234518

AMA Style

Narajewski M, Ziel F. Estimation and Simulation of the Transaction Arrival Process in Intraday Electricity Markets. Energies. 2019; 12(23):4518. https://0-doi-org.brum.beds.ac.uk/10.3390/en12234518

Chicago/Turabian Style

Narajewski, Michał, and Florian Ziel. 2019. "Estimation and Simulation of the Transaction Arrival Process in Intraday Electricity Markets" Energies 12, no. 23: 4518. https://0-doi-org.brum.beds.ac.uk/10.3390/en12234518

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop