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.
89,114 characters · 26 sections · 11 citation commands
Estimating Dynamic Conditional Spread Densities to Optimise Daily Storage Trading of Electricity
Whilst day ahead electricity price forecasting has been a topic of substantial and wide ranging research in terms of methods, the focus has mostly been upon price levels for the delivery periods (usually hourly) in the following day. More recently there has been an interest in density forecasts for the hourly prices, motivated by considerations of risk management, see weron2014electricity, nowotarski2017recent for extensive reviews. In this paper we provide a new formulation with a focus upon price spreads, and specifically we forecast the density functions for the intraday spreads in the day-ahead prices. The optimal operation of storage facilities, e.g. batteries, or load shifting programmes, e.g. demand-side management, over daily cycles depends upon these spreads if they are operated as merchants, arbitraging buying and selling from the wholesale market. Furthermore if risk is a consideration, analysis of the mean-differences in price levels would be inadequate, and we therefore directly estimate the density functions of all hourly spreads in prices at the day-ahead stage. Our specification, estimation and forecasting of these arbitrage spreads is new and computationally-intensive. We then show how these spread densities can support the optimal daily operation of a risk-constrained merchant battery facility.
In contrast to the body of work on gas and other storable commodities, e.g. boogert2011gas, secomandi2018improved, because of the daily periodicity in electricity prices and the predominance of the day-ahead auctions, which typically set all hourly prices for the following day simultaneously, the operational horizon for storage and load shifting is usually episodic on a daily basis. So the continuous time, dynamic optimisation formulations used for gas and other commodity storage operations are less appropriate for electricity, and furthermore the stochastic simplifications generally required in the analysis of these other models would not meet the requirements of adequate fit to the more complex power price dynamics. Thus, based upon day-ahead forecasts for the drivers of electricity prices, such as demand, wind and solar production, gas and coal prices, forecasts for electricity price levels have been proposed from various methods, e.g. nowotarski2015computing, garcia2012forecasting, karakatsani2010fundamental but, apparently, no methods have been developed specifically for forecasting intraday spread densities. Until recently, storage assets, such as pumped hydro storage would regularly store energy at night and discharge at the daily peak demand periods, which were quite predictable. But with the penetration of wind and especially solar generating facilities, the peaks and troughs in prices move around the day and in sunny locations with substantial solar energy, e.g. California, the lowest prices may often be in the middle of the day denholm2015overgeneration. Thus, the expected daily spreads in prices, and their consequent arbitrage opportunities for storage or load shifting, will depend upon wind and solar forecasts, as well as demand and supply considerations. Furthermore the price density functions are non-normal with skewness switching between positive and negative depending upon the dynamics of production of renewable energy gianfreda2017stochastic.
We apply our spread densities formulation to the German market. This is the largest and the main daily reference market for wholesale power in Europe. It is also strongly influenced by wind and solar production, as well as providing a context where batteries and demand-side management are active innovations. The day-ahead auction has been actively researched and closes at noon each day, with the vector of 24 hourly prices for the next day being released an hour later. In this context, a reasonable question might be why is there a need in the morning to forecast these prices, if they are available by 13:00, and apply to the following day? For the operation of spread-based arbitrage facilities in particular, it is clear that the traders need to have a forward plan in order to decide whether they will be entering the auction to buy or to sell at particular hours and thereby formulate their bids and offers into the auction accordingly.
For modelling the spread densities, we adapt the Generalised Additive Model for Location, Scale and Shape (GAMLSS) parametric regression model stasinopoulos2007generalized, which has already been used effectively to form day-ahead densities of price levels in the German context gianfreda2017stochastic. Within this framework, the hourly electricity price spreads form a response variable, whose distribution function varies according to multiple exogenous factors. The GAMLSS framework allows choice from a wide range of distributions, whose moments change according to the exogenous variables specified using (non)linear relationships. The dynamic location, scale and shape parameters ("latent moments" related to the mean, volatility, skewness and kurtosis of price spreads) are therefore explicitly incorporated into the forecasting model.
The paper therefore proceeds by first describing the data and the density estimation process. In section 3 we use the Pinball Loss function to select the best fitting density model with four latent moments. Then in section 4 we undertake a rolling window forecasting evaluation and demonstrate the value of the dynamic, conditional latent moment estimates. Section 5 uses these density functions to devise the optimal daily scheduling of battery storage and demonstrates through backtesting the superior profitability compared to using normal densities. Section 6 concludes.
The German hourly day-ahead electricity price, wind forecast, solar forecast and actual total load data were downloaded from the Open Power System Data website \url{https://data.open-power-system-data.org} for the period of 01.01.2012 - 31.03.2017 (resulting in $t=1, ..., 1917$ time steps). The data is comprised of four German control areas in MW: 50 Hertz, Amprion, TenneT and TransnetBW. The summer time hour change was accounted for by creating hour 02 with interpolation between hours 01 and 03. For the clock change back the (later) 02 hour was deleted.
The German day-ahead total load data is calculated as an average of 4 x 15 min segments following the beginning of each hour. The "actual total load" for each control area (obtained at the end of the 15 min real bidding time in the balancing markets) is averaged and results across four regions are summed up to give the total actual load.
We use the steam coal ARA 1 month forward benchmark index for steam coal (one price per month) and the Germany Gaspool (GPL) natural gas day-ahead forward (one price per day) for gas. The weekly seasonality and holidays are included into a single dummy variable which takes on value 1 for Saturday/Sunday and the following German state holidays: New Year's Day, Good Friday, Easter Monday, Labour Day, Ascension Day, Whit Monday, German Unity day, Christmas Day, Boxing Day and New Years Eve (see Appendix (ref) for further details).
The intraday hourly electricity spot price data displays high volatility, fast changing dynamics and highly skewed distributions (see Appendix Figures (ref) - (ref)). The variance-covariance matrix of the hourly prices displays strong correlations, as demonstrated by highly-positive off-diagonal entries (see Appendix Figure (ref)). Thus, there is no independence between the hourly prices, and so the spread densities need to be modelled directly.
Because the hourly spread data possesses high degrees of skewness (see Appendix Figure (ref)) and kurtosis (see Appendix Figure (ref)), we model the full four parameter distribution of electricity price spreads at each time step (day), $t$. An intraday spread, $Y^{(s)}_{t}$, between two hours, denoted as spread number $s$, of the day-ahead electricity spot prices, is calculated by taking away the later hour price information from the earlier one, resulting in a positive spread if later hour is less expensive and negative otherwise. The full spread data set is comprised of a total possible $\frac{n(n+1)}{2}= \frac{24\times25}{2}=276$ non-duplicate entries, and is denoted by $\mathbf{Y} \in \mathbb{R}^{1917 \times 276}$. Each electricity spread time series data is tested for stationarity using the Augmented Dickey Fuller test (see Appendix Section (ref)) and is confirmed to be stationary at the 1% significance level. We plot example histograms for intraday spreads obtained between hours $00-08, 08-12, 12-16, 16-20$ (see Appendix Figures (ref) - (ref)) which depict high skewness and kurtosis of the spreads.
The spreads between exogenous variables are likewise formed by taking away the later hour values from the earlier ones in all cases except for the gas forward prices, coal forward prices and the dummy variable, for which only the daily, rather than hourly, values are available. Hence the following exogenous variables are considered for modelling $Y^{(s)}_{t}$: (1) spread of the lagged intraday electricity price, (2) gas Gaspool forward daily price, (3) coal ARA forward daily price, (4) spread of wind day-ahead forecast, (5) spread of solar day-ahead forecast, (6) dummy variable taking value of 1 for weekends/holidays, (7) spread of the day-ahead total load forecast, and (8) an interaction load variable, calculated as $Load_{spread} * \frac{1}{2}(Load_{earlierHr} + Load_{laterHr}) = \frac{1}{2}*(Load_{earlierHr}^2 - Load_{laterHr}^2)$ i.e. average of the load for the two hours from which the spread is calculated, weighted by the load spread obtained for those hours. This variable provides interaction of $\Delta Load * {avLoad}$ in order to account for the rate of change in load. The full exogenous variables spread data is stored in design matrix $\mathbf{X} \in \mathbb{R}^{1917 \times 9 \times 276}$, where the first column for each spread number $s=1,...,276$ is a vector of 1s needed for the calculation of an intercept.
The Generalised Additive Model for Location, Scale and Shape (GAMLSS) framework allows modelling each parameter (${\mu}, {\sigma}, {\nu}, {\tau}$) of the response variable's distribution as a function of explanatory variables. The range of possible distributions includes both exponential family (as per Generalised Linear Model (GLM)) and general family distributions, which allow for both discrete/continuous and highly skewed/kurtotic distributions (unlike GLM). The probability density function of response variable for $T$ observations is given by $f_{Y}(y^{(s)}_{t} | \boldsymbol{\theta}^{(s)}_{t}) \sim D(\boldsymbol{\theta}^{(s)}_{t})$, where $s=1,...,276$ is the spread number, $\boldsymbol{\theta}^{(s)}_{t} = [\theta^{(s)}_{t,1}, \theta^{(s)}_{t,2}, \theta^{(s)}_{t,3}, \theta^{(s)}_{t,4}]^T = [\mu^{(s)}_t, \sigma^{(s)}_t, \nu^{(s)}_t, \tau^{(s)}_t]^T$ is a vector of distribution parameters, and $D(\boldsymbol{\theta}^{(s)}_t)$ represents the distribution of spread $s$ on day $t$. We use parametric linear GAMLSS framework which relates distribution parameters to explanatory variables by
where $s = 1, ..., 276$ indicates the number given to each unique spread between two intra-day hours; $k=1,...,4$ specifies the distribution parameter corresponding to ${\mu}, {\sigma}, {\nu}, {\tau}$ respectively; $\boldsymbol{\theta}^{(s)}_k \in \mathbb{R}^T$ is a vector comprised of values for distribution parameter $k$ over $t=1,...,T$ time steps (i.e. $\boldsymbol{\theta}^{(s)}_1 = \boldsymbol{\mu}^{(s)}, \boldsymbol{\theta}^{(s)}_2 = \boldsymbol{\sigma}^{(s)}$, $\boldsymbol{\theta}^{(s)}_3 = \boldsymbol{\nu}^{(s)}, \boldsymbol{\theta}^{(s)}_4 = \boldsymbol{\tau}^{(s)}$); $T$ is the number of observations; $g_{k}(\cdot)$ is the monotonic link function of distribution parameter $k$; $\boldsymbol{\eta}^{(s)}_k \in \mathbb{R}^T$ is the linear predictor vector for distribution parameter $k$; $ \boldsymbol{\beta}^{(s)}_k =[\beta^{(s)}_{0,k},\beta^{(s)}_{1,k}, ..., \beta^{(s)}_{J_k,k}]^T \in \mathbb{R}^{J_k + 1}$ vector of coefficients learnt for parameter $k$; $J_k$ is the number of significant exogenous variables for parameter $k$ obtained at 5% significance level following the estimation and specification steps; and $\mathbf{X}^{(s)}_k \in \mathbb{R}^{T \times J_k+1}$ is the design matrix with each column containing spread data for significant independent variables.\\ Equation (ref) can be re-written for each distribution parameter
The choice of the link function influences how the linear predictor (i.e. systematic component $\mathbf{X}^{(s)}_k \boldsymbol{\beta}^{(s)}_k$) relates to each parameter. For example: a log link function for the standard deviation $g_{2}(\boldsymbol{\sigma}^{(s)}) = log(\boldsymbol{\sigma}^{(s)})$ results in the relationship $log(\boldsymbol{\sigma}^{(s)})= \boldsymbol{\eta}^{(s)}_2 = \mathbf{X}^{(s)}_2 \boldsymbol{\beta}^{(s)}_2$, hence the distribution parameter itself is obtained through a transformation $\boldsymbol{\sigma}^{(s)} = \exp(\mathbf{X}^{(s)}_2 \boldsymbol{\beta}^{(s)}_2)$.
Parameters $\boldsymbol{\theta}^{(s)}_k$ are calculated by maximizing the penalized likelihood using an algorithm which does not require calculation of the likelihood function's cross derivatives, but is a generalisation of the MADAM algorithm used for fitting mean and dispersion additive models, as in rigby1996mean.
The GAMLSS family provides a number of four parameter continuous distributions with pre-determined link functions. In order to model electricity spreads, the parameters of the distribution of choice must meet the following criteria: the standard deviation and kurtosis should take on positive values only, while the mean and skewness should be able to take on negative values since the spread data shows that both parameters can be positive or negative. Therefore distributions with the following link functions for each parameter are considered (see Table (ref)): mean: identity, standard deviation: log or logit, skewness: identity; kurtosis: log or logit.
There are seven continuous distributions within the GAMLSS framework which correspond to the necessary link function specifications: Johnson's SU ($\mu$ the mean), Johnson's original SU (JSU), Skew power exponential type 1 (SEP1), Skew power exponential type 2 (SEP2), Skew t type 1 (ST1), Skew t type 2 (ST2), and Skew t type 5 (ST5). The Box-Cox power exponential and Box-Cox t distributions are not suitable due to only being defined on the positive interval $Y_t \in [0,+\infty)$, while the electricity spread data is highly skewed and kurtotic (see Appendix Figures (ref) and (ref)). The Skew t type 3 and 4 are also unsuitable since they are only defined for positive skewness yet the spread data is often negatively skewed.
The expected value of a random variable $Y_t$ for each distribution in Table (ref) is given by $E(Y_t) = \mu_t + \sigma_t E(Z_t)$, where $z_t = \frac{y_t - \mu_t}{\sigma_t}$ is a normalised value of $y_t$, $Z_t$ is specified for each distribution separately, and $\mu_t,\sigma_t,\nu_t,\tau_t$ are the mean, standard deviation, skewness and kurtosis of the given distribution at time step $t$. See Appendices (ref) - (ref) for details on expected value calculations for each distribution.
In order to build accurate models of the spreads, we begin by analysing which of the seven possible distributions in Table (ref) fits each price spread data the best. We also analyse whether using a single distribution for modelling all of the spread data is a viable possibility and provide reasoning for when this may be beneficial. The analysis is divided into two main steps: (a) simple distribution fit, where each distribution of Table (ref) is fitted to the spreads in the training data set and the Akaiki Information Criterion (AIC) is used to assess the goodness-of-fit; and (b) factor-based distribution fit, where for each spread hour, the candidate distributions resulting from the analysis of step (a) are utilised to build a model using exogenous factors. The fit of these models is assessed using a validation data set and a number of goodness-of-fit measures.
The simple distribution fit is performed using GAMLSS function gamlss <- $\mathbf{y} \sim \mathbf{1}$, which results in a model under the specified distribution and time series $\mathbf{y}$. The analysis is performed using training data, comprised of the first 60% of the spread time series, $\mathbf{Y}_{train} \in R^{1150 \times 276}$. Therefore, for each spread number $s=1,...,276$ and each distribution $i=1,...,7$ of Table (ref), we build a model $\widehat{M}^{(s,i )} \leftarrow \mathbf{y}^{(s)}_{train} \sim \mathbf{1}$, resulting in $\widehat{ \mathbf{M} } \in \mathbb{R}^{276 \times 7}$ models. The AIC criteria of models $\widehat{ \mathbf{M} } ^{(s)} \in \mathbb{R}^{7}$, obtained for spread number $s$, are ranked in ascending order and the distribution corresponding to the model with the lowest AIC criterion is selected as the 'distribution of best fit'.
The results show that Skew t type 5 distribution is selected most often as the distribution of best fit based on this simple selection criterion (see Table (ref)). It is closely followed by the same family Skew t type 1 distribution. Overall, all of the distributions except JSUo are indicated as potentially of the best fit for some spread hours. A more detailed breakdown of the spreads for which each distribution was selected is displayed in Appendix Figure (ref). We use the result to further analyse the six possible distributions with a factor-based distribution fit method.
The results of simple distribution fit indicated that six out of seven continuous four parameter distributions could be used for modelling electricity spread data (see Table (ref)), which are $D^{(1)}$ - JSU, $D^{(2)}$ - SEP 1, $D^{(3)}$ - SEP 2, $D^{(4)}$ - ST1, $D^{(5)}$ - ST2, and $D^{(6)}$ - ST5. Hence we proceed with a factor-based analysis by fitting a regression to each parameter of candidate distributions using exogenous variables within the GAMLSS framework. The training data set is comprised of the dependent $\mathbf{Y}_{train} \in R^{1150 \times 276}$ and independent $\mathbf{X}_{train} \in R^{1150 \times 9 \times 276}$ variables, where for each spread number $s$, under distribution $D^{(i)}, i=1,...,6$, initial equations for the central moments are
where $\bm{\widehat{\beta}}_k =[\widehat{\beta}_{k,0},\widehat{\beta}_{k,1}, ..., \widehat{\beta}_{k,8}]^T \in \mathbb{R}^{9}$ is initial vector of coefficients for distribution parameter $k$, $\mathbf{x}_t =[1,x_{1,t},...,x_{8,t}]^T \in \mathbb{R}^{9}$ is the initial vector of independent variables where $x_1$ is the spread of lagged day-ahead electricity price, $x_2$ is the gas Gaspool forward daily price, $x_3$ is the coal ARA forward daily price, $x_4$ is the spread of wind day-ahead forecast, $x_5$ is the spread of solar day-ahead forecast, $x_6$ is the dummy variable taking value 1 for weekends/holidays, $x_7$ is the spread of the day-ahead total load forecast, $x_8$ is the interaction load variable.
The models are specified using iterative updating of the equations for each moment, where the refinement is performed by deleting the most insignificant variable one-by-one and re-estimating the model, until all variables are significant at 5% (see Algorithm (ref)). This results in $\mathbf{\widehat{M}} \in \mathbb{R}^{276 \times 6}$ models containing the estimated coefficients. Note: the process revealed that some distributions were not suitable for modelling certain spreads within the GAMLSS framework. Convergence was not achieved for some distribution parameters, typically $\tau$ and occasionally $\mu$ (see Appendix Section (ref)) and when this happened the distribution was omitted from candidacy for that spread.
The models estimated using the factor-based distribution fit are analysed in four ways: (1) producing the expected value fit over the training data, (2) producing the expected value fit over validation data (comprised of the next 20% of unseen time series, at data points $t=1151,...,1534$), (3) analysing the goodness-of-fit over the validation data using Root Mean Squared Error, and (4) analysing the goodness-of-fit over validation data using Pinball Loss function measure.
The fitted distribution parameters $\widehat{\boldsymbol{\theta}}^{(i)}_{train} \in \mathbb{R}^{1150 \times 4 \times 276}$ for each distribution $i=1,...,6$ are used to find the fitted expected value, $E(\widehat{\mathbf{Y}}_{train})$, of the training price spread data for $s=1,...,276$ over training data points $t=1,...,1150$ using Equations (ref) - (ref). Illustrative examples of the fitted expected values for spreads between hours 00-08 and 08-12 across the six possible distributions are depicted in Appendix Figures (ref) and (ref) respectively. The true price spread values, $E(\mathbf{Y}_{train})$, are given by blue lines and the fitted, $E(\widehat{\mathbf{Y}}_{train})$, by red lines. The plots show a good fit to the true data across all distributions with slight underestimation in the spread price. In the case of the two example spreads, the spikes in the data seem to be fitted best by different distributions. The spread hour 00-08 spike at time step $t=360$ is only fitted using the SEP2 distribution, while spike of spread 08-12 at time step $t=409$ is fitted best with the ST2 distribution. This supports the analysis of the six possible distributions for candidacy of best fit to individual spread data.
The fitted models $\widehat{\mathbf{M}} \in \mathbb{R}^{276 \times 6}$, containing the estimated coefficients for each spread number $s=1,...,276$ under each distribution $D^{(i)}, i=1,...,6$, are used to forecast the distribution parameters, $\widehat{\boldsymbol{\theta}}^{(i)}_{validate} \in \mathbb{R}^{ 383 \times 4 \times 276}$ over validation data. We note that the same estimated model $\widehat{M}^{(s,i)}$ is used to build the forecasts over validation time series, i.e. each vector of estimated coefficients $\boldsymbol{\widehat{\beta}}^{(s,i)}_k$ for distribution parameters $k=1,...,4$ is re-used at each time step $t$ to make the predictions. Once the distribution parameters are forecasted the expected value of the spread is calculated using Equations (ref) - (ref). Illustrative examples of the forecasted expected values for spread hours 00-08 and 08-12 across the six possible distributions are depicted in Appendix Figures (ref) and (ref) respectively. The true spread prices, $E(\mathbf{Y}_{validate})$, are given by blue lines and the forecasted, $E(\widehat{\mathbf{Y}}_{validate})$, by red lines. The results for spread 00-08 show that all distributions are able to fit the validation data well, however the downward spikes seem to be underestimated across all distributions. The expected value forecast for the spread hour 08-12 shows that some distributions fit the data better than others, for example ST1 does not fit the data as well as ST5.
A common performance measure, Root Mean Squared Error (RMSE), is used to assess the goodness of fit for the forecasted expected value of the spreads over the validation data set. The RMSE is calculated for the forecasted expected values of each spread, $s=1,...,276$, obtained under each distribution $D^{(i)},i=1,...,6$, using
Since the expected value at each time step $t$ is calculated from all of the four forecasted parameters $\big( \widehat{\mu}^{(s,i)}_t, \widehat{\sigma}^{(s,i)}_t, \widehat{\nu}^{(s,i)}_t, \widehat{\tau}^{(s,i)}_t \big)$, the RMSE measure provides a goodness of fit based on the forecast for the entire distribution specification. The results are summarised in Table (ref), and show that ST5 distribution was used to form models corresponding to the lowest error across $\frac{110}{276}*100=40\%$ of the spreads. This is in line with the finding of the simple distribution fit over the spread training data, which indicated that ST5 was chosen as the best distribution most often among suitable four parameter distributions (see Section (ref)). A detailed breakdown of best distribution assignment for each spread number and the corresponding RMSE values are shown in Appendix Figures (ref) and (ref) respectively. The results show that spreads with hours 13.00, 14.00 are the hardest to forecast since they have the highest RMSE error (dark red).
Our model fits the entire four parameter distribution to the spread data at each time step. Therefore we consider it more appropriate to evaluate the goodness-of-fit using an entire fitted distribution as represented by quantiles. The Pinball Loss (PL) function is typically used as the objective function in quantile regression and can also be interpreted as the accuracy of a quantile forecasting model. We adopt this measure as our main performance metric to select the best distribution to each spread. We follow Algorithm (ref) when making an assessment of each predictive density power and outline the steps involved in calculating the performance measure below.
The Pinball Loss Performance Measures are analysed for each spread $s$ and the results are presented in a $24 \times 24$ upper diagonal matrix, containing the best distribution number selected for each intraday spread between the hours indicated by row and column labels (see Figure (ref)). The figure shows that majority of the spreads are fitted best with ST5 distribution (number 6 - purple colour). This is in line with the results obtained from the factor-based distribution fit using the RMSE measure (see Appendix Figure (ref)), and with the simple best distribution fit given by $\mathbf{y} \sim 1$ function (see Appendix Figure (ref)), which both favoured the ST5 distribution.
Table (ref) specifies the number of times each distribution was selected as best from the 276 possible spreads. While ST5 is selected for almost 50% of the spreads, the other 5 distributions were selected nearly with the same proportion for the remaining half of the spread data.
Further to this, we calculate the % difference of the best distribution PL performance measure value vs that of the next best distribution, $\frac{|{\cal{L}}^{(s)}_{1} - {\cal{L}}^{(s)}_{2}|}{{\cal{L}}^{(s)}_{2}} * 100\%$, for each spread number $s$. When ST5 was selected as the best distribution, the average difference of the performance measure was 4.84% compared to, when other distributions were selected as best, the average difference was 1.44% (with approx 1/3 of these cases containing ST5 is the second best). This makes the case for ST5 to potentially be used as a general best fit distribution across all spreads.
Next, we re-estimate models based on the best distributions established for each spread number, analyse performance using Rolling Window forecast technique and evaluate results using Pinball Loss function utilised within a formal significance test framework, based on Diebold-Mariano test.
We perform a careful forecast analysis based on the best distribution selected for each spread number (see Figure (ref)) using rolling window technique, and compare the results to the performance of models obtained with a Normal two-parameter benchmark using Diebold-Mariano test. The estimation horizon (rolling-window size, $T_0$) comprises of 80% of the data ($1534$ observations), which is moved up 1 time step at a time. For each spread number $s=1,...,276$ we re-specify and re-estimate the model over rolling window time frame $T_0$ using the chosen distribution $D^{(s)}$ and create day-ahead forecasts.
The rolling window forecast analysis is performed over the unseen data points $t=1535,...,1917$ and depicted in Figure (ref), with steps outlined in Algorithm (ref). For each spread $s=1,...,276$ a model is specified and estimated over a fixed length horizon of $T_0=1534$ observations starting at time step $t_1$ (green coloured bracket). The equation specification for each distribution parameter is obtained using iterative improvements, where upon each iteration the most insignificant variables are deleted one-by-one (except intercept) until all variables in each equation are significant at 5% level. Following this a 1-step ahead forecast, f$_1$ (green colour), is made using the estimated model and a full predictive density is obtained. The window is moved by 1 time step to $t_2$ (blue colour) and the procedure is repeated. The rolling window analysis is continued until all forecasts for the current spread are created, comprising 383 data points. The procedure was parallelised across 8 cores on an Intel Core i7 2.9 GHz processor and took 216 hours to complete.
The rolling window analysis of Algorithm (ref) was repeated using Normal distribution for all $D^{(s)}$ resulting in benchmark forecasts. We note that the estimation of models was not always convergent within the GAMLSS framework and certain spreads failed to be fitted with our models, with Normal distribution based models or both, which was mostly due to the dataset. For example, the spread obtained between hours 02-03 is dominated by noise, with $\frac{1}{3}$ of the data set comprising of $0$'s, which creates a problem during weight update calculations.
The results are analysed using two stages: (1) obtaining the difference between Pinball Loss performance measures of the forecasts made with our models and those made by the benchmark; and (2) performing an official statistical significance (Diebold-Mariano) test on whether the difference between two Pinball Loss performance measures is statistically significant.
First we obtain the PL performance measures of our models and benchmark using Algorithm (ref), where for the benchmark $D^{(s)}$ we use the Normal distribution for all spreads.
If for a given spread number $s$ convergence issues were experienced during analysis of our models, the forecast for that time step was omitted (and accounted for in PL average calculations). Number of times at least 1 forecast was missing for skew type distribution was 46 times (at most 70 time steps out of 383). If more than 200 out of 383 forecasted time steps were missing, the result was judged as unreliable and forecasts for that spread were treated as unavailable; this happened 4 out of 276 times for spreads between hours: 02-03/12, 03-05, 13-14. The convergence issue was also experienced for the Normal benchmark, and the number of times at least 1 forecast was missing was 5 times (at most 42 time steps out of 383). The number of times forecasts were missing for more than 200 time steps were observed more frequently, which happened 13 out of 276 times for spreads between hours: 01-12/13/14/23; 02-03/23; 05-13; 06-12/13/14/15/16; 07-08. When this happened we judged our model to have a superior result for that spread number $s$.
We report the difference in PL performance measures between those obtained with our models and the benchmark for each spread number $s$ (see Figure (ref), where negative values (blue, green) indicate that our model outperformed the benchmark obtained with the Normal distribution for that spread number $s$). Initially results are obtained to 5 d.p. and these show that for 276 spreads our models forecasted the full density more accurately: 258 times vs 18 times for the Normal distribution. Reducing decimal places to 3, results in the same number of Normal distributions having smaller error values (18 times), however now results show that for 5 spreads the Normal distribution produces as good a fit as the distributions used in our models, while for 253 cases our models produce more accurate forecasts. The overwhelming majority preference for models obtained with skew type distributions is evident through negative values, as expected due to highly skewed and kurtotic nature of spread price data. Next, we seek to establish a formal significance test on the difference between the PL performance measures between our models and benchmark.
In order to draw statistically significant conclusions over the outperformance of the best selected distributions when forecasting unseen data points over the accuracy of forecasts made with the Normal distribution, we use the Diebold-Mariano (DB) test diebold2002comparing. The test is applicable to forecast errors that do not have zero mean, that are not Gaussian, and that may be serially / contemporaneously correlated.
We use a variation of the standard DB test with implementation proposed by harvey1997testing. For each spread number $s$, there are $t=1,...,383$ forecasts produced by two models ($\widehat{M}^{(s,t)}_1$ - estimated with the best chosen distribution and $\widehat{M}^{(s,t)}_2$ - estimated with Normal distribution), which are tested against each other using a one-sided test at 5% significance level. The null hypothesis is the that the two models have the same forecast accuracy, with a one-sided alternative hypothesis that forecasting power of the best distribution outperforms that of the Normal, $H_0 : E(\Delta_{M_1,M_2,t,s}) \leq 0$, where the Loss Differential Series $\Delta_{M_1,M_2,t,s}$ is
where $s$ is the spread number, $t$ is the forecast time step, $\bar{L}_t^{(s,1)}$ average quantile score of model $\widehat{M}^{(s,t)}_1$ obtained with the best chosen distribution at forecast step $t$, $\bar{L}_t^{(s,2)}$ average quantile score of model $\widehat{M}^{(s,t)}_2$ obtained with Normal distribution at forecast step $t$.
The p-values are displayed in Figure (ref) and the results show that out of 276 spreads: (a) the best chosen distributions are significantly better at forecasting the spreads: 161 times at 5% (bright green) and 15 times at 10% (olive green). Note: the distributions which were used to learn these 176 models are: JSU - 9 times, SEP1 - 14 times, SEP2 - 16 times, ST1 - 13 times, ST2 - 18 times and ST5 - 106 times ($\frac{106}{164}*100 = 64.5\%$). (b) the models have the same forecasting power for 48 spreads (i.e. the null hypothesis could not be rejected at 10%). Note: the distributions which were used to learn the corresponding 48 models are: JSU - 9 times, SEP1 - 4 times, SEP2 - 8 times, ST1 - 8 times, ST2 - 4 times and ST5 - 15 times ($\frac{15}{48}*100 = 31.25\%$). (c) the results could not be obtained for the models obtained with best chosen distributions 49 times (white spaces of upper triangular) due to at least 1 forecast missing due to Quantile estimate convergence issues. Note: the distributions which were used to learn these 49 models are: JSU - 4 times, SEP1 - 8 times, SEP2 - 21 times, ST1 - 9 times, ST2 - 5 times and ST5 - 2 times ($\frac{2}{49}*100 = 4\%$) (i.e. this typically happened for distrubitons other than ST5).
We focused on selecting suitable four parameter distributions which fitted best the individual spread data and performed detailed analysis of the forecasting power of such distributions. As the result of the above analysis it was found that: (a) the ST5 distribution was selected most frequently as the best distribution by simple distribution fit based on $\mathbf{y} \sim \mathbf{1}$ and factor-based distribution fit based on both RMSE and PL functions; (b) the ST5 distribution was the most reliable for convergence of quantile estimates using GAMLSS function qFUN (e.g. qST5), especially for the extreme quantiles of $q_1, q_2, q_3, q_{97}, q_{98}, q_{99}$ where other distributions such as SEP1, SEP2 often failed. Quantile estimates are important for our research because they are utilised in statistical testing for comparing performance of two models and in the Value-At-Risk calculations used in trading strategy optimisation. For example, the PL calculations revealed that for 49 spreads at least one forecast step (out of 383 steps) failed to extract 95 quantiles from the estimated model. Upon further examination it was revealed that out of the 49 occurrences only 4% had ST5 as the underlying distribution for which the model was estimated. This supports our claim that ST5 is a reliable distribution for the quantile estimates; (c) the best distribution for each spread was chosen based on the PL performance measure calculated over the validation data. Each spread had 6 possible distributions from which the best distribution was chosen and the one with the lowest score was taken as the best distribution for that spread. Further analysis reveals that on 1/3 of occasions when other distributions than ST5 were selected as 'best', the ST5 was the second best distribution, which on average was only worse by 1.44%. However for the spreads where the ST5 was the best distribution, the PL performance measure was on average better by 4.84%, which points to the possibility that the best distribution did not have a significantly difference performance when compared to ST5. Therefore we conclude that if one wishes to use a single distribution across all spreads, the ST5 distribution forms a robust choice. We continue our analysis using a more detailed approach, where individual spreads have established distributions of best fit assigned to them.
In order to validate the need for dynamic modelling, we display the evolution of fitted moments for four example spreads obtained for the first rolling window (i.e. first 1534 data points used for specification and estimation phase), with their associated distributions: 00-08 (ST1); 08-12 (ST1); 12-16 (ST1); 16-20 (ST5). \\ Evolution of latent 4 central moments throughout four years i.e. examining how selected 4 spreads behave throughout each year. We selected four years: 2012, 2013, 2014, 2015 in order to depict the evolution and the changing dynamics of the moments with time. The evolution of each distribution parameter is plotted on separate graphs (see Figure (ref) for the evolution of $\hat{\boldsymbol{\mu}}$, Figure (ref) for evolution of $\hat{\boldsymbol{\sigma}}$, Figure (ref) for the evolution of $\hat{\boldsymbol{\nu}}$, Figure (ref) for evolution of $\hat{\boldsymbol{\tau}}$). The results show that the mean is fitted in line with what would be expected for the true spread price, where the spreads between 08-12 hours tend to be positive (i.e. later hour is cheaper), while the 16-20 spreads are negative, i.e. later hour is more expensive. The standard deviation is highest for the spreads between night-time (less busy) and early morning / afternoon (green and blue lines). While the skewness tends to be positive for the 08-12 hours, and negative for the 00-08 hours.
Next, we plot the true $E(\mathbf{Y})$, vs fitted $E(\widehat{\mathbf{Y}})$, expected values of spreads over the four years, selecting different spread hours to the ones used above, in orer to show a variety of underlying distributions used: 00-09 (SEP2), 08-11 (SEP1), 11-19 (ST2), 16-22 (ST5). Figure (ref) shows the true evolution of the spreads, compared to the expected values produced by the estimated models (see Figure (ref)), note slightly smaller scale. The fitted expected values follow the true pattern throughout each year, for example the 16-22 hour spread tends to have more negative values in the summer time (i.e. electricity at earlier hour is less expensive) and positive values in the winter time (i.e. electricity at earlier hour is more expensive). It can be seen that the fitted spread values are slightly under-fitted as indicated by the difference in plot scale.
Evolution of the 4 latent central moments throughout a day demonstrates how distribution parameters change throughout a day and throughout different times of the year. We plot the four parameters for 276 spreads on 01 Jan 2015, 01 Mar 2015, 01 June 2015 and 01 September 2015 i.e. one plot per day of the season. The spreads were plotted sequentially starting with 23 spreads for midnight hour with all other hours of the day, $00-01, 00-02,...,00-23$, followed by 22 spreads of hour 01 with all other hours, continuing on until the last spread between hours 22 and 23. A total of 276 spreads are plotted for each of the four selected days of year 2015. The results show periodic spikes due to points where the spread moves between hour 00 with all other, hour 01 with all other etc. We show the changing dynamics of the distribution parameters throughout seasons of the year (see Figure (ref)).
Plots detailing the same intra-day information but only for spreads of 4 chosen example hours of the day with all other hours (00-; 08-; 12-; 16-) are given below for clarity. Effectively, for example, the Figure showing the mean for different days of the year 2015 (see Figure (ref)) has sub-plots (a-d) which demonstrate zoomed in sections of Figure (ref) depicted by the grey dotted box. Note: the sub-plots $x$ values range reduces when going from (a) to (d) since midnight hour 00 spreads with all hours of the day, but 16 hour only spreads with 7 other hours.
Variations in the size and sign of explanatory variable coefficients for different spreads of the same day illustrate the varying impact of drivers (see Tables (ref) and (ref) which show the estimated coefficients for selected 4 spreads: 00-08; 08-12; 12-16; 16-20). The coefficients were extracted from the model estimated for the first rolling-window frame i.e. $t=1,..., 1534$. The displayed values for the coefficients are all significant at 5% thus forming the equation for that moment. Missing values indicate that the independent variable was not significant and thus was omitted from the equation specification for that moment (table column).
The coefficients correspond to their associated independent variables: $\beta_0$ coeff is the intercept; $\beta_1$ coeff for spread of the lagged day-ahead electricity price; $\beta_2$ coeff for gas Gaspool forward daily price; $\beta_3$ coeff for coal ARA forward daily price; $\beta_4$ coeff for spread of wind day-ahead forecast; $\beta_5$ coeff for spread of solar day-ahead forecast; $\beta_6$ coeff for dummy variable taking value of 1 for weekends/holidays; $\beta_7$ coeff for spread of the day-ahead total load forecast; $\beta_8$ coeff for an interaction load variable.
Overall the signs and significances of the coefficients are intuitive. In particular, wind and solar production spreads have negative effects on the mean and skewness of spreads for the morning and afternoon spread pairs shown. Recall that the spread is defined as the former minus the later hours and so a higher wind and solar production spread will generally reduce the the average spreads and also the skewness. This is consistent with the effects of wind and solar production on price levels reported in gianfreda2017stochastic and elsewhere.
The average Root Mean Squared Errors for the expected values of spreads forecasted with models using the best and Normal distributions over the rolling-window forecasting horizon are given in Table (ref). Each score is calculated by averaging RMSE values across all forecasted spreads (note: if at least 1 time step forecast was missing for a given spread, the forecasted time series was omitted from RMSE calculation). A detailed breakdown of the RMSE values for the best and Normal distribution models is given in Appendix Figures (ref) and (ref) respectively. The RMSE values are in line with, but slightly higher overall, than the ones found for the validation data set (see Figure (ref)). The RMSE based forecasting power evaluation concludes that the best distribution and Normal distribution models produce compatible results, without significant differences.
The battery operation trading schedule is optimised by maximising the trading profit using forecasted spread densities for the day-ahead hourly electricity prices, while meeting a risk criterion of making a profit per MWh traded of more than $c$ per daily cycle with 95% confidence. We consider $c$ to be the round trip transaction cost for a storage facility in charging and discharging. This covers the technical efficiency loss of the battery between the charging and discharging of energy, as well as the transmission, distribution, trading, balancing, levies and other use of system costs for a facility seeking to operate in the wholesale market. Estimates of these costs vary widely in practice and so for comparison we use a sensitivity analysis approach with $c = 5, 10, 15$ Euro/MWh. We impose an assumption of a finite horizon daily operation in the day-ahead forecasting analysis, which constrains the opening and closing battery charge levels to be equal each day. We investigate the optimal opening/closing level of charge level by performing calculations for a number of initial battery levels, $b = 0, 0.1 ,0.2,...,0.9$, and selecting the one with the highest profit. Therefore, the schedule is optimised with respect to initial battery levels, $b$, of the total capacity of a nominal 1 MWh battery, which could trade in the wholesale market. The battery is assumed to be fully (dis)chargeable within 1 hour of trade execution.
The realised Profit and Loss (P&L) for a day's trade, $PNL_t$, executed using spread hour $s$, is calculated as per Equation (ref) and is based on realised spread value and total round-trip cost.
We backtest the trading schedule over approximately 1 year of data (383 days, time steps $t=1535,...,1917$) and report for each starting battery level: P&L over the backtest period, $PNL$ (Eq. (ref)), average P&L over the backtest period, $\overline{PNL}$, (Eq. (ref)), standard error of the P&L average, $s^{\overline{PNL}}$ (Eq. (ref)), number of trades which resulted in a loss after all costs are taken into account, $n_l$, (Eq. (ref)), total monetary value resulting from loss days, $l$, (Eq. (ref)), and average loss, $\overline{l}$, (Eq. (ref)).
Figure (ref) (a) depicts a hypothetical trade where the battery is charged to 50% level at the start of the day. The example trade consists of a charge at 00 midnight from 50% to full capacity of 100% (green line - 1 hour charge), followed by a holding period (grey) until discharge at hour 13 (red line) down to 50% charge level, thus completing a single trade for the day (consisting of two legs of the spread), returning the battery to the initial 50% charge level. An alternative trade is displayed in sub-plot (b) where the initial battery level is 90%. This charge level is held until discharge at hour 13 down to 0% charge level, followed by a re-charge at hour 14 back to 90% charge, thus completing a single trade for the day.
Initially, the optimal trade to be performed on each day is established using an exhaustive state space search over all available trade executions, performed using Algorithm (ref) (note, the same procedure is repeated for Normal two parameter $(\mu_t,\sigma_t)$ distribution estimated models to benchmark the results). At each test point time step, $t=1,...,383$, i.e. 1 day of trading, we consider $s=1,...,276$ potential trades based on forecasted spreads, where for each trade we examine with 95% confidence the possibility of making a profit of at least $c$ Euro/MWh.
First, we calculate the forecasted expected value of the spreads at each time step $t$ using Equations (ref) - (ref) (note: $E(\widehat{Y}^{(s)}_t) = \hat{\mu}^{(s)}_t$ for distributions JSU and ST1, see Appendix (ref) for calculating expected values of the other distributions used). Occasionally the forecasted kurtosis was very large and when this happened it was capped at the value of 100 in order to successfully use the equations for expected value calculation. Next we establish whether the expected value $E(\widehat{Y}^{(s)}_t)$ of a spread number $s$, is positive (later price is lower, therefore discharge then charge) or negative (later price is higher, therefore a profitable trade would require charging first, followed by discharging). If the expected value of the spread was missing, the calculations for that spread were skipped (note this happened a total of 1795 for models based on the best chosen distributions and 4159 for models based on Normal benchmark distribution, out of $276\times383=105,708$ possible times). For each case, we calculate the critical value corresponding to 95% confidence interval. If forecasted expected value of the spread $s$ on day $t$ is positive, $E(\widehat{Y}^{(s)}_t) > 0$, we access the 5$^{th}$ quantile, $q_5$, since the body of the distribution is to the right of this critical value (later price is lower, hence discharge then charge). If $E(\widehat{Y}^{(s)}_t) < 0$ we access the 95$^{th}$ quantile so that body of the distribution is to the left of the critical value (i.e. later price is higher, hence charge then discharge).
Next, we calculate the forecasted profit for each spread at time step $t$, only if the critical value obtained from the quantile estimation exceeds the cost of $c$ Euro/MWh. The forecasted profit is found as the weighted difference of the forecasted expected value and the total roundtrip cost. The trade corresponding to the maximum P&L value at time $t$ is selected as the optimal trade. The realised profit is calculated according to the realised outcome for that spread (see Equation (ref)) and the profit & loss analysis is performed as per Equations (ref) - (ref) for each trade and with respect to the benchmark Normal type distribution. Note 1: we do not trade on any day, for which the forecasted profit is below the roundtrip cost $c$ Euro/MWh. Note 2: if a quantile fails to be estimated for a given model, the critical value $\hat{q}_{xx}$ gets set to 0 (for the best chosen distributions models, this happened at forecast time steps, $t = 251, 264, 304, 349$).
The roundtrip cost of 5 Euro/MWh is considered first, the results of which are reported in Table (ref) (note: models based on the best chosen distribution did not get used for trading on 2 days, while those using the Normal distribution did not trade on 0 days).
The results of 5 Euro/MWh roundtrip cost trading indicate that our method and the benchmark approach both have the same predictive power, with average P&L results for each initial battery level not having significantly different results at 95% level. The best initial battery level transpires to be 0% charge, which is explained by the fact that typically electricity would be cheaper during the night due to lack of demand (which makes it cheaper to charge the battery), and thus making the most profit by discharging at some point during peak demand in the day. The total P&L value across all battery levels of Normal type models is better by approximately $\frac{44384-43856.5}{43856.5}*100= 1.2$% over the backtest period, however this does not reflect the number of encountered loss days and their total monetary value. The ratio of total number of loss days is $\frac{71}{15}=4.7$ times more for the Normal type, with the monetary value of these losses $\frac{143.8}{19.9}=7.2$ times higher over models based on the best chosen distributions.
Overall, the 5 Euro/MWh cost basis is not large enough to reveal the advantage of using our method when making trading decisions, however this becomes obvious when the cost basis is increased to more realistic scenarios of 10 and 15 Euro/MWh (see below).
The P&L results for 10 Euro/MWh roundtrip cost trading scenario are reported in Table (ref) (note: models based on the best chosen distribution did not get used for trading on 57 days for b=0 and 49 days for all other starting battery levels, while Normal did not trade on 62 days).
The higher roundtrip cost of 10 Euro/MWh, reveals the advantage of using Skew type and similar distribution based models: the average P&L is significantly better than Normal model's at 95% confidence, and the best initial battery charge level is again 0% charge level (i.e. fully discharged). The P&L over the tested period across all battery levels is $\frac{29469-25099.2}{25099.2}*100 = 17.4$% higher using our approach over the Normal. The ratio of total number of loss days for Normal to Skew models is $\frac{262}{80}=3.3$, with the monetary value of losses being $\frac{482.1}{108.9} = 4.4$ times higher for the Normal types.
The P&L results for 15 Euro/MWh roundtrip cost trading scenario are reported in Table (ref) (note: models based on the best chosen distribution did not get used for trading on 223 for b=0 and 210 days for all other starting battery levels, while Normal did not trade on 255 days (out of the total tested period of 383 days).
The average P&L is significantly better than Normal model's at 95% confidence, and the best initial battery charge level is again 0% charged. The P&L over the tested period across all battery levels is $\frac{15075.4 - 9062.8}{9062.8}*100 = 66.3$% higher for the Skew and similar distribution types selected as best, over the Normal. The ratio of total number of loss days for Normal to Skew models is $\frac{229}{123}=1.9$, with the monetary value of losses being $\frac{411.5}{140.8} = 2.9$ times higher for the Normal types.
We have demonstrated the value of detailed, computationally-intensive modelling of intra-day power price spread densities using a flexible four parameter distributional form, generally the skew-t. This allows the dynamic conditional parameter estimates to follow stochastic evolutions driven by exogenous factors, most importantly the day ahead demand, wind and solar forecasts. These forecasts fit well in backtesting and support the optimal daily scheduling of a storage facility, operating on a single cycle. The model outperforms baseline comparisons to a normal density model. The model specification and validation process is computationally intensive, and whilst modelling simplifications could be introduced, accuracy is important. Overall, this formulation and application shows the merits of a computationally intensive approach to accurate specifications and these are likely to be more attractive in practice than methods based upon analytical simplifications of the stochastic price processes. The optimal choice of spreads to trade do vary daily and the need to utilise forecasts in well-specified models is evident, as is the delicate balance between expected profits and risk. Furthermore, the algorithmic nature of the modelling presented would lead naturally to the potential for algorithmic trading by battery asset owners, which may be more economical to small enterprises than outsourcing their trading to larger service providers. In terms of optimisation, we have indicated the value of the methodology in supporting optimal storage on a one cycle per day basis. Further extensions to two or more cycles per day is clearly possible and would most likely further endorse the value of accurate spread density specifications.