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.
75,342 characters · 27 sections · 93 citation commands
Deep Distributional Time Series Models and the Probabilistic Forecasting of Intraday Electricity Prices
{0.15cm} {0.15cm} \pagestyle{empty}
\pagestyle{plain} \setcounter{equation}{0}
Deep models with rich feature vectors are often very effective in problems which require accurate forecasts goodfellow+bc16. These include financial applications, such as predicting equity risk premiums and returns feng2018,gu2020,gu2020autoencoder and bond returns bianchi2020. Another financial application where deep models have high potential is the forecasting of intraday electricity prices. Electricity prices exhibit a strong and complex nonlinear serial dependence, quite unlike security prices, and accounting for this is key to obtaining accurate forecasts nowotarski2018,manner2019. Shallow neural networks (NNs) amjady2006,mandal2007, and more recently deep neural networks (DNNs) lago2018,ugurlu2018, have been shown to capture these nonlinearities well and produce accurate point forecasts. However, it is the accurate forecasting of the entire distribution of prices---variously called probabilistic, density or distributional forecasting---that is important for both market operators and participants. Yet, to date, probabilistic forecasts of electricity prices using DNNs are rare. In this paper, we propose a number of time series probabilistic forecasting models that exploit and extend state-of-the-art deep models, and apply them to data from the Australian market.
Day-ahead wholesale electricity markets operate throughout the world, including in the U.S. and Europe. In these markets, generators and distributors bid for sale and purchase of electricity at an intraday resolution in an auction one day prior to transmission. The auction clearing price is widely called the electricity spot price; see kirschen2018 for an introduction to such markets. Accurate price forecasts at an intraday resolution, one or more days ahead, are central to both the efficient operation of the market and profitability of participants. Particularly important are probabilistic forecasts of price, not just the mean, variance or other moments. This is because overall profitability of market participants is strongly affected by prices in the tails, which are very heavy in most wholesale markets.
Recurrent neural networks (RNNs) are DNNs tailored to capture temporal behavior, and are suitable for forecasting nonlinear time series goodfellow+bc16. However, RNNs typically have a very large number of hidden weights and are difficult to train and tune. Therefore, we use a variant of RNNs called echo state networks (ESNs) Jae2007,LukJae2009, that are flexible and employ a form of regularization that makes them scalable to long series and computationally stable. We build statistical time series models based on ESNs using two approaches. The first extends that of chatzis2011echo and mcdermott+w17,mcdermott2019deep, who use ESNs within statistical models, with the output layer coefficients of the hidden state vector estimated using Bayesian methods. In our work we include a shrinkage prior for these coefficients to provide additional regularization, along with three different additive error distributions for the output layer: Gaussian, skew-normal and skew-$t$.
Our second approach is the main methodological contribution of the paper. It uses the implicit copula of the time series vector from a Gaussian probabilistic ESN of the type described above. By an “implicit copula” we mean the copula that is implicit in a multivariate distribution and that is obtained by inverting the usual expression of Sklar's theorem as in nelsen06. This implicit copula is both a deep function of the feature vector, and also a “copula process” with the same dimension as the time series vector. We combine our proposed copula with a non-parametrically estimated marginal distribution for electricity prices, producing a time series model that captures the complex serial dependence in the series. An accurate estimate of the marginal distribution ensures “marginal calibration”, which is where the long run average of the predictive distributions of the time series variable matches its observed margin gneiting+br07,gneiting2014. Importantly, the entire predictive distribution from the copula model is a deep function of the feature vector. We note that our copula model extends the deep distributional regression methodology of KleNotSmi2019 to deep distributional time series and ESNs.
In both our deep time series models the feature vector includes a rich array of past series values and possibly other variables. To regularize these, ESNs use sparse and randomly assigned fixed weights for the hidden layers of the DNN. For each of $K$ random configurations of weights, Markov chain Monte Carlo (MCMC) is used to estimate the statistical model and compute Bayesian predictive distributions. The probabilistic forecasts are then ensembles of these predictive distributions over the $K$ configurations of weights. We show in our empirical work that this ensemble provides for accurate uncertainty quantification.
We use our deep time series models to forecast intraday electricity prices in the Australian National Electricity Market (NEM). The NEM is one of the earliest established wholesale markets (in 1998), comprises approximately 1% of all Australian economic output, provides publicly available data, and has a design that is typical of other day-ahead markets. It has five regional price series for which forecasting has been much studied; see ignatieva2016, smith2018, manner2019 and han2020 for overviews. Using a feature vector with lagged prices from all five regions, we compute predictions at the half-hourly resolution for a 24 hour horizon over an eight month validation period during 2019. serinaldi2011,gianfreda2018 and narajewski2020 show the `generalized additive models for location, scale and shape' (GAMLSS) methodology of rigby+s05 applied to time series allows for the accurate modeling and forecasting of electricity prices, and we employ this as a benchmark. Using contemporary metrics, the probabilistic ESN with skew-$t$ errors produces more accurate point and probabilistic forecasts than ESNs with either Gaussian or skew-normal errors, and also the GAMLSS benchmark. However, the deep copula model produces substantially more accurate probabilistic forecasts, including in both tails, and its forecasts have superior coverage. Thus, marginal calibration also improves calibration of the (conditional) predictive distributions.
Participants in the NEM are provided with probabilistic forecasts of electricity demand by the system operator. Recent studies ziel2016,shah2020 suggest that using accurate demand forecasts may further improve time series forecasts of the distribution of price. An advantage of deep models is that additional predictors are easily included in a flexible fashion as extra elements in the feature vector. We do so here using three quantiles of the 24 hour ahead demand forecasts, and find this increases forecast accuracy of the upper tail of the price distribution, but only when using the deep copula model. Accurate forecasting of the upper tail of electricity prices impacts the profitability of participants substantially christensen2012.
The paper is organized as follows. Sec. (ref) provides an overview of electricity markets and price forecasting, with a focus on the NEM. Sec. (ref) outlines ESNs, their probabilistic extension using additive disturbances, and Bayesian methods for their estimation and prediction. Sec. (ref) outlines the implicit copula process and the proposed deep distributional time series model. Sec. (ref) compares the deep time series and benchmark model forecasts, Sec. (ref) considers the inclusion of demand forecast data, and Sec. (ref) concludes. The Web Appendix (“WA” hereafter) provides key algorithms, computational details and additional empirical results.
Wholesale electricity markets include the European Power Exchange, mutiple regional markets in the U.S. (such as the PJM interconnection and the Southwest Power Pool), and national markets in many countries including Australia, Chile and Turkey. While the designs of these markets differ, they are largely “day-ahead” markets where generators and distributors place bids for the sale and purchase of electricity at an intraday resolution up to one day prior to transmission (or “dispatch”). The market is cleared at a wholesale spot price that reflects the marginal cost of supply at each intraday period. Prices also vary at different geographic reference nodes, creating multiple related price series. Markets are overseen by system operators, which match generation with short-term demand forecasts, impose constraints to ensure system stability (i.e. avoid load-shedding or blackouts), and enforce any price caps; see kirschen2018 for an overview of wholesale markets.
Central to the operation of day-ahead markets is the intraday electricity spot price. For market participants accurate short-term forecasts of the price are key to profitability. Because electricity is a flow commodity with a high cost of storage, arbitrage opportunities are limited. This fact, along with the complexities of transmission and that short-run demand is inelastic with respect to price, means that prices exhibit unique stylized characteristics; see knittel2005,karakatsani2008,panagiotelis2008 and weron2014 for summaries of these. From a time series perspective, this includes strong and complex nonlinear serial dependence, while from a distributional perspective prices have very heavy tails, skew and often multiple modes that correspond to different regimes janczura2010 and economic equilibria smith2018 .
Many methods have been used for short-term forecasting of electricity spot prices; see weron2014 and nowotarski2018 for recent overviews of point and probabilistic forecasting methods, respectively. In the machine learning literature, while shallow neural networks (NNs) have long been popular for forecasting electricity prices, DNNs have the potential to produce more accurate forecasts. For example, lago2018 and ugurlu2018 both found that RNN models provide more accurate point forecasts than a range of benchmark models. However, most previous usages of shallow and deep NNs have focused on point forecasts of prices, and only in a few cases also quantify predictive uncertainty using bootstrap or other Monte Carlo methods rafiei2016.
In contrast, a number of other methods have been used to construct probabilistic forecasts misiorek2006,panagiotelis2008,huurman2012,bunn2016. One particularly promising avenue is to extend distributional regression methods to time series forecasting. For example, gianfreda2018 and narajewski2020 do so for German electricity prices using exogenous covariates, and serinaldi2011 does so for Californian and Italian electricity prices using historical prices and other variables. These papers report more accurate probabilistic forecasts. In Sec. (ref) we develop a copula-based model that exploits the accuracy exhibited by deep models within a distributional time series forecasting setting, thereby combining the advantages of both approaches.
Copulas have become increasingly popular in time series models of electricity prices because they can capture complex dependence, while also allowing for highly flexible margins. smith2012, ignatieva2016, manner2016, pircalabu2017 and manner2019 use low-dimensional copulas to capture cross-sectional dependence between regional prices, price spikes and other energy series in interconnected power systems, while smith2018 use high-dimensional copulas to capture both serial and cross-sectional dependence jointly for multiple regional prices. However, these studies all employ copulas that are very different to the copula process proposed here.
Electricity generation in the NEM is an important component of economic activity, with 19.4bn Australian dollars of turnover during the 2018-2019 financial year AERstate. Operations in the NEM are managed by the Australian Electricity Market Operator (AEMO), and since April 2006 it has had five regions which correspond to the power systems in the states of New South Wales (NSW), Queensland (QLD), Victoria (VIC), South Australia (SA) and Tasmania (TAS). Separate prices are set at a central location (or “node”) in each region, although they are dependent because the state-based power systems are interconnected by high voltage direct current lines.
Participating utilities place bids for the purchase (by distributors) and sale (by generators) of electricity at five minute intervals up to one day prior to dispatch. The trading price is the average clearing price for this auction over six consecutive five minute periods, so that it is observed at a half-hourly resolution. Re-bidding of prices is allowed before dispatch, although not the amount of energy. Price forecasting in the NEM has been studied extensively, with contributions by higgs2009, panagiotelis2008, nowotarski2013, ignatieva2016, rafiei2016, smith2018, apergis2019 and manner2019 among others. Intraday forecasts over a horizon of 24 hours are used by market participants to develop effective strategies for bidding, re-bidding and managing risk. In our study, we employ the half-hourly trading prices (measured in Australian dollars per MW/h) observed during 2019 in the five regions of the NEM. Prices can be negative for short periods (when it is more cost-effective to sell into the market at a loss, rather than ramp down generation temporarily) although the floor price is -\$1,000. There is also a maximum price that is adjusted annually on 1 July, which was \$14,500 prior to 1 July 2019 and \$14,700 afterwards.
Fig. A in the WA plots prices in NSW during four weeks corresponding to the four seasons, illustrating the strong heterogeneity based on the time of day, day of the week and season. Table B in the WA provides a summary of the prices in each region. These have extreme positive skew, so that we follow most previous studies and work with the logarithm of price $Y=\log(\mbox{Price}+1001)$. The addition of 1001 accounts for the minimum price, at which $Y=0$. It is straightforward to construct density forecasts of the nominal price from those of $Y$ using the Jacobian of the transformation. Fig. B in the WA plots histograms of the five price series, along with (bounded) kernel density estimates (KDEs), showing that even the distribution of $Y$ is both asymmetric and heavy-tailed for each series. Figure (ref) gives boxplots of $Y$ for NSW (which is the region with the most energy demand), broken down by hour of the day. It shows that the entire distribution varies substantially over the day. Equivalent plots for the day of the week, and month of the year, (see the WA) also reveal similar heterogeneity in the distribution of prices.
RNNs represent complex temporal relationships between variables by allowing cycles and sequences in their hidden layers. However, RNNs are typically computationally expensive to estimate and can be numerically unstable for long series PasMikBen2013, both of which are issues here. One way to address this is to employ a variant called an Echo State Network (ESN) that decreases the number of weights that need to be trained Jae2007. The key idea of ESNs is to only adapt the output layer in training, while keeping the weights of recurrent and input connections fixed yet randomly assigned. Such an approach is called a `reservoir computing' method LukJae2009, because it establishes a multiple linkage hidden reservoir which can be of much higher dimension than the input, providing a lot of flexibility at a lower computational cost. However, classical ESNs and many of its extensions rarely consider uncertainty quantification, and in this section we outline probabilistic ESNs that do so. The approach is extended in Sec. (ref) to allow for marginal calibration of the data distribution using a deep copula construction.
We first specify the ensemble ESN proposed for spatio-temporal data by mcdermott+w17. For series with highly nonlinear dependence, these authors demonstrate that allowing the response equation to depend on the hidden states quadratically in an ESN can increase predictive accuracy, which we also find for the electricity price series examined here. Let $\{Y_t\}$ be a stochastic process, then the ESN with Gaussian disturbances $\varepsilon_t\sim\mathcal{N}(0,\sigma^2)$ takes the following form for $t\geq 1$:
Here, $\text{\boldmath$x$}_t$ is an $n_x$-dimensional feature vector that includes unity for an intercept term, $\text{\boldmath$h$}_t$ is a $n_h$-dimensional hidden state vector, `$\circ 2$' denotes the element-wise square of a matrix, $V,U$ are $(n_h\times n_h)$ and $(n_h\times n_x)$ matrices of hidden layer weights, and $g_h$ is an activation function, set here to $g_h(x)=\tanh(x)$. The parameter $0<\kappa\leq 1$ and is known as the “leaking rate”. Similar to mcdermott+w17 we set $\kappa=1$ after having checked that predictive performance is not improved by setting $\kappa<1$, so that $\tilde \text{\boldmath$h$}_t=\text{\boldmath$h$}_t$. The feature vector $\text{\boldmath$x$}_t$ contains past values of both the response and other series, which we specify in Sec. (ref). The response $Y_t$ is a nonlinear function not only of $\text{\boldmath$x$}_t$, but also of all previous values $\text{\boldmath$x$}_1,\ldots,\text{\boldmath$x$}_{t-1}$, so the stochastic process is not Markov.
The constant $\lambda_V$ is the largest eigenvalue of $V$, and $0<\delta<1$ a scaling parameter, so that $(\delta/\lambda_V)V$ has spectral radius less than one. A spectral radius greater than one can result in unstable behavior in the latent states LukJae2009. We follow mcdermott+w17 and set $\delta=0.35$ and $n_h=120$, although we found the forecasting results to be insensitive to variations in these settings.
The elements of the matrices $V=\{v_{il}\},U=\{u_{ij}\}$ are assumed to be random, and distributed independently from mixtures of a uniform distribution and a point mass at zero. If $\mathcal{U}(a,b)$ denotes a uniform distribution over domain $(a,b)$, $\mathcal{B}(\pi)$ denotes a beta distribution with mean $\pi$, and $\delta_0$ is the Dirac function at zero, then the elements
We follow mcdermott+w17 and set $a_v=a_u=\pi_v=\pi_u=0.1$ after having checked the predictive performance under several alternative settings.
Consider $T$ time series observations $\text{\boldmath$Y$}=(Y_1,\ldots,Y_T)'$ of the stochastic process with corresponding $(T\times n_x)$ matrix of feature values $X=[\text{\boldmath$x$}_1|\cdots|\text{\boldmath$x$}_T]'$. Denote $\text{\boldmath$\xi$}=\lbrace V,U,\kappa,\delta\rbrace$, $H_{\xi}(X)=[\text{\boldmath$h$}_1|\cdots|\text{\boldmath$h$}_T]'$ as the $(T\times n_h)$ matrix of hidden state values, $B_{\xi}(X)=[\text{\boldmath$\iota$},H_{\xi}(X)|H_{\xi}(X)^{\circ 2}]$, $\text{\boldmath$\iota$}$ as a vector ones, and $\text{\boldmath$\beta$}=(\beta_0,\text{\boldmath$\beta$}_1',\text{\boldmath$\beta$}_2')'$. Then (ref) can be written as the linear model
Given $\text{\boldmath$\xi$}$, $X$ and $\text{\boldmath$h$}_0=\bm{0}$, the hidden state matrix $H_{\xi}(X)$ is known without error as the hidden state vectors can be computed recursively. Only $\text{\boldmath$\beta$}$ and $\sigma^2$ (which we refer to as model parameters) require estimation, for which we use their Bayesian posterior distribution. Differing from mcdermott+w17, we regularize $\text{\boldmath$\beta$}$ by adopting the shrinkage prior
where $\mathcal{IG}$ denotes an Inverse Gamma distribution. We found a ridge prior with $P(\tau^2)=\tau^2 I$, and hyper-prior $\tau^2\sim\mathcal{IG}(\tilde a,\tilde b)$ with $a=b=\tilde a=\tilde b=0.001$ to work well. The posterior of the model parameters $\text{\boldmath$\vartheta$}=(\text{\boldmath$\beta$}, \sigma^2,\tau^2)$ of this regularized linear model can be computed using the standard MCMC sampler at Algorithm 1 in Part A of the WA.
Most implementations of ESNs only draw a single set of weights from (ref). Here, we follow mcdermott+w17,mcdermott2019deep and simulate $K=100$ matrices $\{V^{k},U^{k}\,; k=1,\ldots,K\}$ from (ref). An ensemble is then used to integrate over $U,V$ when constructing the probabilistic forecasts. If $\text{\boldmath$\xi$}^{k}=\lbrace V^{k},U^{k},\kappa,\delta\rbrace$, then the density forecast of $Y_{T+h}$ at time $T$ for each point $h=1,\ldots,h_1$ in the forecast horizon is the ensemble
where the subscript notation indicates $f_{T+h|T}$ is conditional on the filtration at time $T$.
The density $p^k$ in (ref) is the Bayesian posterior predictive density computed for configuration $\text{\boldmath$\xi$}^k$ as follows. Let $X_{(t)}\equiv[\text{\boldmath$x$}_1|\cdots|\text{\boldmath$x$}_t]'$, then
which is an integral over any unobserved feature values and the model parameter posterior.\footnote{It is implicit that all densities in the integrand of (ref) are conditional on weight configuration $\text{\boldmath$\xi$}^k$.} From (ref), the first term in the integrand is the density of a $\mathcal{N}(\text{\boldmath$h$}_{T+h}'\text{\boldmath$\beta$}_1+(\text{\boldmath$h$}_{T+h}')^{\circ 2}\text{\boldmath$\beta$}_2,\sigma^2)$ distribution, where $\text{\boldmath$h$}_{T+h}$ is computed through the recursion of the hidden state equation (therefore $\text{\boldmath$h$}_{T+h}$ is a deep function of $\text{\boldmath$x$}_1,\ldots,\text{\boldmath$x$}_{T+h}$). The outer integral in $\text{\boldmath$\vartheta$}$ can be evaluated by averaging over draws from the posterior $p(\text{\boldmath$\vartheta$}|X,\text{\boldmath$y$})$ obtained from running the MCMC sampler. However, plugging in the posterior mean $E(\text{\boldmath$\vartheta$}|X,\text{\boldmath$y$})$ computed from the Monte Carlo sample is much faster, yet can be almost as accurate, so we follow this approach.\footnote{Implementation requires running the MCMC sampler and computing the posterior mean a total of $K$ times, once for each hidden weight configuration $\text{\boldmath$\xi$}^1,\ldots,\text{\boldmath$\xi$}^K$.}
In our empirical work $\text{\boldmath$x$}_t$ contains past values of both the focal price and the other four regional prices, so that $\text{\boldmath$x$}_1,\ldots,\text{\boldmath$x$}_{T+1}$ are observed at or before time $T$. However, some elements of $\text{\boldmath$x$}_{T+2},\ldots,\text{\boldmath$x$}_{T+h}$ are unobserved, and the integrals over these feature vectors in (ref) are with respect to their unobserved elements only. The integrals are computed in a Monte Carlo manner by simulating all five series values sequentially from their predictive distributions. Algorithm 2 in Part B of the WA simulates from the ensemble density (ref) for all five regional price series and over each time point in the forecast horizon.
Adopting Gaussian disturbances to the output layer at (ref) is inconsistent with the strong asymmetry in the empirical distribution of (the logarithm of) prices. Thus, we also consider the skew-$t$ distribution of Azz2003 with location zero for the disturbances, with density $p(\varepsilon_t;\omega^2,\alpha,\nu) = \frac{2}{\omega}t_\nu \left(\frac{\varepsilon_t}{\omega}\right) T_{\nu+1}\left( \frac{\alpha\varepsilon_t}{\omega} \sqrt{\frac{\nu+1}{\nu+x_{\varepsilon_t}^2}}\right)$, where $t_\nu$ and $T_{\nu}$ are the student $t$ density and distribution functions. The parameters $\omega^2,\alpha,\nu$ are scale, skew and degrees of freedom parameters, respectively, and we write $\varepsilon_t\sim\mathcal{ST}(0,\omega^2,\alpha,\nu)$. When $\nu>30$ the density is effectively that of a skew-normal distribution, whereas when $\alpha=0$ it is that of a $t$ distribution. We fix values of $\nu$, and use the prior at (ref) for $\text{\boldmath$\beta$}$, so that the model parameters are $\{\text{\boldmath$\beta$},\tau^2,\omega^2,\alpha\}$.
Bayesian estimation of the model parameters uses a conditionally Gaussian representation of the skew-$t$. This introduces latent variables $\zeta_t,w_t$, and expresses $\varepsilon_t$ conditional on these values as $\varepsilon_t=\psi\zeta_t+\epsilon_t$, where
$\mathcal{G}$ denotes a Gamma distribution, $\mathcal{TN}_{[0,\infty)}$ is a Gaussian truncated to $[0,\infty)$, and $\{\alpha,\omega\}$ are re-parameterized as $\{\psi,\sigma\}$ with $\psi=\alpha\sigma$ and $\sigma^2=\omega^2-\psi^2$. Integrating out $\zeta_t,w_t$ recovers the skew-$t$ distribution $\varepsilon_t\sim\mathcal{ST}(0,\omega^2,\alpha,\nu)$. Following FruPyn2010, the hyperpriors are $\psi\sim\mathcal{N}(0,D_0)$, $\sigma^2\sim\mathcal{IG}(c_0,C_0)$ and $\tau^2\sim\mathcal{IG}(b_0,B_0)$, with constants $D_0=1$, $c_0=2.5$, $b_0=1$, $B_0=0.005$, $C_0=0.5s_y^2$ and $s_{y}^2$ denoting the sample variance of the response $Y_t$.
Given $\text{\boldmath$\xi$}$, the re-parameterized model parameters $\text{\boldmath$\vartheta$}=\{\text{\boldmath$\beta$},\tau^2,\psi,\sigma^2\}$ are estimated using their Bayesian posterior computed by an MCMC sampler that generates $\text{\boldmath$\zeta$}=(\zeta_1,\ldots,\zeta_T)'$ and $\bm{w}=(w_1,\ldots,w_T)'$, similar to those proposed by panagiotelis2008 and FruPyn2010. As with the Gaussian probabilistic ESN, prediction is based on Monte Carlo draws from the ensemble density at (ref) produced using a minor adjustment of Algorithm 2. Details are given in Part B of the WA.
While the deep time series models in Sec. (ref) provide probabilistic forecasts they have two drawbacks: (i) the feature vector only affects the mean of the response equation at (ref), and (ii) the density forecasts are not calibrated in any manner. In this section, a copula model is outlined that is a deep distributional time series model in which the feature vector affects the entire predictive distribution and the probabilistic forecasts are marginally calibrated.
gneiting+br07 discuss different forms of calibration of density forecasts, including predictive marginal calibration, which is defined as follows. For $t=T+1,\ldots,\infty$, assume a future observation of $Y_t$ has true distribution $H_{t|T}(y)$ and forecast distribution $F_{t|T}(y)$, where the subscript indicates that both distributions are conditional on the filtration at time $T$. Then, if $\bar H(y)\equiv\lim_{h_1\rightarrow \infty}\left\{\frac{1}{h_1}\sum_{h=1}^{h_1}H_{T+h|T}(y)\right\}$ and $\bar F(y)\equiv\lim_{h_1\rightarrow \infty}\left\{\frac{1}{h_1}\sum_{h=1}^{h_1}F_{T+h|T}(y)\right\}$, the forecast distributions are called marginally calibrated if and only if $\bar F(y)=\bar H(y)$.
gneiting+br07 highlight that because $H_{t|T}$ is unknown for $t>T$, marginal calibration can be assessed in practice by comparing the empirical distribution function over the forecast horizon to the average of the corresponding distributional forecasts. In a study with a moving window (such as ours) and forecast origins $T=T_{\mbox{\tiny start}},\ldots,T_{\mbox{\tiny finish}}$, there are a total of $N=T_{\mbox{\tiny finish}}-T_{\mbox{\tiny start}}+1$ sets of forecasts. Averaging over these and the forecast horizon of length $h_1$, gives
where $\mathcal{I}(X)$ is an indicator function equal to one if $X$ is true, and zero otherwise. The closer $\hat F$ is to $\hat H$, the greater evidence for predictive marginal calibration.
We adopt a copula model for the joint distribution of $\bm{Y}_{(t)}=(Y_1,\ldots,Y_t)'$, conditional on $X_{(t)}=[\text{\boldmath$x$}_1|\cdots|\text{\boldmath$x$}_t]'$ and weight configuration $\text{\boldmath$\xi$}$, with the density decomposition
where $\text{\boldmath$y$}_{(t)}=(y_1,\ldots,y_t)'$ and $\text{\boldmath$u$}_{(t)}=(u_1,\dots,u_{t})'$. The deep copula process has $t$-dimensional density $c_{\mbox{\tiny ESN}}(\text{\boldmath$u$}_{(t)}|X_{(t)},\text{\boldmath$\xi$})$ specified below in Sec. (ref). The density $p_Y$ and corresponding distribution function $F_Y$ in (ref) are time invariant, and are estimated non-parametrically from the training data. This ensures in-sample marginal calibration, which will extend to predictive calibration whenever the series is stationary.
We stress that even though in (ref) the distribution $Y_t|\text{\boldmath$x$}_t$ is assumed marginally invariant with respect to $\text{\boldmath$x$}_t$, the distribution $\bm{Y}_{(t)}|X_{(t)}$ is still related to the matrix of feature vector values $X_{(t)}$ though the joint distribution. A consequence is that for $t>T$ the predictive density of $Y_t$ can be heavily dependent on the feature vector $\text{\boldmath$x$}_t$, as outlined in Sec. (ref).
At (ref) we employ an (implicit) copula process with density $c_{\mbox{\tiny ESN}}$ constructed from the joint distribution of a second stochastic process $\{\tilde Z_s\}$ that follows the Gaussian probabilistic ESN at (ref) with $\text{\boldmath$\beta$}$ integrated out under the prior (ref). We call $\tilde Z_s$ a “pseudo-response” because it is not observed directly, but is introduced only for specification of its implicit copula. The $t$ observations $\tilde{\bm{Z}}_{(t)}=(\tilde Z_1,\ldots,\tilde Z_t)'$, are conditionally distributed
with $\tilde{\text{\boldmath$z$}}_{(t)}=(\tilde z_1,\ldots,\tilde z_t)'$ and the matrix $B_{\xi}(X_{(t)})$ is specified at (ref), but without an intercept term (i.e. the first column of ones) because level is unidentified in a copula.
The implicit copula of a Gaussian distribution is called a Gaussian copula, and is constructed for (ref) by standardizing the distribution. Let $\bm{Z}_{(t)}=(Z_1,\ldots,Z_t)'=\sigma^{-1}S\tilde{\bm{Z}}_{(t)}$, where $S=\mbox{diag}(\psi_1,\ldots,\psi_t)$ is a diagonal scaling matrix with elements $\psi_s=(1+ \text{\boldmath$b$}_s' \text{\boldmath$b$}_s/\tau^2)^{-1/2}$, and $\text{\boldmath$b$}_s$ is the $s$-th row of $B_{\xi}(X_{(t)})$. Then $\bm{Z}_{(t)}|X_{(t)},\sigma^2,\tau^2,\text{\boldmath$\xi$}\sim \mathcal{N}(\bm{0},R)$ where
has ones on the leading diagonal, and is a function of $(X_{(t)},\text{\boldmath$\xi$},\tau^2)$. Then the copula has density $c_{\mbox{\tiny ESN}}(\text{\boldmath$u$}_{(t)}|X_{(t)},\text{\boldmath$\xi$},\tau^2) = \frac{\phi(\text{\boldmath$z$}_{(t)};\bm{0},R)}{\prod_{s=1}^t \phi_1(z_s)}\,,\label{eq:copdens}$ where $z_t=\Phi_1^{-1}(u_t)$, $\text{\boldmath$z$}_{(t)}=(z_1,\ldots,z_t)'$, $\phi_1$ is the standard normal density, and $\phi(\cdot\,;\bm{0},R)$ is the density of a $\mathcal{N}(\bm{0},R)$ distribution. Notice that $\sigma^2$ does not feature in $c_{\mbox{\tiny ESN}}$, which is because implicit copulas are scale free.
Because $c_{\mbox{\tiny ESN}}$ is conditional on $X_{(t)}$, this is a Gaussian copula process on the feature space Wilson2010. Because the copula captures the dependence structure in $\{Z_t\}$, $\bm{x}_t$ contains past values of this process (rather than $\{Y_t\}$), and the equivalent processes for the other four regional electricity prices. These can be computed easily as $Z_t=\Phi^{-1}(F_Y(Y_t))$. We write the copula density $c_{\mbox{\tiny ESN}}$ also as a function of the unknown parameter $\tau^2$ and employ it in the copula model at (ref).
Given a configuration $\text{\boldmath$\xi$}$, the only unknown copula parameter is $\tau^2$, for which we adopt the Weibull prior of KleKne2016 with scale parameter $b_{\tau^2}=2.5$. Direct estimation using the likelihood at (ref) is difficult because evaluation of $c_{\mbox{\tiny ESN}}(\text{\boldmath$u$}|X,\text{\boldmath$\xi$},\tau^2)$ requires inversion of the matrix $R$, which is computationally infeasible for all but small sample sizes. This problem is solved by expressing the likelihood conditional on $\text{\boldmath$\beta$}$ as follows. For a sample of size $T$, denote the observations as $\text{\boldmath$y$}\equiv \text{\boldmath$y$}_{(T)}$ and feature matrix as $X\equiv X_{(T)}$. Then the conditional likelihood is obtained by a change of variables from $\bm{y}$ to $\bm{z}=(z_1,\ldots,z_T)'$, with elements $z_t=\Phi_1^{-1}(F_Y(y_t))$, so that \[ p(\text{\boldmath$y$}|X,\text{\boldmath$\beta$},\text{\boldmath$\xi$},\tau^2)= p(\text{\boldmath$z$}|X,\text{\boldmath$\beta$},\text{\boldmath$\xi$},\tau^2)\prod_{t=1}^T \frac{p_{Y}(y_t)}{\phi_1(z_t)}\\ =\phi\left(\text{\boldmath$z$};\bm{0},SB_\text{\boldmath$\xi$}(X)\text{\boldmath$\beta$}, S^2\right)\prod_{t=1}^T\frac{p_{Y}(y_t)}{\phi_1(z_t)}\,, \] which can be evaluated in $O(T)$ operations because $S$ is diagonal. Part C of the WA provides an MCMC sampler to generate draws from the augmented posterior of $(\text{\boldmath$\beta$},\tau^2)$, so that $\text{\boldmath$\beta$}$ is integrated out in a Monte Carlo manner and direct computation of $R$ is avoided.
As with the probabilistic ESNs in Sec. (ref), an ensemble is used to integrate over the distribution of $U,V$, so that the density forecast of $Y_{T+h}$ is again given by (ref). The ensemble components $p^k$ (i.e. the Bayesian posterior predictive densities at (ref)) are derived from the copula model where $\text{\boldmath$\vartheta$}=\{\text{\boldmath$\beta$},\tau^2\}$. To do so, the first term in the integrand of (ref) is obtained by a change of variables from $y_{T+h}$ to $z_{T+h}=\Phi_1^{-1}(F_Y(y_{T+h}))$, so that
where $\mu_{T+h}=\frac{1}{\psi_{T+h}}\text{\boldmath$b$}_{T+h}\text{\boldmath$\beta$}$, $\psi_{T+h}=(1+\frac{1}{\tau^2}\text{\boldmath$b$}_{T+h} \text{\boldmath$b$}_{T+h}')^{-1/2}$ and $\text{\boldmath$b$}_{T+h}=(\text{\boldmath$h$}_{T+h}',(\text{\boldmath$h$}_{T+h}^2)')$ is a row vector. Notice that the entire density at (ref) is a nonlinear function of the feature vector $\text{\boldmath$x$}_{T+h}$ via $\text{\boldmath$b$}_{T+h}$, so that $Y_{T+h}$ is not marginally invariant of $\text{\boldmath$x$}_{T+h}$ in the predictive distribution. As in Sec. (ref), the posterior mean $E(\text{\boldmath$\vartheta$}|\text{\boldmath$y$})$ is computed from the Monte Carlo sample and plugged in for $\text{\boldmath$\beta$},\tau^2$ in (ref). Because $\text{\boldmath$x$}_t$ comprises past values of $\{Z_t\}$ for both the focal price, and the other four regional prices, the integrals over $\text{\boldmath$x$}_{T+2},\ldots,\text{\boldmath$x$}_{T+h}$ in (ref) are computed by simulating values for $Z_{T+t}$ sequentially for all five regions. Algorithm 4 in Part C of the WA generates iterates from the ensemble density for the copula model.
We illustrate the accuracy of our deep model forecasts using half-hourly price data observed during 2019 in all five regions. Competing models were fit using a moving window of three months training data, starting with the period 1 February--30 April. (We start at February, rather than January, because the input vector $\bm{x}_t$ defined below includes lagged values of the series.) This window was advanced monthly until the period 1 September--30 November, resulting in eight fits per model. Forecasts were constructed for a eight month evaluation period from 1 May to 31 December; a total of 245 days. The origin $T$ was advanced every half hour, and probabilistic forecasts constructed for a horizon of 24 hours (i.e. $h_1=48$ half-hours).\footnote{Because there are 245 days in the evaluation period, this means there was a total of $245 \times 48$ forecast origins. At each forecast origin, a total of $48 \times 5=240$ probabilistic forecast densities were constructed for each model. Thus, there were a total of 2,822,400 probabilistic forecast densities for each model.} When constructing forecasts during each month, the model parameter estimates obtained from the three preceding months training data were used. The reason each model was only refit eight times (rather than for every one of 11,760 half-hourly forecast origins) was to reduce the computational burden of the study. A forecast horizon up to 24 hours is necessary for the bidding and re-bidding process in the NEM.
The following four deep time series models were applied to each of the five price series:
For the RNN, RNNST and RNNSN, the feature vector contains the $n_x=271$ elements \[ \text{\boldmath$x$}_t=(1,\underbrace{\bm{Y}_{\mbox{\tiny{ALL}},t-1}, \bm{Y}_{\mbox{\tiny{ALL}},t-2},\ldots, \bm{Y}_{\mbox{\tiny{ALL}},t-48}}_{\text{prices in previous 24 hours}}, \underbrace{\bm{Y}_{\mbox{\tiny{ALL}},t-96},\bm{Y}_{\mbox{\tiny{ALL}},t-144}, \ldots,\bm{Y}_{\mbox{\tiny{ALL}},t-336}}_{\text{prices at same hour 2 to 7 days previously}})'\,, \] where $\bm{Y}_{\mbox{\tiny{ALL}},t}= (Y_{\mbox{\tiny{NSW}},t}, Y_{\mbox{\tiny{QLD}},t}, Y_{\mbox{\tiny{SA}},t}, Y_{\mbox{\tiny{TAS}},t}, Y_{\mbox{\tiny{VIC}},t})$ is the vector of all five regional prices at hour $t$ and the first entry is the intercept. For RNNC, the $\text{\boldmath$x$}_t$ contains values of the transformed series $Z_t=\Phi_1^{-1}(F_Y(Y_t))$ for each of the five regions at the same lagged time points as above. The choice of inputs is motivated by previous studies panagiotelis2008,higgs2009 that identify strong dependence in Australian electricity prices with those in the previous 24 hours, and at the same time during the previous seven days. Price dependence is induced by serial dependence in both supply-side factors and in electricity consumption smith2000. Complex cross-sectional serial dependence in price is also well-documented panagiotelis2008,higgs2009,ignatieva2016,han2020 and is caused by inter-regional trade in electricity between the five regions.
While the deep time series models above are trained separately for each price series, the inclusion of lagged prices from all regions as elements of $\text{\boldmath$x$}_t$ forms five-dimensional multivariate deep times series models. From each model multi-step ahead forecast densities for the five regional prices are constructed jointly via simulation using Algorithms 2 and 5 outlined in the WA. At each forecast origin, $2000$ Monte Carlo draws from (ref) are obtained for each half hour $h=1,\dots,48$ and region $r=1,\ldots,5$.
The density $p_Y$ was estimated for each series using the bounded KDE implemented in MATLAB function `ksdensity', with bounds set equal to the admissible prices during the forecast period. These were a floor price of -\$1000 ($Y_{T+h}\geq 0$), and a maximum price of \$14,500 prior to 1 July ($Y_{t+h}\leq 9.582$) and \$14,700 on or after that date ($Y_{t+h}\leq 9.596$). These densities are given in Fig. B in the WA. This ensures that the predictive densities from the copula model at (ref) are constrained to the same range, which is important for the accurate forecasting of prices near the bounds. In contrast, density forecasts from the other models are unconstrained, and can have substantial mass outside the bounds. Thus, we simply truncate their predictive densities to the admissible price region and normalize.
In addition to our four proposed deep time series models, following serinaldi2011 and gianfreda2018 we also employ GAMLSS as a benchmark. In this framework the response $Y_t$ has a specified parametric distribution, where the parameters are modeled as flexible functions of covariates. We used the R package `gamlss', but found the software unstable when the ESN output layer terms were used as covariates, and also for many of the candidate distributions supported by the package. However, the four parameter Johnson's SU distribution with $\text{\boldmath$x$}_t$ as linear covariates worked well, and we include forecasts from this model. serinaldi2011 also found the JSU distribution to work well for electricity prices, which are a challenging series with which to calibrate models due to the frequent price spikes. Details on our GAMLSS implementation and experiments are given in Part D of the WA.
For point forecasts, the mean absolute error (MAE) and root mean square error (RMSE) are used to measure accuracy. For a univariate probabilistic forecast distribution $F$ there are a growing number of measures of accuracy, many of which are listed by nowotarski2018, and we compute the following four. The first is the quantile score function $\mbox{QS}_\alpha(F^{-1}(\alpha),y)=2({\cal I}(y<F^{-1}(\alpha))-\alpha)(\alpha-y)$, where ${\cal I}(A)=1$ if $A$ is true, and zero otherwise. This is also called the “pinball loss” raiffa61,gneiting2011. While the entire quantile score function is of interest, we also consider its value at $\alpha=0.05$ and $\alpha=0.95$ as measures of lower and upper tail accuracy, respectively. The second measure is the continuous ranked probability score (CRPS) gneiting+br07 defined as $\mbox{CRPS}(F,y) =\int_0^1 \mbox{QS}_\alpha(F^{-1}(\alpha),y)d\alpha$ for observed value $y$, where the integral is computed numerically. It is a proper scoring rule which measures overall probabilistic forecast accuracy.
The third measure is a loss function proposed by fissler2016 that is based on the “value-at-risk” and “expected shortfall”, which are popular measures of financial tail risk for low values (i.e. financial losses). To measure tail risk associated with high values (which is the primary source of price risk in electricity markets) we employ the restatement of the loss function for the upper tail given by nolde2017. For $Y\sim F$ and $\alpha$ close to 1, let $\mbox{EL}(\alpha)\equiv E(Y|Y\geq F^{-1}(\alpha))$ be the “expected longrise” (analogous to the expected shortfall measure of lower tail risk). Then from nolde2017, the loss function is
where $G_1(x)=x$ and $G_2(x)=G_3(x)=\exp(x)$. We compute this score for $\alpha=0.975$, and lower values indicate more accurate upper tail forecasts of the price distribution. The fourth and last measure is the empirical coverage of the predictive distributions at the 95% level.
For each of the first three metrics, we report the weighted average of their values over the evaluation period and over regions, with weights equal to the proportion of total electricity consumption in each region.\footnote{These weights are 0.3687 (NSW), 0.2355 (VIC), 0.2818 (QLD), 0.0624 (SA) and 0.0516 (TAS).} We call these “system-weighted” average metric values, and they reflect the differing importance of price forecasting accuracy across regions in the NEM. For the fourth metric, the coverage is with respect to the price distributions in all five states.
Fig. (ref) presents histograms of the out-of-sample observations of $Y$ for each region during the evaluation period. Also plotted are the average predictive densities $\hat F$ for the five forecasting methods, allowing a visual comparison of out-of-sample marginal calibration. The copula model (RNNC) produces forecasts that are most accurately calibrated, followed by RNNSN and RNNST. However, GAMLSS and RNN both exhibit poor predictive marginal calibration. The exceptionally wide intervals for RNN are caused by extremely inaccurate density forecasts at horizons close to 24 hours ahead as characterized further below.
To visualize the predictive densities, Fig. (ref) plots the 95% prediction intervals for forecasts made during the month of July 2019. Results from the five methods for forecasts 1, 6 and 24 hours ahead are given in separate panels. RNN forecasts are negatively biased, and for some days lack “sharpness” gneiting+br07. The latter is because a high variance is necessary to capture the price spikes in the training sample using a distribution with symmetric thin tails for the errors at (ref). In contrast, skew-normal disturbances (RNNSN) can produce a heavier upper tail through right skew, but at the expense of overly sharp densities. De-coupling of the level of asymmetry from the kurtosis through skew-$t$ disturbances (RNNST) allows for asymmetric predictive densities with heavy tails. RNN, RNNSN and RNNST all have probabilistic forecasts that are nearly homoscedastic, because they are ensembles of homoscedastic Bayesian predictive densities. In contrast RNNC allows for a much heavier upper tail for electricity prices and heteroscedasticity; the latter is because it is a distributional model. Last, the performance of the benchmark GAMLSS model (which is also a distributional model) declines substantially over the forecast horizon.
To assess forecast accuracy we compute the metrics outlined in Sec. (ref) for each of the five methods. Table (ref) reports out-of-sample point forecast accuracy, and the asymmetric deep time series models RNNSN, RNNST and RNNC are more accurate than RNN and GAMLSS in terms of MAE at all horizons, but not by the RMSE metric at short horizons. At horizons of 12 or more hours, RNNC produces the most accurate point forecasts. However, the main objective is accurate probabilistic forecasting, and Table (ref) reports out-of-sample metrics for this. RNN is the poorest method by all measures, which is unsurprising given the asymmetry and heavy tails in even logarithmic prices. GAMLSS is also very poor overall. RNNST is more accurate than RNNSN, demonstrating that allowing for heavy tails increases accuracy over-and-above allowing for asymmetry. However, it is the deep copula model RNNC that provides clearly superior forecasts at all points in the forecast horizon, both in terms of the entire density (as measured using CRPS) and the upper and lower tails (as measured by the quantile scores and upper tail loss function).
Fig. F in the WA plots the full quantile score functions (pinball loss). The poor calibration of the RNN and GAMLSS probabilistic forecasts is clear, as is the superiority of the asymmetric deep time series models (in particular RNNC). To compare the probabilistic forecasts for the two most accurate models (RNNC and RNNST) Fig. (ref) plots the expected shortfall ($E(Y|Y<F^{-1}(\alpha))$ for $\alpha=0.025$) and expected longrise ($E(Y|Y>F^{-1}(\alpha))$ for $\alpha=0.975$) at 12:00 during July 2019 for both models. The RNNC produces predictive forecasts that exhibit substantial heteroscedasticity and are sharper (i.e. lower spread between the longrise and shortfall) than the RNNST. These differences in the two deep time series models are because RNNC is a distributional time series model, whereas RNNST is not. Finally, Table (ref) reports the coverage of the predictive distributions at the 95% level, and there is a substantial difference between the methods, with the predictive distributions from RNNC exhibiting superior coverage at all horizons.
We now show how to extend the deep time series models to include demand forecasts as additional information in a flexible fashion, allowing for greater price forecast accuracy.
Demand for electricity is almost perfectly inelastic to price in the short term because individual users face fixed tariffs. As demand varies over time, the spot price traces out the average supply curve. This explains the occurrence of price “spikes” because electricity supply curves are typically kinked, and for short periods of time demand can exceed the location of this kink; see geman2006,clements2015,smith2018 and ziel2016. An implication is that including accurate demand forecasts in a price forecasting model may increase its accuracy, particularly for the upper tail.
However, incorporating demand forecasts into most existing time series models for price is difficult for two reasons. First, multiple aspects of probabilistic demand forecasts beyond the point forecast are likely to increase the accuracy of probabilistic price forecasts. Second, demand forecasts are likely to be related to future prices in a highly nonlinear fashion, particularly when combined with past price values. However, both difficulties can be addressed by incorporating multiple summaries of the demand forecast distribution as additional elements of the feature vector of a deep time series model.
Methods for accurate short term probabilistic forecasting of demand are well established, with examples for NSW demand provided by smith2000 and cottet2003. AEMO makes demand forecasts publicly available, which are updated every half hour by an automated system. The forecasts are of the 10th, 50th (i.e. median) and 90th percentiles of demand at a half-hourly resolution over a horizon of one week. We employ demand forecasts for NSW that are exactly 24 hours ahead, to forecast NSW prices 24 hours ahead. Thus, only demand forecasts truly available at the forecast origin $T$ were used. We label the three percentile forecasts of NSW demand as $\mbox{D}10,\mbox{D}50$ and $\mbox{D}90$. Fig. D in the Web App. plots actual demand against $\mbox{D}50$ during 2019, showing the high degree of point forecast accuracy. However, empirical coverage during 2019 suggests some minor miscalibration, with 12.1%, 62.5% and 95.6% of demand observations falling below $\mbox{D}10,\mbox{D}50$ and $\mbox{D}90$, respectively.
The forecasting study was extended, where the NSW demand forecasts were included in the feature vectors of our two best performing deep time series models, RNNST and RNNC, for forecasting NSW price. For RNNST \[ \text{\boldmath$x$}_t=(1,Y_{\mbox{\tiny{NSW}},t-1}, Y_{\mbox{\tiny{NSW}},t-2},\ldots,Y_{\mbox{\tiny{NSW}},t-48}, Y_{\mbox{\tiny{NSW}},t-96},Y_{\mbox{\tiny{NSW}},t-144}, \ldots,Y_{\mbox{\tiny{NSW}},t-336}, \mbox{D}10_t,\mbox{D}50_t,\mbox{D}90_t)'\,, \] so that $n_x=58$. The feature vector for RNNC contained the equivalent lagged values of the transformed series $Z_{\mbox{\tiny{NSW}},t}=\Phi_1^{-1}(F_Y(Y_{\mbox{\tiny{NSW}},t}))$, and the three demand forecasts. For comparison, the two models without the demand forecasts (so $n_x=55$) are also used.
Table (ref) reports the average value of the forecast metrics over the evaluation period. Including the demand forecasts as features decreases the accuracy of the point forecasts and the CRPS metric for both models. However, upper tail forecast accuracy increased significantly for the RNNC. This is important in this application because short periods of very high prices often dominate the profitability of participants in the NEM. To illustrate, Figure (ref) plots the density forecast where the inclusion of the demand forecast information greatly increased accuracy (other examples are given in Fig. G of the WA).
Finally, we also experimented with employing the ESN at (ref) without the quadratic term in the output layer, but found that this did decrease forecast accuracy significantly. Moreover, including ridge priors with different values of $\tau^2$ for the linear and quadratic terms in the output layer had little effect on forecast accuracy; see Part F of the WA for further details.
This paper makes contributions to both the time series and energy economics literatures. We propose new deep time series models that exploit the reservoir computing techniques found in ESNs for high frequency time series. Recasting a quadratic ESN as a statistical time series model, with a shrinkage prior for regularization of the output layer weights, allows for uncertainty quantification and probabilistic forecasting. Our first approach is to allow for different error distributions for the output layer, extending the Bayesian methodology of mcdermott+w17,mcdermott2019deep. However, our main methodological contribution is the proposal of a new deep {\em distributional} time series model. This is obtained by constructing the implicit copula of a Gaussian probabilistic ESN, and extends the deep distributional regression method of KleNotSmi2019 to time series. This copula is a deep process on the feature space, allowing for highly adaptive nonlinear serial dependence, and generalizing existing echo state Gaussian processes chatzis2011echo. When combined with a nonparametric estimate of $p_Y$, it allows for marginal calibration. The entire density forecast is a function of the feature vector through (ref), unlike with the other deep time series models. Our empirical work also suggests that accurate in-sample marginal calibration also results in more accurate calibration of both the out-of-sample marginal and predictive distributions.
While the distribution of intraday electricity prices have a sizable predictable component, these time series are complex ignatieva2016,manner2019. Key to their accurate modeling and forecasting is to capture jointly three features: nonlinear serial dependence, high levels of asymmetry and kurtosis, and strong time-variation in the distribution. All our proposed deep time series models account for the first feature, while allowing for skew-normal or skew-$t$ disturbances also accounts for the second. However, only the deep copula model allows for all three, producing the most accurate probabilistic forecasts. Recent models ziel2016,shah2020 suggest that including demand forecasts may further increase probabilistic forecasting accuracy. We find this to be the case for the upper tail of the 24 hour ahead NSW electricity prices when including demand forecasts as features in the deep copula model.
Last, we highlight two areas for future work. First, estimation of our deep time series models using variational methods, rather than MCMC, has the potential to speed up computations. Second, there are other applications of our deep time series models. One example is macroeconomic forecasting, where large regularized time-varying parameter models are popular bitto2019,carriero2019,huber2020. The deep time series models provide a flexible alternative, as they can incorporate feature vectors of many economic variables regularized through the reservoir structure of the ESN, plus Bayesian shrinkage of the output layer coefficients. Moreover, the copula model can produce time-varying and asymmetric probabilistic forecasts, which is a feature of this problem. \singlespacing