Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
71,855 characters · 16 sections · 45 citation commands
Probabilistic Forecasting in Day-Ahead Electricity Markets: Simulating Peak and Off-Peak Prices
In the last few decades since the deregulation of electricity markets it has become increasingly important to capture uncommon features of electricity prices such as nonstorability, which makes electricity prices really volatile ( see weron2014electricity). In this paper we use different time series models to forecast electricity by simulation and then evaluate those forecasts using various criteria with different properties. We believe it is crucial to take into account the dependency structures in order to properly forecast multivariate time series. The innovation in this paper is that we include dependency structures in some of the multivariate forecasting models and in one of the forecast evaluation criteria to show that the incorporation of the dependency structures substantially improves electricity price forecasts. The electricity prices we model and forecast in this paper are the peak and off-peak price series based on the German-Austrian day-ahead price. These time series are important for derivatives trading.
As mentioned above, electricity prices show special characteristics which are usually classified in the relevant literature ( see weron2014electricity and ziel2016forecasting). Specifically, these properties are i) mean reverting behavior; ii) seasonal behavior; iii) time dependent volatility; iv) price spikes; and v) cross-period effects (e.g. night hours influence day-time hours even though they take place on the same day). All these aspects are known in the literature, but there is no electricity price forecasting model which incorporates all of them. For instance, karakatsani2008forecasting cover all the above effects except interaction effects, ziel2015efficient consider all effects except price spike effects. We propose electricity price models which incorporate all the said effects into a probabilistic electricity price forecasting framework.
A two-step approach is used to forecast prices. In the first step the conditional mean model is estimated: the mean must be properly estimated so that the residuals have a zero mean. Therefore, in the mean equation all the seasonal properties must be included. Accordingly, uniejewski2016automated and ziel2018day propose mean equations with autoregressive, non-linear effects and seasonal effects. Once the conditional mean model is properly estimated we proceed to estimate the residuals, which must have a zero mean so the models differ in the structures of the standard deviation. We consider mean reverting jump diffusion models (MRJD) such as the model included by seifert2007modelling and ioannou2018effect applied to electricity prices. The MRJD model is an Ornstein–Uhlenbeck (OU) process proposed by uhlenbeck1930theory. Unlike weron2008market and cartea2005pricing, where the jump component is first estimated and then an OU process is assumed in the continuous part; we first estimate the mean model and then we assume a MRJD structure in the residuals. weron2014electricity offers a good review of MRJD models applied to electricity price forecasting. As mentioned above, our interest is the dependency structure between different time series, which we include by assuming correlated jump occurrence processes, a procedure we believe has never been used before. To obtain a correlated jump we focus on the bivariate Bernoulli process, proposed by dai2013multivariate.
Once the models are estimated electricity prices are simulated and forecast, then those forecasts and their paths are evaluated using different criteria. In this article we use four different criteria: mean absolute error (MAE), mean square error (MSE), pinball score (PB, also known as quantile loss) and energy score (ES). The first two are the most widely used in the literature of forecasting evaluation; for instance, keles2012comparison apply MSE to evaluate the electricity price forecasts from a model including spikes and other structures such as ARIMA and GARCH. voronin2014hybrid use the MSE and MAE criteria to evaluate the performance of different electricity price forecasts in the NORDPOOL market. In this paper we focus more on the PB and ES as we are interested in the performance capturing the whole distribution and how the different models capture dependency structures. The PB has been applied by maciejowska2016hybrid, dudek2016multilayer, and juban2016multiple, all involving an electricity price forecasting competition with the PB used to check performance, as the objective was to approximate the forecast distribution. The ES has not been applied to electricity price forecasts so far, weron2019electricity. However, it has been applied a few times in the energy forecasting context, e.g. in pinson2012evaluating for wind power forecasting. The ES is built up as per gneiting2007strictly and then applied to our time series. We pay more attention to this score because it takes into account dependency structures. As mentioned above, our contribution is to include correlation structures in the models as well as in the evaluation. Then, to check whether the differences between the forecasting performances of the models in pairs are significant or not, the Diebold-Mariano (DM) test is applied.
The rest of the paper is organized as follows: Section (ref) explains the data and highlights the relationships to the derivative markets, Section (ref) introduces the models, Section (ref) explains the estimation methods and how the forecasts are generated, Section (ref) describes the evaluation criteria, Section (ref) discusses the results, and Section (ref) summarizes our results, and outlines the most important facts.
As mentioned in the introduction, we focus on off-peak and peak price series from the EPEX market because they are relevant for derivative trading, especially future products. On the European Energy Exchange (EEX) different future products for electricity with cash settlement for the German/Austrian delivery zone are traded. They are base, off-peak and peak price products (also known as Phelix) traded at EEX. The underlying of these products are based on the hourly German/Austrian EPEX day-ahead electricity prices. The Phelix base product is simply calculated as the mean of all hourly EPEX prices in the delivery period. For example, the underlying of Phelix base week future contracts are calculated as the mean of the 168 hourly prices from Monday 0:00-1:00 to Sunday 23:00-24:00. For Phelix peak products the underlying is the mean of the day-ahead price from the $9^{\text{th}}$ hour of the day to the $20^{\text{th}}$ (12 hours in total) on Monday to Friday. Thus, for Phelix peak week futures contracts the underlying are computed as a mean of the $5\times12=60$ hourly mean prices for the peak hours from Monday to Friday. The remaining $168-60=108$ hours would be the underlying for Phelix off-peak week future products. However, Phelix off-peak products are only available for longer delivery periods (month, quarter and annual) and are rather illiquid. Therefore, the primarily focus for traders is on the Phelix base and peak products.
As traders focus on base and peak products it makes sense to concentrate on forecasting the corresponding underlyings. However, the fact that definition of the Phelix peak products depends on the day of the week makes the modeling a bit cumbersome. Intuitively, it makes sense to model and forecast the daily base price (the mean of the 24 hourly prices) and the daily peak price (the mean of the 12 prices 8:00-9:00 to 19:00-20:00). Of course, for trading Phelix peak products a forecast for Saturday and Sunday peak prices is not relevant. Nonetheless, it is more convenient to model the peak price in the above mentioned manner to preserve the time series structure. However, the base and peak time series are partially based on the same prices, in fact the peak prices. But, from the modeling perspective it is more convenient to have less correlated data. This linear dependency can be reduced easily by modeling the daily peak and off-peak prices as they are computed based on completely different hourly prices. If we are interested in a base price, we may obtain it directly by averaging the daily off-peak and peak prices. Hence, it is completely sufficient to model the base and peak prices for trading purposes. Thus, we proceed to analyze the above mentioned time series henceforth.
Finally, we would like to mention that it would be more informative to have a model for the 24 hourly electricity prices than just a model for the peak and off-peak prices. The problem with these models in the considered probabilistic forecasting setup are the computational burden, as there would be too many variables to estimate and we would not be able to optimize the models. However, these forecasts would not add any information regarding derivatives markets because hourly prices are not traded in these markets.
The considered electricity price data starts on 1 st January 2014 and ends on 31 st December 2017. It is measured in EUR/MWh. To calculate the higher moments and the dependencies, use the following notation;
$$ \text{m}_{i,j} = \mathbb{E}\left[\left(\frac{Y_{d,1}-\mu_{Y_{d,1}}}{\sigma_{Y_{d,1}}}\right)^i \left(\frac{Y_{d,2}-\mu_{Y_{d,2}}}{\sigma_{Y_{d,2}}}\right)^j\right], $$ where $Y_{d,1}$ and $Y_{d,2}$ refer to off-peak and peak time series with their means $\mu_{Y_{d,1}}$ and $\mu_{Y_{d,2}}$ and standard deviations $\sigma_{Y_{d,1}}$ and $\sigma_{Y_{d,2}}$. We show below the sample statistics (of the input data), but these may not be good estimators for the corresponding statistical counter-parts. However, under some mixing assumptions (e.g. weakly periodically stationary) the sample mean/variance/skewness/kurtosis/etc converge to the corresponding counterpart. Additionally, we would like to point out the fact that if the time series are bounded, and this is our case, then all moments exist. The sample descriptive statistics for both time series are shown below:
As expected, Table (ref) shows that the mean and the standard deviation are higher in the peak time series. As the volatility is higher the range for the peak series is higher than that of the off-peak time series. The correlation shows quite a high positive linear relationship between the two time series. The skewness shows that the off-peak series is clearly asymmetric and that the peak series is slightly asymmetric. The coskewness coefficients show how the variance of one time series and the mean of the other are related. As observed in Table (ref), the relationship between the off-peak central variance and the peak central mean is stronger than the other way round; in the case of m$_{2,1} = -0.59$, this means that the higher the value of the peak series the lower the variance of the off-peak series. In view of these results it can be concluded that none of the time series follows a normal distribution pattern.
Figure (ref) shows the off-peak and peak time series through our sample. The first two years are used only for estimation purposes and the last two years are first predicted and then used as observations of the following rolling windows. How the rolling windows are developed is explained in Section (ref). Figure (ref) is divided in two to emphasize this aspect. It is observed in Figure (ref) that the volatility was higher at the beginning of 2017 and also at the end of the year. As can be observed in Table (ref) and in Figure (ref), the volatility is higher, and so is the mean in the peak series compared to the off-peak figures. However, generally the trend in the graphs is quite similar, as shown by the correlation coefficient. In both cases there is evidence of volatility clustering and spikes.
The histograms and density functions in Figure (ref) show the distribution of the two time series. It may be observed that both series have heavy tails and the asymmetry is more pronounced in the off-peak series. In both cases there is evidence of spikes, which are rare events where the price is extremely low or high. Regarding the scatter plot, a strong and complex correlation between the two is confirmed, which leads us to include correlation structures in our models. In the models we propose, the correlation is not included only in the continuous part of the variation but also in the jump occurrence process, as the depicted graphs show. From the scatter plot it is also possible to observe the bivariate density, which shows how the scatter plot is distributed. In Figure (ref), the darker colors show the higher quantiles of the distribution. The bivariate distribution confirms the intuition of the scatter plot, where the darker areas are those where there are more points. In both graphs - the scatter plot and the bivariate distribution - one may observe that the spread is higher in lower values than in higher values.. The correlation coefficient for the values when the off-peak price is lower than 30€ is 0.65, while when the off-peak price is higher than 30 the correlation coefficient is 0.84. This is an example of the complicated dependency structure.
The models that we analyze in this paper are two step models. In the first step we estimate the mean equation and in the second we study the residuals from the previous step.
For the sake of simplification we define $\bm{Y}_d = (Y_{d,1}, Y_{d,2})'$ as the bivariate vector of the off-peak and peak prices, so index $1$ corresponds to the off-peak price and index $2$ to the peak price.
In this subsection we introduce the conditional mean model that we assume. The mean equation is based on the mean models proposed in uniejewski2016automated and ziel2018day. To calculate the mean equation we assume a model with autoregressive structure with exogenous variables (ARX) for the peak and off-peak series. The ARX model was shown to perform really well in forecasting electricity prices in uniejewski2016automated and ziel2018day. We consider the mean model for the two time series as:
where $i=1,2$ and $\text{DoW}_d^j$ is a day of the week dummy of day $j$ at day $d$ such that e.g. $\text{DoW}_d^1$ is $1$ if $d$ falls on a Monday, $\text{DoW}_d^2 = 1$ if $d$ is on Tuesday etc. The residuals are $\epsilon_{d,1}$ and $\epsilon_{d,2}$, and by construction the mean of the two terms must be 0. The model has in total $p=1 + 2\times 8 + 3\times7 = 38 $ parameters with corresponding parameter vector $ \bm{\beta}$. Obviously, model (ref) is a linear model that can be written as
where $\bm{X}_{d,i} $ and $ \bm{\beta}_i$ are $p$-dimensional.
The error terms are considered to be distributed as:
where $\bm{\epsilon}_d=(\epsilon_{d,1},\epsilon_{d,2})'$, $\bm{0} = (0,0)'$ and $\bm{\Sigma}$ is the covariance matrix of $\bm{\epsilon}_d$.
Model (ref) covers the major characteristics of electricity prices, especially mean reverting properties, seasonal structure, and cross-period effects. Only volatility and price spikes are not captured by the structure assumed. Hence, for all the remaining models we consider the same mean equation, but modify the error model (ref) to capture the missing effects.
In this subsection we explain the ARX-IJ model. We consider MRJD in each residual independently. This is the standard OU process applied in electricity price forecasting, and has been applied several times, e.g. in keles2012comparison and widely discussed in weron2014electricity. Jump diffusion models are accurate for capturing price spikes as observed in the tails of Figures (ref) and (ref).
After Euler discretization the model is written as follows:
where $\bm{\epsilon}_{d,cont}$ is the continuous part of the error term and $\bm{\epsilon}_{d,jump}$ is the jump component. $\lambda_i$ is the probability of jumps, $\text{Ber}$ is the Bernoulli distribution, $\mu_i$ is the mean size of the jump, $\gamma_i$ is the standard deviation of the jump and $\sigma_i$ is the standard deviation of the continuous part all of them defined for $i=1,2$. We need the terms $-\lambda_{1} \mu_{1}$ and $-\lambda_2 \mu_2$ in the continuous term in order to ensure that the mean of $\epsilon_{d,1}$ and $\epsilon_{d,2}$ is 0, as it must be by construction. In this model we assume that the Bernoulli random variables $b_{d,1}$ and $b_{d,2}$ are independent. The conditional error term $\bm{\epsilon}_d | \bm{B}_d = \bm{\epsilon}_{d,cont} + \bm{B}_d \bm{\epsilon}_{d,jump}|\bm{B}_d$ is distributed as follows:
For the unconditional distribution of $\bm{\epsilon}_d$ first note that with $\text{Var}[XY] = \mathbb{E}[X]^2\text{Var}[Y] + \text{Var}[X] \mathbb{E}[Y]^2 + \text{Var}[X]\text{Var}[Y]$ the following holds:
Thus, as a result of the independence of all occurring random variables, it holds that
as $\text{Cov}[b_{d,1} \epsilon_{d,jump,1}, b_{d,2} \epsilon_{d,jump,2} ] =0$.\\
However, it is clear that $\bm{\epsilon}_d$ does not follow a bivariate normal distribution pattern.
The next model that we introduce, the ARX-BiJ model, is related to the previous one as it is based on an MRJD structure, but in this case the jump component is assumed to be bivariate ( more precisely bivariate Bernoulli). This dependency structure in the jump is one of our contributions to the literature. We further assume that the jump sizes can be correlated. We write the model as
As mentioned above, we now assume that the arrivals of the jumps are bivariate $\text{\textbf{Ber}}_2(\bm{P})$ distributed with probabilities $\bm{P}$. In this case, $p_{0,0}$ is the probability of no jump, $p_{1,0}$ is the probability of a jump occurring only in the off-peak component, $p_{1,1}$ is the probability of there being a jump in both components at the same time, and $p_{0,1}$ is the probability of a jump occurring only in the peak component. Therefore, $ \lambda_2 = p_{0,1}+p_{1,1}$ is the total probability of jumps in the peak component, and $\lambda_1= p_{1,0}+p_{1,1}$ is the equivalent probability in the off-peak series. The condition $p_{1,0}+p_{1,1}+p_{0,1} + p_{0,0}=1$ must hold. Unlike the previous model, in this case $b_{d,1}$ and $b_{d,2}$ are not independent: they must coincide in no jump with probability $p_{0,0}$ and in jump with probability $p_{1,1}$. For the bivariate Bernoulli setting we follow dai2013multivariate. Assuming a bivariate jump process. we capture dependency structure in the continuous component as well as in the jump component.
This subsection presents the ARX-BiJ-$\mu_d$ model. This model is very similar to the previous one, but in this case the mean of the jump is assumed to depend on the price observed previously. In order not to make things too tedious for the reader we only note those points that differ from the previous model, i.e.:
In this model we seek to capture the effect of the previous price on the mean of the jump component.
The next model (ARX-GARCH) considers bivariate constant conditional correlation GARCH (CCC-GARCH) structures, as first introduced by bollerslev1990modelling. We follow silvennoinen2009multivariate in their implementation:
where $z_{d,1}$ and $z_{d,2}$ are independent white noises with a standard deviation of 1. The parameters of the GARCH structure must fulfill $\alpha_{0,i},\alpha_{1,i},\alpha_{2,i}>0$ and $\alpha_{1,i}+\alpha_{2,i}<1$ conditions in order for the time series to be stationary. In this model we assume that the correlation between the two time series is constant and there are no cross-dependencies between the volatility series. Structures of this type are often used in the literature to forecast multivariate time series, for instance zanotti2010hedging and higgs2009modelling apply CCC-GARCH models in electricity markets.
Our last model (ARX-BiJ-$\mu_d$-GARCH) includes CCC-GARCH structures in the continuous component of the model described in Equation (ref). This is our most complex model, and it is a combination of the ARX-BiJ-$\mu_d$ and ARX-GARCH models:
where all the components are assumed to be distributed as in the previous subsections. As it is the most complex model, it has the largest number of parameters to estimate. The model is able to capture all the aspects mentioned above.
For the estimation we assume that there are $D$ observations available. We denote the resulting price vectors and regression matrix by $\mathbb{Y}_i = (Y_{1,i}, , \ldots, Y_{D,i})'$ and $\mathbb{X}_i = (\bm{X}_{1,i}', , \ldots, \bm{X}_{D,i}')'$, corresponding to regression equation (ref).
To estimate the ARX model (Equation (ref)) we apply two different estimation methods: OLS\footnote{In order to avoid perfect collinearity in the OLS estimation we drop the interaction between Wednesday and the previous observations for both time series.} and elastic net. This gives us two different estimations and therefore two different forecasts, which we note as ARX-OLS and as ARX-enet, respectively.
Using the OLS estimator the estimated values are: $$ \widehat{\bm{\beta}}^{\text{OLS}}_i = \operatorname*{arg\,min}_{\bm{\beta}\in\mathbb{R}^p} \left[ \| \mathbb{Y}_i - \mathbb{X}_i'\bm{\beta} \|_2^2\right], $$
The second estimation method applied to estimate Equation (ref) is the elastic net, introduced by zou2005regularization, which is very similar to OLS but has quadratic and linear penalties. However, in defining the elastic net estimator it is crucial to consider the corresponding scaled OLS problem. Hence, we introduce $\widetilde{\mathbb{Y}}_i $ and $ \widetilde{\mathbb{X}}_i$ as a scaled response vector and scaled regression matrix. We require them to be scaled in such a way that any column has a zero mean and standard deviation of 1.
Given the scaled OLS problem, the scaled elastic net estimator is given by the optimization problem $$ \widehat{\widetilde{\bm{\beta}} }_i^{\text{enet}} = \operatorname*{arg\,min}_{\bm{\beta}_\in\mathbb{R}^p} \left[ \| \widetilde{\mathbb{Y}}_i - \widetilde{\mathbb{X}}_i'\bm{\beta}||_2^2 + \lambda\left(\frac{1-\alpha}{2} ||\bm{\beta}||_2^2 + \alpha||\bm{\beta}||_1\right)\right], $$ where $\lambda$ and $\alpha$ are tuning parameters that characterize the penalty term $\lambda\left(\frac{1-\alpha}{2} ||\bm{\beta}||_2^2 + \alpha||\bm{\beta}||_1\right)$. We receive the (unscaled) elastic net estimator $\widehat{\bm{\beta}}_i^{\text{enet}}$ simply by rescaling $\widehat{\widetilde{\bm{\beta}}}_i^{\text{enet}}$. If $\alpha=1$ the estimation method is equivalent to the lasso penalty developed by tibshirani1996regression, and when $\alpha=0$ it is equivalent to the ridge penalty first introduced by hoerl1970ridge. The lasso estimator has the property of sparsity, which means that for certain values of $\lambda$ the resulting solution sets irrelevant parameters to zero while keeping relevant parameters at non-zero. The lasso estimation enjoys some popularity in electricity price forecasting: see ziel2016forecasting, gaillard2016additive,steinert2019short, narajewski2019econometric, amongst others.
The elastic net can be seen as an augmented data lasso shrinkage with some ridge elements. Like the lasso, the elastic net has automatic sparsity property for $\lambda > 0$. In the lasso estimation we apply cross validation (CV) to select the tuning parameter $\lambda$. It takes into account the number of observations, the number of parameters, the variance and the correlation, according to hebiri2013correlations. However, uniejewski2016automated conclude that better forecasts are achieved when the elastic net method is consider, thus we incorporate a ridge penalty to the lasso estimation method. uniejewski2016automated apply elastic net and lasso estimation methods (amongst others) in the electricity price forecasting context. Their findings suggest that $\alpha = 0.5$ is a good choice for applications, so we apply it in this paper as well. We choose $\lambda$ by 10-fold block-CV. In order to control weekly seasonal dependency structure we consider block-CV, the block length is 7 so that weekly seasonality is taken into account. This estimation method is often used in time series analysis, see. e.g. racine2000consistent.
In Table (ref) we show the percentage of times each one of the estimated parameters in Equation (ref) is not equal to 0, i.e. the percentage of times parameters are included in the model. As can be observed, the autoregressive components are relevant either for the peak series or for the off-peak series. The day of the week dummies (mainly weekend effects) are relevant as well. Regarding the iteration between the day of the week and previous observation, the number of times estimated parameters are included in the models decreases.
We compare the results for the two estimation techniques. The elastic net provides better forecasting results, so in estimating the second step the residuals are calculated via the elastic net.
In Figure (ref), the scatter plot of the residuals for two different days are shown, these are the residuals from the first of the rolling windows and the last one. We would like to point out firstly that the two scatter plots are very much alike which means that the ARX model performs similarly in both cases. Furthermore, the complex correlation structure is also observed, and is similar to the one in Figure (ref). The density of the scatter plot may be observed since the darker the color, the greater the number of points clustered.
In the second step residuals are estimated by maximum likelihood\footnote{In the ARX-IJ model, $\rho$ is estimated in a third stage as $\rho=\text{Cor}(\epsilon_{d,1},\epsilon_{d,2})$.} using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. This algorithm is a quasi-Newton method for non-linear optimization. In order to apply the maximum likelihood estimation we need a log-likelihood function which changes depending on the structure of the residuals assumed. That log-likelihood function is based on the assumed distribution for each of the different models explained in Section (ref). The maximum likelihood estimation is initialized as follows:
The above procedure to calculate the initial values ensures that whenever we include a new component, if the estimated parameters of this component are not zero, then the component improves the log likelihood function performance.
Once all the parameters are estimated, off-peak and peak time series are simulated. In our case we predict $H=7$ horizons and for each of horizon $M=16000$ paths are simulated. All the simulations (sometimes called ensemble) can be seen as multivariate probabilistic forecasts as they approximate well the underlying distribution of the forecasts. All relevant properties can be derived from these paths. It is possible to analyze only the marginal properties of each predicted horizon. The estimation and simulation process is repeated $N=731$ times, as mentioned in Section (ref), via a rolling window. The first estimation is made using the first two years of the data, then the next $H=7$ days are predicted with $M$ paths in each horizon. Then, the estimation sample is shifted one day forward and the process is repeated.
In this section we introduce the evaluation criteria used to assess probabilistic forecasting. In total we use four different criteria: MAE, MSE, PB, and ES. We first explain the scores that we use, then briefly introduce the DM test used to test whether differences in forecasting performance are significant or not, taking the models in pairs. To compute the DM test it is necessary to define a loss function. We therefore introduce each criterion with the corresponding loss function. As mentioned in Section (ref) we forecast seven horizons, and all the four criteria are independently evaluated for each of the horizons.
For evaluating the point forecasts, we consider the popular MAE and MSE measures. The MAE is a strictly proper forecasting criterion for the median and the MSE a strictly proper evaluation criterion for the mean, with "strictly proper" here meaning that only the perfect model minimizes the corresponding criterion. Therefore we define $\widehat{Y}_{d,i}^{\text{med}}$ as a median forecast and $\widehat{Y}_{d,i}$ mean for day $d$ and volatility series $i$ derived from the sample counterparts of the $M$ simulated paths. The MAE and MSE are defined using the absolute error (AE) and the squared error (SE). Thus, with
equation (ref) is often used in the literature for mean forecasts but it is not proper from a statistical point of view. It is proper when we have symmetry in the sample as in this case median is equal to mean. We define,
for $ i=1,2$. Hence, we can evaluate the point forecasts for the off-peak and peak price separately.
The two criteria introduced above are the most widely used in the literature, but we are more interested in the marginal properties of the models, and so we use the PB.
The PB measures the distance for each quantile. Therefore, as for the median and mean forecasts, we define $\widehat{Y}_{d,q,i}$ as a forecast for the $q$-quantile day $d$ and time series $i$. We get these quantile forecasts by taking the sample quantile of our $M$ simulated paths. The pinball loss function is computed as follows,
where $q\in \bm{Q}$ is a quantile, in our case $\bm{Q}= \{Q_q\}_{q\in \{1,\ldots, K\}} = \{ 0.01,0.02,\ldots,0.99 \}$ with $K=99$. $\widehat{Y}_{d,q,i}$ stands for the estimated price of quantile $q$ on day $d$ and in time series $i$ and $Y_{d,i}$ is the observed value at day $d$ and time series $i$. To calculate the PB of quantile $q$, series $i$ and N days, we proceed as follows: $$ PB_{q,i} =\frac{1}{N}\sum_{d=1}^N PB_{d,q,i} \ \ \text{for} \ \ i=1,2 \ \ \text{and} \ \ q\in Q, $$ Thus, the PB for $N$ days is computed by averaging across the $\bm{Q}$ quantiles, $$ PB_{i} =\frac{1}{K}\sum_{q=1}^K PB_{Q_q,i} \ \ \text{for} \ \ i=1,2 \ \ \text{and} \ \ q\in Q, $$ where $K$ is the number of quantiles. Note that if the distance in the quantile grid $\bm{Q}$ converges to $0$ then the PB converges to the probabilistic forecasting evaluation measure CRPS (continuous ranked probability score), which is strictly proper with respect to the (marginal) distribution of $Y_i$. For further information on the PB score, see steinwart2011estimating.
When applying PB we may observe all the marginal properties of the forecasts, given that we may observe the forecasting performance of the different models in each quantile. Therefore, the first criterion introduced in this section is merely a special case of PB when $q=0.5$. With this score we can compare how the different models capture spikes, as we can analyze the behavior in the tails.
The last criterion that we use is the ES, which is a generalization of the CRPS. The ES is the only score that takes into account the dependency structures. This score is applied to all the variables at the same time in order to take in the correlation. In both cases, in the modeling and evaluating we pay close attention to the dependency structure as this is the main contribution of our paper. The loss function of the ES is computed as follows:
$\bm{Y}_d^{[m]}$ for $m=1,\ldots,M$ is the predicted $m^{th}$ path of the multivariate data for day $d$. Our estimator for the ES is based on gneiting2007strictly and is divided in two parts: the ED and EI. The ED measures the mean Euclidean distance from each one of the paths to the observed value, so this part is measuring marginal properties. On the other hand, the EI measures within the path dependencies, i.e. how well the path are aligned with each other. The EI appears with a negative sign which means we are interested in spread amongst trajectories. Equation (ref) is the ES for day $d$, thus the ES for N days is calculated as, $$ \text{ES} = \frac{1}{N} \sum_{d=1}^N \text{ES}_d. $$
The next step is to check whether or not the differences between the forecasting performances with each criterion are significantly different from zero. For the significance test we use the DM test, which compares models in pairs. As mentioned above, we need a loss function such as the ones written above to apply the DM test. Let $L_d$ denote the loss function of a certain model; the loss differential between models $\mathbb{A}$ and $\mathbb{B}$ is defined as $\delta_{d,\mathbb{A},\mathbb{B}}=L_{d,\mathbb{A}}-L_{d,\mathbb{B}}$. The only required assumption is for the loss differential to be covariance stationary. To apply the DM test we compute
$$ \frac{\widetilde{\delta}_{\mathbb{A},\mathbb{B}}}{\sigma_{\widetilde{\delta}_{\mathbb{A},\mathbb{B}}}} \sim \mathcal{N}_1 (0,1) $$
where $\widetilde{\delta}_{\mathbb{A},\mathbb{B}}=\frac{1}{N}\sum_{d=1}^N \delta_{d,\mathbb{A},\mathbb{B}}$ and $\sigma_{\widetilde{\delta}_{\mathbb{A},\mathbb{B}}}$ is the standard error which we estimate by the corresponding sample counterpart. For further information on the DM test, see diebold2002comparing and diebold2015comparing.
In this section we assess the forecasting performance of each model, using the criteria introduced in Section (ref). Then, the DM test is applied to check whether the differences between the models in pairs are significant or not\footnote{DM test results for all criteria and all horizons are available upon request.} for each criterion. The forecasting horizon (H) is 7, which means that for each rolling window we get the forecasts for the following 7 days. In order to show what the different trajectories look like, we present trajectories for the ARX-enet and ARX-BiJ-$\mu_d$ models at the end of the Section.
In the following graphs we show how each one of the models performs using all the four criteria from Section (ref). According to all the evaluation criteria the optimal scores zero, so the lower the scored value the better the forecasting performance. In all four criteria, the relative performances compared to the ARX-OLS model are shown, i.e. the performance of each one of the models is divided by the ARX-OLS values for each one of the horizons Score(model)/Score(ARX-OLS).
Figure (ref) shows the results for the MAE criterion for the two time series. As explained above, the MAE criterion is a special case of the PB that compares performance on the median. As observed in Figure (ref), the differences between the models are not too big: the only clear result is that the performance of the ARX-OLS model is significantly poorer. Overall, results for the ARX-GARCH and ARX-IJ can be said to be slightly better but, the forecasts are not significantly better according to the DM test. It is important to underline that a comparison between the OLS and elastic net estimation methods reveals that the latter gives significantly better forecasting results.
Figure (ref) shows the forecasting performance at the mean. A look at Figure (ref) and the DM test suggests that bivariate jump structures are not effective in capturing mean behavior, as the forecasting performance is poor according to the MSE criterion. At the same time, the Figure shows that the forecasts of the ARX-IJ are better than other models and the DM test confirms that the differences are significant, with the exception of the first horizon, where the difference between ARX-IJ and ARX-enet is not significant for either of the time series. The latter means that for the mean forecast, it is important to introduce jump structures; but dependency structures are not highly relevant. The superiority of the ARX-IJ forecasting performance becomes greater when the horizons are increased. According to this criterion, elastic net forecasts are significantly better than the OLS forecasts when a simple ARX model is simulated.
As shown in the previous section, the PB takes into account the whole distribution of the forecast paths quantile by quantile. As shown in Figure (ref), in the case of the peak time series the ARX-BiJ-$\mu_{d}$ models forecast outperforms the other models except in one case ( in H2 the PB of the ARX-BiJ model is lower). On the other hand, in the off-peak seriesfor the first horizons, the ARX-BiJ forecasts are seen to be the best, but after the $6^{\text{th}}$ horizon the ARX-enet has the best forecasting performance. This is curious because the ARX-enet does not take into account heavy tails and a model with jumps would be expected to capture tail behavior more efficiently. However, it must be underlined that the differences between these models are not significant according to the DM test. The DM test only concludes that in both time series the forecasting performance is significantly poorer for the ARX-OLS and the ARX-BiJ-$\mu_{d}$-GARCH models. As with the MAE and MSE criteria, with the PB score the elastic net estimation method provides better forecasts than the OLS .
Figure (ref) shows $PB(model)-PB(OLS)$ for the first forecast horizon quantile by quantile, that is all the models are compared to the ARX-OLS. Focusing on the model ARX-BiJ-$\mu_{d}$-GARCH, which is the most complex model, it may be observed that it is the best model in the first and the last quantiles, but in the middle quantiles it has the worst forecasting performance. This is even more pronounced for the peak time series. The behavior of the forecasts with ARX-BiJ-$\mu_{d}$ and ARX-BiJ models is similar in both time series. In the off-peak time series, we observe that in the middle quantiles the forecasting performance is the best for the ARX-enet and the ARX-GARCH. Furthermore, it is surprising that the ARX-IJ performs so poorly in the highest and lowest quantiles because unlike ARX-OLS, ARX-enet and ARX-GARCH models, the ARX-IJ model considers spikes which captures behavior on the tails.
So far we have distinguished between the two time series because dependencies are not taken into account in the criteria mentioned above. It may be observed for the three previous criteria that the off-peak time series has a lower error term, which means that forecasts are more accurate according to all three criteria. This fact is expected as the off-peak time series is less volatile and therefore easier to predict. The ES takes into account dependency structures, which are the key feature of this paper.
We start analyzing the ES performance by showing the two parts defined on Equation (ref). The ED measures the euclidean distance, and only takes into account marginal properties. This measure is similar to the MSE measure but taking the two time series at the same time. Thus, as in Figure (ref) there are no big differences between the two times series performance, we are expecting a similar performance with the ED measure. The EI measures within the path dependency, i.e. how well the trajectories are aligned with each other.
Figure (ref) shows the forecasting performance for the two parts of the ES. According to the ES criterion a model has a better forecasting performance when the estimated value of the score is lower. Thus, on the ED lower values are needed while on the EI higher ones are better. Focusing on the ED, poorer performance of the models with bivariate jumps is observed (similar to the MSE criterion), and the ARX-IJ model performs the best. Regarding the EI, it is shown that including spikes in the models increases the spread of the trajectories considerably, thus increasing the estimated value of the EI. Moreover, incorporating bivariate jumps increases the spread even more. The differences on the ED are lower than the differences on the EI, and even if the EI value is multiplied by 0.5 the difference is still bigger. Consequently, models containing bivariate jumps perform better on the ES.
Figure (ref) shows the ES for our seven models and our seven horizons. In Figure (ref) the relative performance against the ARX-OLS model is shown and in Figure (ref) the shape of the ARX-OLS model is depicted. It is observed that the best forecasts are made by the ARX-BiJ-$\mu_d$, which makes sense given that it incorporates a complex dependency structure. One would expect the ARX-BiJ-$\mu_d$-GARCH to capture all features of the time series better, but it may the case that there are too many parameters to estimate: with our starting values we have hit a local maximum\footnote{We have tried different starting values but the ones used in the estimation have the best results. Other starting values could be used but this would greatly increase a lot the computation time. As 731 estimations are made in total, increasing the number of starting values would make estimation unfeasible.}. Another reason might be structural breaks in the dataset, leading to poor forecasting performance. As already mentioned in the previous paragraph, in Figure (ref) is shown that that the models with bivariate jumps generally provide better forecasts with this criterion, which means that the dependency structures are properly included by considering bivariate jumps. This is an important result as our objective is to efficiently capture the complex correlation structure between our price series.
As may be observed in Table (ref), with ES it is not possible to distinguish between the OLS and enet estimation methods, while with the other three criteria the elastic net estimation method procures significantly better forecasting performance. This happens due to the fact that the performance with the ED is better with the ARX-enet model while the opposite occurs with the EI. Besides, it is shown that the ARX-BiJ-$\mu_{d}$ model provides significantly better forecasts, followed by the ARX-BiJ model. As already noted, according to the ES criterion models with bivariate jumps offer significantly better results, which means that our adjustment on the jump diffusion models helps to catch the dependencies efficiently. It is clear that according to this score the assumption of no constant mean of the jump helps with forecasting accuracy. The improvement in incorporating no constant mean to the ARX-BiJ model was not significant when applying the other three criteria. The forecasting performance of the ARX-GARCH model is weaker than expected. This could be because CCC-GARCH structures are more focused on symmetric effects.
In figures (ref) and (ref) we show the first 100 trajectories\footnote{It is not possible to see all the 16000 different trajectories at the same time.} for days 2017-12-12 and 2017-12-16; these days refer to the horizon 0, i.e. when the models are run and the following 7 days are predicted. Only paths for the ARX-enet and ARX-BiJ-$\mu_d$ models are depicted\footnote{It would be tedious for the reader to show paths for all the seven models.}. It may be observed that the trajectories are random and that the observed trajectory or true trajectory is located between the paths most of the times. The shape of the trajectories with the ARX-enet model is more compacted and extreme events are not well captured, some real values are out of the 100 trajectories. On the other hand, it may be observed that the with the ARX-BiJ-$\mu_d$ model the spread of the trajectories is bigger while the trajectories in the middle are more dense. Remember some of the observed values are out of the range of 100 bivariate trajectories of the ARX-enet model, e.g. Figure (ref). It may be observed in both figures that the correlation between peak and off-peak time series is high as both moves in the same direction. However, in practice it might be highly relevant if price spikes occur at peak and off-peak prices together or not. The ES is the only considered criterion that takes into account the full distribution of the bivariate data. Thus, it is the only criterion which can discriminate for the performance of such a double spike event. It can analyze the whole distribution and it is possible to have a full picture of the predictions. According to the ES, the ARX-BiJ-$\mu_d$ outperforms the ARX-enet model.
The relevance of the ES can be highlighted from the practitioners point of view as well. Imagine there is a trader who manages a portfolio with two assets, peak and off-peak futures contracts. If the portfolio manager is interested in the mean of the portfolio, then the portfolio manager could focus on the MSE criterion. However, the mean behavior is also included in the ES criterion as full distribution is taken into account. On the other hand, portfolio managers are often interested in the value at risk (VaR) with respect to a certain probability of the portfolio. Then, only the ES is suitable for such an evaluation purpose, as the portfolio return is a weighted sum of two dependent random variables which depends on the full bivariate distribution. Note that if the peak and off-peak prices would be independent from each other then the resulting distribution could be derived by convolution, and the evaluation by the pinball score on a dense grid would be sufficient. Unfortunately, the prices are highly dependent. Therefore, practitioners should focus on the ES to make decisions regarding the VaR of the portfolio.
We would like to underline that none of the models used in this research is perfect. The ARX-BiJ-$\mu_d$ model that we claim to be the best one according to the ES criterion is not perfect as the model can be improved with respect MSE results. If one of the models would be superior to the other ones should be the best model according to all the criteria.
Proper modeling and forecasting in electricity markets is crucial for all participants. In this paper we focus on off-peak and peak time series. These time series are important for trading in derivatives markets, hence in this paper we try to improve already existing models to achieve accurate forecasts. Participants in derivatives markets can adjust their trading positions and evaluate trading strategies properly when there are accurate forecasts available.
The objective of this paper is to properly incorporate and evaluate the complex dependency structure in bivariate analysis. We believe that it is highly relevant for the forecasts to preserve the correlation structure from the original time series. In the literature, so far, MRJD models have been applied to time series, assuming independence between them. Our approach is to include bivariate jump occurrences in the MRJD model. We then need to assess whether these correlation structures have been properly included or not. To that end, we need a criterion which takes dependencies into account. In our case we use the energy score, but this is not the only criterion we apply: we also use the MAE, MSE, and the pinball score. Additionally, we apply the DM test to compare the models in pairs for each horizon.
It may be observed that models with bivariate jumps do not forecast better according to the MAE and MSE criteria. However, with the pinball score criterion, where the distribution of the forecasts is assessed, the performance is slightly better when a bivariate jump structure is considered in the model. Focusing on the ES score it is shown that models containing correlated jumps perform significantly better compared to those models without them. Nonetheless, the most complex model, which features bivariate jumps, no constant jump size, and CCC-GARCH structure, does not outperform the forecasts of the same model without CCC-GARCH. The way we have chosen the initial values for the ARX-BiJ-$\mu_d$-GARCH model ensures that the in-sample performance is better when CCC-GARCH structures are incorporated. However, the out-of-sample performance is poorer than expected.
For further research it might be interesting to develop dependency structure models considering hourly data and conduct a multivariate analysis with 24 variables. The problem with 24 variables is that the number of parameters to be estimated increases substantially. Here, vine copulas might help to tackle the above mentioned problem.
Peru Muniain is grateful for financial support from Dpto. de Educaci\'{o}n, Universidades e Investigaci\'{o}n del Gobierno Vasco under research grant IT-783-13, from Dpto. de Educaci\'{o}n, Pol\'{\i}tica Ling \"{u}\'{\i}stica y Cultura del Gobierno Vasco through Beca Predoctoral de Formaci\'{o}n de Personal Investigador no Doctor and research grant EGONLABUR from the same department. Peru Muniain also acknowledges financial support from the Spanish Ministry of Economics and Competitiveness (ECO2015-64467-R MINECO/FEDER).