Econometric Modeling of Regional Electricity Spot Prices in the Australian Market

1. Melbourne Business School, University of Melbourne, 200 Leicester Street, Carlton, VIC 3053, Australia; Email: mike.smith@mbs.edu
2. Department of Information, Risk, and Operations Management, McCombs School of Business, University of Texas at Austin, 1 University Station, B6500, Austin, TX 78712-0212, U.S.A; Email: Tom.Shively@mccombs.utexas.edu

Corresponding Author

This work was partially supported by Australian Research Council Future Fellowship FT110100729.

Econometric Modeling of Regional Electricity Spot Prices in the Australian Market

Abstract

Wholesale electricity markets are increasingly integrated via high voltage interconnectors, and inter-regional trade in electricity is growing. To model this, we consider a spatial equilibrium model of price formation, where constraints on inter-regional flows result in three distinct equilibria in prices. We use this to motivate an econometric model for the distribution of observed electricity spot prices that captures many of their unique empirical characteristics. The econometric model features supply and inter-regional trade cost functions, which are estimated using Bayesian monotonic regression smoothing methodology. A copula multivariate time series model is employed to capture additional dependence — both cross-sectional and serial — in regional prices. The marginal distributions are nonparametric, with means given by the regression means. The model has the advantage of preserving the heavy right-hand tail in the predictive densities of price. We fit the model to half-hourly spot price data in the five interconnected regions of the Australian national electricity market. The fitted model is then used to measure how both supply and price shocks in one region are transmitted to the distribution of prices in all regions in subsequent periods. Finally, to validate our econometric model, we show that prices forecast using the proposed model compare favorably with those from some benchmark alternatives.

Key Words: Bayesian Monotonic Function Estimation, Intraday Electricity Prices, Copula Time Series Model.

JEL: C11, C14, C32, C53.

## 1 Introduction

During the past two decades, traditional vertically integrated electricity power systems have been replaced with wholesale markets in many countries and regions. While the design of such markets varies, day-ahead and bid-based markets are common. In such markets, bids are placed by generators, distributing utilities and third parties at an intraday resolution one day prior to consumption (or ‘dispatch’). These auctions are overseen by network management organizations, which schedule the dispatch of electricity on a least cost basis, subject to demand forecasts and system security requirements. The spot price is the highest priced bid dispatched to meet demand at a point in time. The effective modeling and forecasting of spot prices is important for utilities to operate in a market profitably, and also for the management of a market; see Kirschen & Strbac (2004) for an introduction to contemporary electricity markets.

In North America, Australia and New Zealand, markets also feature nodal pricing, where there are different prices in different regions within a market. Electricity is traded between regions via high voltage interconnectors, and increased synchronization of regional prices is often considered symptomatic of increased market efficiency. Moreover, even when nodal pricing is not a market feature, the construction of interconnectors between adjacent power systems allows for trade between markets with different prices. In this study, we propose a new approach for the modeling of the joint distribution of regional electricity spot prices in a bid-based interconnected market. A statistical model is proposed that is motivated by an economic spatial equilibrium pricing model. We aim to show that exploiting the structural relationships between regional supply, interconnector flows and regional prices suggested by the economic model helps explain some of the unique empirical features of electricity prices. It also provides a framework from which density forecasts can be constructed and event studies undertaken.

Electricity is a flow commodity that cannot be stored economically. Demand has a strong predictable component — for example, see Harvey & Koopman (1993), Ramanathan et al. (1997), Smith (2000), Taylor et al. (2006), Clements et al. (2016) and many others — and is almost perfectly inelastic to the wholesale price in the short run because most consumers face fixed tariffs. Supply comprises generators with very different marginal costs of production and ‘ramping’ times (i.e. the time required to change the amount generated). Trade between regions is also constrained by interconnector capacity. These features ensure that time series of spot prices exhibit unique characteristics, including sharp spikes, asymmetry, periodicity, and cross-sectional and serial dependence; all of which have been documented widely (Knittel & Roberts 2005; Panagiotelis & Smith 2008; Karakatsani & Bunn 2008b; Clements et al. 2017). In response to this complexity, a variety of nonlinear statistical models have been developed for modeling and forecasting intraday spot prices; see Weron (2014) for an extensive overview. Of these, regime switching regression and time series models have been particularly successful (Huisman & Mahieu 2003; Weron et al. 2004; Haldrup & Neilsen 2006; Karakatsani & Bunn 2008b; Janczura & Weron 2010; Bunn et al. 2016), where different latent regimes correspond to different price distributions under different economic equilibria. Yet less attention has been given to the multivariate modeling of regional prices, which is what we undertake here.

Our approach is motivated by an extension of the spatial equilibrium model of DeVany & Walls (1999). This is a network economic model (Nagurney 1999) that relates each regional supply curve (often called a ‘stack’ in the power systems literature) and inter-regional transmission cost functions, to prices in all regions. It features three inter-regional price equilibria. The first is where prices are synchronized across regions, up to transmission costs, the second is where prices deviate due to temporal ramping constraints, and the third is where prices deviate due to interconnector capacity constraints; the latter scenario is often called ‘congestion pricing’ in the power systems literature (Bakirtzis 2001). However, this economic model cannot be implemented directly in practice. One reason is that it is written in terms of directed electricity flows between all pairwise combinations of regions, which are unobservable in wholesale pools. Another reason is that aggregate cost and supply functions are difficult to construct for many markets.

We therefore consider a statistical approximation, where pairwise flows are replaced with observed directed flows on inter-regional interconnectors, and cost and supply functions are estimated using the Bayesian monotonic function estimation methodology of Shively et al. (2009). We extend the approach of these authors to allow the regression disturbances to follow a mixture of three normals, with moments corresponding to those of the price distributions under the three theoretical equilibria. To allow for additional serial and cross-sectional dependence in prices we employ a multivariate time series model based on a copula construction of the type discussed by Biller (2009), Smith (2015) and Smith & Vahey (2016). Such dependence is likely to arise due to strategic bidding, the impact of intermittent renewable generation, the effect of climatic conditions on the power system, and other market imperfections. The marginal mean of each price series is determined by the monotonic regression means, and the mean-corrected prices are modeled using their empirical distribution functions, so that each price series distribution is marginally nonparametric. Both serial and cross-sectional dependence is then captured by a high-dimensional Gaussian copula, with parameter matrix equal to the autocorrelation matrix of a stationary vector autoregression (VAR) process. The VAR has a long lag length, but with many null coefficient matrices, so that it is still parsimonious. The autocorrelation matrix of the VAR can be readily estimated using sparse matrix computations, and provides an estimate of the Gaussian copula parameter matrix.

Copulas are popular tools for multivariate modeling because they allow dependence to be captured separately from the other features of the distribution, which are in turn captured by arbitrary marginal models. They are used increasingly in the energy modeling literature, with applications as diverse as accounting for spatial dependence in renewable generation (Papaefthymiou & Kurowicka 2009), hedging of electricity futures (Liu et al. 2010) and modeling energy spreads (Westgaard, 2014; Ch.4). Copulas are also attractive for the multivariate modeling of regional electricity prices, because they allow for extreme levels of asymmetry and heavy tails through an appropriate choice of marginal model. For example, Smith et al. (2012), Wang, Cai & He (2015), and Ignatieva & Trück (2016) all use various low-dimensional copulas to capture cross-sectional dependence between electricity prices in different regions, while Manner, Türk & Eichler (2016) use low-dimensional copulas to capture the cross-sectional dependence between price spike incidence in different regions. However, our use of copulas differs from these authors in that we use a much higher dimensional copula to capture both serial and cross-sectional dependence. Such a copula-based time series model preserves the nonparametric marginal distributions of prices, including the high levels of asymmetry, heavy tails and regression function means.

We outline how to construct the joint predictive density of regional prices from the fitted model using forecasts of electricity load, which are readily available in practice; for example, see Clements et al. (2016) and references therein. The network economic model motivates a different multivariate time series model for each supply region, each producing a separate joint predictive distribution of prices in all regions. We combine these distributions to produce an ensemble forecast distribution. Ensemble methods for combining predictive densities from different models have proven both popular and useful in the fields of meteorology (Sloughter et al. 2010) and macroeconometrics (Mitchell & Hall 2005; Jore et al. 2010).

We apply our model to the Australian National Electricity Market (NEM) using half-hourly data. This market is the world’s longest interconnected power system, with five regions that each have their own price setting mechanisms. The market design shares features with many other wholesale markets, and the price series exhibit the empirical characteristics often observed in such markets. Estimates of the supply functions are highly nonlinear, and in line with those from previous studies; for example, see Geman & Roncoroni (2006). We show how the fitted models can be used to assess the impact on prices in all regions of a supply-side shock in just one region, such as the loss of a major generating unit. We also compute the generalized impulse response (Koop, Pesaren & Potter 1996) of the nonlinear multivariate time series model to assess how regional price shocks are transmitted to the distribution of prices in all regions in subsequent periods. In addition, to validate our model we show that its predictions are more accurate than two naïve benchmarks.

We note here that our approach has a number of features that are similar to those found in some other recent statistical models of electricity prices. These include quantile regression and time series models (Bunn et al. 2016; Jónsson et al. 2014), which are popular for this problem because they produce flexible forecast densities. By employing nonparametric marginal distributions, this is also the case for the forecast densities produced by our copula multivariate time series model. However, it is considerably more complex to extend quantile regression to the multiple time series case we focus on here. Regression or time series models that include ‘fundamental’ supply and/or demand-side variables have also proven popular for forecasting electricity prices; see, for example, Karakatsani & Bunn (2008a; 2008b), Ferkingstad et al. (2011) and Gonzalez et al. (2012). Our regression models with monotonic effects in supply and inter-regional flows also employ fundamental supply-side variables, motivated by the structural relationships in the network economic model. Moreover, by exploiting an economic model, our study is similar in spirit to that of Ziel and Steinert (2016), who construct dynamic demand and supply curves based on a time series model fit using auction data. Last, constructing ensemble forecast densities from five separate multivariate statistical time series models for regional prices, mirrors the long-standing practice of point forecast combination used previously in the energy literature; for example, see De Menezes et al. (2000).

The rest of the paper is organized as follows. Section 2 outlines the network economic model, and the statistical model for prices. Section 3 introduces the NEM and the data. Section 4 discusses estimation of the monotonic regressions, and the copula model. The density forecasting methods are also discussed here. Section 6 contains the empirical analysis, including the two event studies and the validation study, while Section 7 concludes.

## 2 Modeling Electricity Prices

### 2.1 Spatial Equilibrium Pricing

We follow DeVany & Walls (1999) and consider a spatial equilibrium model of price formation. We denote the independent price-setting regions in an inter-connected electricity market as nodes . Let be the set of origin-destination (O/D) pairs, between which there is an energy flow of from node to node . In general, the flow is not directly observable in an electricity network, because it consists of the sum of flows down all different pathways from node to node in the power system. Nevertheless, the electricity demanded at node , and the electricity supplied at node , are both observed. Under the assumption of zero transmission loss, the feasibility conditions are

 dj=∑i∈ΩQi,j and bi=∑j∈ΩQi,j,

and the total energy generated and consumed is equal with .

In the absence of load-shedding, demand for electricity in most bid-based markets is almost perfectly inelastic with respect to price at any instant in time (Kirschen 2003). This is because wholesale prices are not passed through directly to the consumer, who largely face fixed tariffs. We denote the monotonically increasing inverse supply function at node as . In a perfect market, the following spatial equilibrium conditions (Nagurney 1999, Chap. 3.1) hold to clear the market for all O/D pairs :

 Si(bi)+Ci,j(Qi,j){=πjifQi,j>0≥πjifQi,j=0, (2.1)

where is the demand price at node , and is the cost of flow . The monotonically increasing cost function measures the wheeling costs due to transmission losses, congestion in the transmission system and possibly aspects of market design. In the above, we refer to the origin node as the supply region.

Flows between nodes are constrained by capacity constraints on the interconnectors, so that there is also a constraint on flow . In the presence of this constraint, the spatial equilibrium conditions (Nagurney 1999, Chap 3.3) at Equation (2.1) are now

 Si(bi)+Ci,j(Qi,j)⎧⎪ ⎪⎨⎪ ⎪⎩=πjif0

for all O/D pairs. The first of the two inequalities occurs when there are no gains from trade because the supply price plus wheeling costs at node , , exceeds the price received at node . For example, this can occur in practice when baseline generators are supplying electricity at node at a low price because there are additional costs to ramping down the generators; see Wolak (2007) for a discussion of the impact of ramping costs on market clearing conditions. The second inequality occurs when trade is constrained because of capacity constraints on inter-regional interconnectors. In practice, in many bid-based markets, this can result in a price spike at node , because to meet demand at this node electricity has to be supplied locally, and the inverse supply function typically has a first derivative that is monotonically increasing. For a depiction of a stylized inverse supply function, see Geman & Roncoroni (2006).

### 2.2 Statistical Model of Regional Prices

In bid-based markets, electricity prices are set at equally-spaced points in time at an intra-day resolution. It is difficult to employ the spatial equilibrium at Equation (2.2) to directly model these prices for a number of reasons. First, as flows between nodes increase, the energisation of an electricity network is governed by Kirchoff’s voltage law (Wood & Wollenberg 1995; Krischen & Strbac 2004). This determines how the currents (and therefore the power flows) distribute themselves through a network, in a manner that is not necessarily efficient economically. Second, in a centralized bid-based wholesale market, while flows are observable on individual interconnectors, flows between individual node pairs cannot be distinguished.111Note that can be observed in theory in markets with bilateral trades, although these are rare in practice, and do not occur in the Australian NEM. Third, the upper bounds also depend on any other flows that share transmission pathways. Fourth, extension of the static model at Equation (2.2) to a dynamic situation is difficult, in part due to differing ramping times and costs for different generation capacity.

Therefore, we instead use the spatial equilibrium model in Section 2.1 to motivate a statistical model for electricity prices. Each supply region produces a separate multivariate model for prices observed at all price-setting nodes in the network. We outline the model for a given supply region below, and it has two main components. The first is a series of regression models for prices, each of which has a finite mixture distribution for the disturbances. The second is a copula to capture both cross-sectional and serial dependence in prices.

Let be a set of directed arcs that connect the O/D pair that corresponds to the physical interconnectors between the two nodes. If there are no such interconnectors, then . Let be the flow on directed arc , and be a monotonically increasing transmission cost function on the bounded interval , with . In our statistical model, the cost function for the unobservable flow is replaced by the sum of the costs of flows along the arcs in , so that . At time , let be the price in region , the supply in region , and the flow on directed arc . Then, the equilibrium conditions at Equation (2.2) motivate the following regression for observations at times :

 πj,t = Si(bi,t)+∑a∈Ei,jca(va,t)+ϵi,j,t = ηi,j,t+ϵi,j,t.

There are separate regressions for each combination , so that there are regressions in total. We refer to the regression above as the region supply regression for prices in region . The random variable is marginally distributed as a mixture of three distributions, with distribution function

 Fϵi,j(ϵ)=3∑l=1ωi,j,lGl(ϵ;αi,j,l,σi,j,l). (2.4)

Here, the weights are such that , and is a distribution function with mean and standard deviation . The marginal mean is therefore nonzero, and the marginal mean of prices is .

The adoption of this mixture distribution is motivated by the three equilibria in the spatial model. It is also consistent with previous empirical analyses that detect two or three regimes in time series of prices; for example, see Karakasani & Bunn (2008b) and Janczura & Weron (2010). To identify the first two moments of the three mixture components, we assume the following parameter constraints:

• Component 1: and are unconstrained (baseline case);

• Component 2: and  (lower mean and higher variance);

• Component 3: and  (higher mean and higher variance).

Component 1 corresponds to the case where in Equation (2.2). Component 2 corresponds to the case where , which can occur when high costs to ramping down large baseline generators produce lower (even negative) prices with higher variability in comparison to Component 1. Component 3 corresponds to the case where , and prices are both higher and more volatile in comparison to Component 1. This last component includes situations that lead to price spikes; something that is observed empirically in many bid-based electricity markets, including the Australian market. While these moment constraints are motivated by the features of the three equilibria, they also identify the likelihood so that the problem of label-switching (e.g. see Frühwirth-Schnatter 2006; Sec. 3.5) does not occur in their estimation.

A multivariate model for the impact of supply in region on prices in all regions over times can be constructed using a copula function. Copulas are popular tools to account for cross-sectional dependence in non-Gaussian models (Patton 2006), and also serial dependence in time series (Smith et al. 2010, Smith 2015). They allow the adoption of arbitrary marginal models, followed by the modeling of dependence in a second separate step. Here, we employ the regression models above as the marginal models for the elements of . We employ a Gaussian copula (Song 2000) to capture both cross-sectional and serial dependence in prices over-and-above that explained by supply and interconnector flows, as discussed below.

For each supply region , we use a copula representation of the distribution of all observed prices . That is, where the joint distribution function is , with a -dimensional copula function; for example, see Nelsen (2006; p.43). The copula function is evaluated at the transformed data vector , where , is uniformly distributed on , and . The values are often called probability integral transformed (PIT) observations in the time series literature (Rosenblatt 1952), or ‘copula data’ in the copula literature when computed using estimates of and . We later show in Section 4.2 how such copula data can be used to estimate the copula function . The density function of is obtained by differentiating its distribution function as Here, is also a density function on , commonly called the ‘copula density’. Each marginal density is a mixture of three components , with .

If denotes the standard normal distribution function, then a Gaussian copula has the density

 k{\tiny Ga}(ξ;Ω)=|Ω|−1/2exp{−12w′(Ω−1−In)w}, (2.5)

for a vector with , and a correlation matrix as dependence parameters. It is easy to show that is distributed , and that all sub-vectors of also have Gaussian copula densities; for example, see Song (2000). Here, we assume that corresponds to the dimensional correlation matrix of a stationary multivariate time series . In particular, we follow Biller (2009) and Smith & Vahey (2016) and employ a stationary Gaussian VAR() for the latent process. The matrix is a block Toeplitz autocorrelation matrix, with th block (Lütkepohl 2006; p.30). The parameter matrix captures cross-sectional linear correlation in the latent time series process, and captures serial linear correlation at lags . However, when combined with highly non-Gaussian margins, the copula framework also allows for nonlinear dependence in the price series.

When employing a Gaussian copula we set the copula density . In this case, it is possible to show that the multivariate series of the PIT values is a (strongly) stationary time series on the unit cube with Markov order (Smith 2015). Similarly, the series is also strongly stationary, where and . Measures of dependence for this time series can be computed from the matrices of the Gaussian copula, and forecast distributions constructed via simulation. We discuss computation of these, along with estimation the Gaussian copula model, further in Section 4.2.

## 3 Australian National Electricity Market

### 3.1 The Market

The Australian National Electricity Market began operating in December 1998, although regional markets were adopted several years earlier. It consists of regions, which coincide with the five adjacent Australian states of New South Wales (NSW), Victoria (VIC), Queensland (QLD), South Australia (SA) and Tasmania (TAS); with the latter region joining in May 2005. All sales of electricity go through a wholesale pool, which is managed by the Australian Electricity Management Organization (AEMO). While some limited generation and consumption of electricity does occur independently at remote locations disconnected from the transmission grid, there are no bilateral trades utilizing the grid. In 2011 there were 305 registered generators (Australian Energy Regulator 2011). While generation capacity based in VIC and SA was almost entirely privately owned during this period, approximately 90% and 70% of capacity in NSW and QLD, respectively, was owned by public corporations; along with all capacity in TAS. Coal and gas fired generators made up the vast bulk of available capacity, with some hydroelectric capacity in TAS and on the NSW/VIC border, and a minor wind and solar capacity located mainly in SA. Table 1 provides a breakdown of the registered capacity during 2011 by region and source.

AEMO operates a separate price setting mechanism in each region. Generators bid for the supply of electricity into the pool one day ahead. Each bid consists of 48 half-hourly prices in Australian dollars per megawatt hour ($/MWh) and quantities in megawatt hours (MWh). Bids are ordered by price, and the highest marginal price of generation capacity dispatched (i.e. employed to meet demand) is computed for every five minute interval. The half-hourly spot price is equal to the average of these six five minute prices, and is the price at which all sales are made during that half-hour. Re-bidding of amount (but not price) is allowed up to five minutes before dispatch, and is widely practised. Prior to 1 July 2010 there was a price cap of$10,000 per MWh, after which it was increased to $12,500. There is a floor price of -$1,000. Negative prices occur in the NEM for short periods of time, either as a result of the high cost of ramping down baseline coal generators, or due to strategic bidding behavior where participants are able to exploit transmission constraints (AER 2011). For an introduction to the price setting mechanism see AEMO (2010).

There is extensive inter-regional trade using six high voltage interconnectors. Figure 1 shows these, and the twelve corresponding directional flow variables, which have O/D pairs as listed in Table 2. From these the sets of directed arcs for this network can be derived; for example, and . Note that of the twenty sets of directed arcs between the five regions, twelve are empty sets. AEMO transmits electricity between regions, within interconnector capacity constraints, with the objective of equalizing prices up to the cost of transmission. When the interconnectors are at capacity, prices in each region often differ substantially.

### 3.2 The Data

To estimate our statistical model we employ half-hourly observations on regional loads, spot prices, and interconnector flows and losses, made between 7 February 2010 and 13 February 2011. We use the publicly available AEMO dispatch data222We are careful to employ the dispatch, rather than pre-dispatch, data because it is from this that market prices are computed. from the website www.aemo.com.au. Data on the supply variable is computed for region from the relationship

 bi= load in region i+xi−mi+interconnector loss % adjustment. (3.1)

Here, are total exports from region to other regions, are total imports into region . For example, using the labels for interconnector flows given in Table 2, exports from NSW are and imports into NSW are . Load in each region is the total demand for electricity, and includes any distribution losses within the region.333Our load variable is labelled ‘TotalDemand’ in the dispatch dataset. This does not include any ‘normally-off’ local scheduled loads, which are often exactly zero and are typically excluded from the definition of demand by AMEO. The interconnector loss adjustment arises from inter-regional transmission line losses.444Our loss adjustment variable for each region is the ‘Allocated Interconnector Losses’ for each region in the dispatch dataset.

#### 6.3.2 Price Impulse

The impact of price shocks not caused by changes in supply, such as unexpected changes to demand or bidding behavior, are captured by the multivariate time series . The impact of such a shock can be measured by comparing the predictive distribution of prices for subsequent periods with, and without the shock. This is a popular approach, and for nonlinear multivariate time series analysis, it is often called a generalized impulse response analysis; for example, see Koop, Pesaran & Potter (1996). For supply region , let be the value of the time series at time with th element , where is the price shock in region on the logarithmic scale.101010For example, if this is a 200 $/MWh increase at time , then the increase on the logarithmic scale is . Then we compute the forecast distribution of with and without the shock over a horizon of half-hours ahead. This is computed by simulating Monte Carlo iterates from the copula multivariate time series model for future values . To illustrate, we consider the impact of a$200 increase in prices (ie. an ‘impulse’) in VIC on prices in all regions using the fitted copula model with margins given by the VIC supply regressions. We consider the impulse over the two hour period 03:00–04:59 on 19 May 2010, denoted as half-hours and . Figure 8 plots estimates of the predictive distributions of , both without the shock (blue lines) and with the shock (red lines). The panels are arranged as a matrix, with rows corresponding to prices in different regions, and columns corresponding to predictive distributions and 96 half-hours ahead. Prices are on the logarithmic scale, and the density estimates are computed using a kernel method with locally adaptive bandwidth (Shimazaki & Shinomoto 2010). All predictive distributions have a long right tail, which has been omitted here for presentation purposes. The non-smooth nature of the presented densities is not due to either Monte Carlo sample error,111111We use 50,000 Monte Carlo iterates to construct each density highly accurately. nor due to inefficiencies in the kernel density estimates. Instead, it is because in the marginal models the nonparametric empirical distribution functions were employed as part of the model, and these are non-smooth.121212Alternatively, smooth predictive densities can be readily obtained by employing smooth parametric models — such as regressions with skew t disturbances — for the margins.

The response to the price shock is strongest at shorter horizons and in the same region as the price impulse (ie. in VIC). The bottom row of panels show the large impact on the predictive distribution of price in VIC, with the impulse increasing the predictive mean by 868.53, 438.54, 403.02, 519.33 and 317.05 $/MWh at horizons and 96 half-hours ahead, respectively. Prices in QLD are almost completely unaffected by the impulse in VIC, which is consistent with the low pairwise dependence between prices in these two regions depicted in Figure 6. While it initially appears in Figure 8 that the impulse has limited effect on prices in NSW, SA and TAS, that is not the case. The impulse substantially accentuates the right hand tail of the predictive distributions in these regions, making price spikes more likely in the forecast horizon and increases the means. To illustrate, the price impulse in VIC increases mean predictive prices one day () ahead in NSW, SA and TAS by 5.12, 26.03 and 123.46$/MWh, respectively.

### 6.4 Validation Study

To validate our fitted model we undertake a limited forecasting study. A number of different models were fit to hourly data and used to forecast hourly prices.131313Hourly, rather than half-hourly, data are used to reduce the computational burden of the validation study. However, the models fitted to hourly and half-hourly data are very similar, so that the results are unlikely to differ meaningfully. We construct forecasts with a daily expanding window with 100 forecast origins from 24 October 2010 to 31 January 2011. From each origin, we construct forecasts of the five regional prices over a horizon of one week at the hourly resolution.141414In total, we produce forecasts for each of the five regional prices at hours, so that a total of separate price forecasts are made. Forecasts were obtained from each of the following methods:

• Naïve 1: The price at the same time on the last observed day is used as a benchmark point forecast.

• Naïve 2: The average price at the same time of the day over the training sample is used as a second benchmark point forecast.

• Fundamental: The ensemble of the marginal expectations from our regression models, , is used as a point forecast. This is computed assuming the values of supply and flows in the forecast period are known, so that , where the marginal expectation comprises constants defined in Section 2.2. This approach treats the five price series as both serially and cross-sectionally independent, conditional upon supply and flow observations.

• Copula & Fundamental 1: The ensemble distribution at Equation (5.2), with its mean used as a point forecast. This forecast distribution is computed assuming the values of supply and flows in the forecast period are known. The approach differs from that above in that it exploits any serial and cross-sectional dependence captured by the copula model.

• Copula & Fundamental 2: The same ensemble distribution and point forecast as above, but with supply and flows in the forecast period determined via hour-by-hour optimizations as outlined in Section 5.2.

• Copula: The copula multivariate time series model, but without the regressions outlined in Equation (2.2). Instead, the marginal distribution of each price series is modeled only by its empirical distribution function, so that no structural information is exploited.

To summarize the overall level of accuracy, we consider forecasts of the demand-weighted (log) price across all five regions,

 π\tiny DWt=5∑l=1(dl,t∑rj=1dj,t)πl,t,

forecasts of which are constructed from the five regional forecasts . The point forecast accuracy is measured using the mean absolute forecast error (MAFE). The means are computed over all forecasts, broken down by different ranges of the forecast horizon, and reported in Table 5. For example, the MAFE at horizons of 4, 5 and 6 hours ahead is 0.983 for the ‘Copula’ method, and is an average of AFE values.151515The corresponding measures of accuracy of the individual regional price forecasts, and also for very high prices (which we define to be prices above the 95th percentile), are reported in Tables 1–6 of the Supplementary Material.

A number of insights can be drawn from the results. First, the two naïve approaches are dominated in most circumstances by the statistical methods, highlighting the value of a model-based approach. Second, the extension of the ‘Fundamental’ method to include multivariate serial dependence via the copula model (ie. the ‘Copula & Fundamental 1’ method), improves the accuracy up to a horizon of 2 days ahead. Third, for horizons of 3 or more days, the forecasts from the ‘Fundamental’ method dominate all alternatives, highlighting the potential value of incorporating the structural information in forecasts at a long horizon. Fourth, forecasting interconnector flows (and supply) using optimization (ie. the ‘Copula & Fundamental 2’ method), reduces price forecast accuracy. Nevertheless, the forecasts dominate the two naïve approaches for short horizons of 1 and 2 hours. The limited accuracy of this approach suggests that it would be worthwhile considering alternative approaches to forecasting interconnector flows. Fifth, the multivariate copula model without structural information (ie. the ‘Copula’ method) performs particularly well for horizons of up to 2 days, but is dominated over longer horizons by the ‘Fundamental’ method that incorporates such structural information. The high accuracy of the copula multivariate time series model at short horizons here mirrors that documented by Smith & Vahey (2016) for macroeconomic variables. Last, we note that only point forecast accuracy is considered here. Given the asymmetry in the logarithm of prices, it is worthwhile to extend the study to consider density forecasts from the statistical models and their accuracy.

## 7 Discussion

In this paper we use a spatial equilibrium model of price formation to motivate multivariate statistical models for a vector of regional electricity prices. There is a separate multivariate model for each supply region, and the estimates and predictions from each model can be combined into an ensemble. Key features of each multivariate statistical model include marginal regressions with monotonic functions, and a high-dimensional Gaussian copula to capture additional serial and cross-sectional dependence. In estimating each regression model we employ a finite mixture of Gaussians for the disturbances, which is motivated by the three price equilibria in the economic model. However, a nonparametric disturbance using an infinite mixture model of the type popular in the Bayesian literature (Hjort et al. 2010) could also be used. When estimating the copula models, we use a two stage estimator. First, the marginal distributions of the mean-corrected prices are estimated using their empirical distribution functions. Second, the copula parameters are estimated using maximum likelihood. Joint estimation of a Gaussian copula and nonparametric marginals is far from straightforward (e.g. see Rosen & Thompson 2015), and would be even more difficult joint with the monotonic functions. For this reason two-stage estimators are the most common approach for copula models with complex margins (Joe 2005).

Our econometric model of electricity prices has two main novel features. First, it is a multivariate model that incorporates structural relationships of inter-regional price formation. A number of recent statistical models of electricity prices employ fundamental variables; for example, see Karakatsani & Bunn (2008a) and González et al. (2012). However, our model is the first study of which we are aware where the fundamental variables and their functional forms are motivated by a network economic model. The second novel aspect is that we employ a copula model for additional dependence in electricity prices. Smith, Gan & Kohn (2012) and Ignatieva & Trück (2016) recently employ a 5-dimensional copula to capture the dependence between prices in the NEM, while Manner, Türk & Eichler (2016) employ a 4-dimensional copula to capture the dependence between price spike incidence in four regions in the NEM. However, these low-dimensional copulas capture only cross-sectional dependence, whereas we employ a parsimonious -dimensional copula to capture both cross-sectional and serial dependence, which offers substantial advantages here. These include the ability to model the marginal distributions flexibly, so that the forecast densities are substantially more realistic than those produced from overly simplistic parametric distributions, such as the log-normal. While quantile regression or time series models (e.g. Bunn et al. 2016) also produce similarly flexible forecast densities, they are harder to extend to multiple price series than the copula model. And while we focus on the Gaussian copula here, it is possible to further extend our analysis to vine copulas as in Smith (2015).

Electricity prices in the Australian market are highly dependent across regions, but do not follow the law of one price. This is in part due to the unique physical aspects of energy generation, such as interconnector constraints and differing ramping times for generators; for example, see Clements (2017) for a discussion of the effect of interconnector constraints on prices in the NEM, and Higgs, Lien & Worthington (2015) for a discussion on the role of production capacity and generation mix. However, prices are also likely to be affected by strategic bidding and the exercise of market power by utilities; for example, see theoretical work by Bunn & Gianfreda (2010) and a recent empirical analysis by Apergis, Baruník & Lau (2017). Our empirical findings are consistent with these observations, identifying both strong cross-sectional and serial dependence in prices over-and-above that induced by supply-side relationships in the marginal regressions. We illustrate the usefulness of our model in two event studies. The first is a supply-side shock equivalent to the outage of a major generator, while the second is a price impulse in one region. We find that such shocks have an impact on prices in all regions in the NEM over a horizon of up to one week. This highlights the importance of employing a multivariate model in studies of regional prices in wholesale markets with nodal pricing. Our validation study illustrates that our model also shows promise when used to forecast prices. In comparison to simple benchmarks, the fundamental model that accounts for supply-side factors provides improved forecasts at longer horizons, while the inclusion of the time series copula model increases accuracy at horizons under two days.

Finally, we note two avenues for future research. First, the extent of inter-regional dependence between the variance, skewness and kurtosis of prices, along with the incidence of extreme prices (ie. price spikes), has attracted much current interest; for example, see Lindström & Regland (2012), Aderounmu & Wolff (2014), Manner, Trük & Eichler (2016), Apergis, Baruník & Lau (2017) and Han, Kordzakhia & Trück (2017). Smith & Vahey (2016) show that the Gaussian copula model employed here can produce predictive distributions with time-varying variance, skewness, kurtosis and tail probabilities. The extent to which our proposed model captures inter-regional dependence in these moments and tail probabilities is an interesting question. Second, while we document the accuracy of the point forecasts in a validation study, we do not document the density forecasts from our statistical models. Given the asymmetry in the (log) prices, measuring the accuracy of these is also of great practical interest.

Appendix

This appendix provides further details on the monotonic smoothing method in the presence of a mixture of normals outlined in Section 4.1. Dropping the regression subscripts and for notational convenience, let , , , and . Also, to keep the notation manageable, we will assume there are no cost functions in the model. The cost functions are estimated similarly to (and joint with) the supply function .

To estimate , we follow Shively et al. (2009) and define to consist of the columns of basis vectors that correspond to the nonzero elements of . Therefore, the spline in Equation (4.1) evaluated at the observed values is . To make the model analytically tractable for use in an MCMC algorithm, we re-parameterize to give , where , , and is the lower triangular matrix defined in Section 4.1. Note that the constrained prior for the discussed in Section 4.1 induces a prior for constrained to the region , where .

The posterior distribution for has density

 f(M,ω,α,σ2,J,γ|π)∝f(π|M,α,σ2,J,γ)f(M,ω,α,σ2,J,γ)

where

 f(π|M,α,σ2,J,γ)=3∏l=1{∏t:Mt=l(2πσ2l)−1/2exp[−12σ2l(πt−wJtγJ−αl)2]},

is the t-th row of , and

 f(M,ω,α,σ2,J,γ)=Pr(M|ω)f(ω)f(α)f(σ2)f(γ|J)Pr(J)

is the prior distribution defined in Section 4.1. Letting

 s(π,M,α,σ2,J,γ)=−log[f(π|M,α,σ2,J,γ)],

and following Shively et al. (2011), we introduce a scalar latent variable such that

 f(M,ω,α,σ2,J,γ,z|π)∝e−zI(z>s(π,M,α,σ2,J,γ))f(M,ω,α,σ2,J,γ). (A1)

The sampling scheme below is used to generate Monte Carlo iterates from the posterior distribution, augmented with this latent variable. For notational purposes, let , , , , and represent , , , and , respectively, without the -th element.

Step 0: Start with some initial values of , ;

Step 1: Generate conditional on , ;

Step 2: For , generate conditional on , , , , , ;

Step 3: Generate conditional on ;

Step 4: For , generate conditional on , , , , ;

Step 5: For , generate