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 787120212, 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 interregional trade in electricity is growing. To model this, we consider a spatial equilibrium model of price formation, where constraints on interregional 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 interregional 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 crosssectional 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 righthand tail in the predictive densities of price. We fit the model to halfhourly 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, dayahead and bidbased 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 bidbased 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 crosssectional 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 interregional transmission cost functions, to prices in all regions. It features three interregional 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 interregional 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 crosssectional 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 meancorrected prices are modeled using their empirical distribution functions, so that each price series distribution is marginally nonparametric. Both serial and crosssectional dependence is then captured by a highdimensional 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 lowdimensional copulas to capture crosssectional dependence between electricity prices in different regions, while Manner, Türk & Eichler (2016) use lowdimensional copulas to capture the crosssectional 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 crosssectional dependence. Such a copulabased 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 halfhourly 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 supplyside 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 demandside 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 interregional flows also employ fundamental supplyside 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 longstanding 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 pricesetting regions in an interconnected electricity market as nodes . Let be the set of origindestination (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
and the total energy generated and consumed is equal with .
In the absence of loadshedding, demand for electricity in most bidbased 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 :
(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
(2.2) 
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 interregional interconnectors. In practice, in many bidbased 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 bidbased markets, electricity prices are set at equallyspaced points in time at an intraday 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 bidbased wholesale market, while flows are observable on individual interconnectors, flows between individual node pairs cannot be distinguished.^{1}^{1}1Note 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 pricesetting 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 crosssectional 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 :
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
(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 bidbased 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 labelswitching (e.g. see FrühwirthSchnatter 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 crosssectional dependence in nonGaussian 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 crosssectional and serial dependence in prices overandabove 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
(2.5) 
for a vector with , and a correlation matrix as dependence parameters. It is easy to show that is distributed , and that all subvectors 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 crosssectional linear correlation in the latent time series process, and captures serial linear correlation at lags . However, when combined with highly nonGaussian 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 halfhourly 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 halfhourly 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 halfhour. Rebidding 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 interregional 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 halfhourly 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 data^{2}^{2}2We are careful to employ the dispatch, rather than predispatch, 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
(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.^{3}^{3}3Our load variable is labelled ‘TotalDemand’ in the dispatch dataset. This does not include any ‘normallyoff’ 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 interregional transmission line losses.^{4}^{4}4Our loss adjustment variable for each region is the ‘Allocated Interconnector Losses’ for each region in the dispatch dataset.
Table 3 provides summaries of halfhourly price and supply for the five regions, broken down into periods of peak (09:00–20:00) and offpeak (20:30–08:30) demand. Mean prices during peak periods are approximately double those during offpeak periods, except for prices in TAS where the price difference is less. This is likely because TAS is the only region where supply is dominated by hydroelectric capacity, which often acts to smooth prices. Price spikes occur during periods of transmission congestion and peak demand. However, inaccurate demand forecasts, unanticipated outages, and strategic bidding can also cause prices to spike or fall heavily into negative territory. To illustrate, Figure 2 plots VIC regional prices on a logarithm scale. There were 32 observations of extreme prices above 500$/MWh, and 14 of negative prices. The time series modeling of these, and other stylized empirical features, has been discussed widely in the economics, engineering and forecasting literatures; for example, see Karakatsani & Bunn (2008b), Panagiotelis & Smith (2008), González et al. (2012), Weron (2014), Manner et al. (2016) and references therein. Figure 3 provides pairwise scatterplots of the five price series on the logarithmic scale. While prices are positively dependent, there are frequent substantial deviations between prices in different regions at all price levels.
Table 2 reports the mean and maximum halfhourly directed energy transmissions for the six interconnectors. Both QLD and VIC are major exporters of electricity, while NSW and SA are predominantly importers. The nominal capacity constraints of the directed flows are also given, although actual constraints vary with direction of flow and over time due to variation in local network thermal ratings and voltage/reactive power limits.
4 Estimation
Aggregate supply curves are typically difficult to construct directly from publicly available data, including for the Australian NEM. Moreover, the cost functions are unknown in general because they are a feature of the econometric model that accounts for multiple aspects of the cost of interregional transmission. Therefore, we estimate both the supply and cost functions in Equation (2.2) using the Bayesian monotonic semiparametric regression methodology of Shively et al. (2009). These authors demonstrate the efficiency of their approach for Gaussian disturbances, and we extend it here to the case where the disturbances are distributed as a mixture of three Gaussians. The fitted regressions form marginal models, conditional on which the Gaussian copula model is estimated using maximum likelihood. Such a twostage approach is used widely in the copula literature and is only slightly less efficient than full maximum likelihood (Joe 2005), yet much simpler to compute. We outline both estimation stages in the following two subsections, although we discuss the Markov chain Monte Carlo sampling scheme used in the first estimation stage in the Appendix.
4.1 Monotonic Regression Models
We first consider the specification of the function in Equation (2.2). Without loss of generality, we normalize observations on the covariate to , and then approximate the function using the quadratic regression spline
(4.1) 
Here, are fixed ‘knots’ placed along the domain of the independent variable , such that and . We set , which is large enough to allow for a high degree of flexibility. However, unrestricted estimation of the coefficients results in a function estimate that has high local variance; ie. is nonsmooth. We therefore follow the popular approach of placing a point mass probability at zero on the coefficients, and estimate the function as a Bayesian model average; e.g. see Smith & Kohn (1996). We define if , and if , and assume these values are a priori independent with , for . We set in our empirical work, but note that the results are insensitive to a wide range of values. Given these priors, the final function estimate is obtained through model averaging over .
Monotonicity of is ensured using linear constraints on the coefficients imposed through the prior. Following Shively et al. (2009), let consist of the elements of corresponding to those elements of that are nonzero. Then linear restrictions on the elements of required to impose monotonicity can be written as , where is a lower triangular matrix that depends on and the knots, and means each element of the vector is nonnegative. For example, if , then the three linear constraints , and ensure that is monotonically increasing. Given and , the prior for is a distribution, constrained to the region . Following Shively et al. (2009), we set and . The functions , , in Equation (2.2) are modeled similarly.
To account for the mixture distribution for in Equation (2.2), we follow the common Bayesian approach (FrühwirthSchnatter 2006) of introducing latent indicators, where if observation for regression is from mixture component . Thus, , and is the posterior probability that observation is from component for regression . To simplify the development of the sampling scheme in the Appendix, we assume each component in Equation (2.4) is Gaussian, so that the marginal distribution of is a mixture of three Gaussians. Although we note that when computing the copula data, we later employ nonparametric estimates of the marginal distributions.
For the mixture weights we assume a Dirichlet prior, where , so that a priori and . The component means , constrained to the region as discussed in Section 2.2. We set , which is uninformative relative to the scale of the data. Last, the component variances have a uniform prior distribution on , constrained so that and as in Section 2.2, and with .
The posterior mean is used as the point estimate of , along with similar estimates for the functions , , all of which are computed using the Monte Carlo iterates from the posterior. Posterior means are also used as point estimates for the parameters , and .
4.2 Copula Model
Estimation of the copula model uses the copula data. To compute these values, we first employ the marginal models fitted with the assumption that the disturbances follow a mixture of three normals, and set
where parameter values with hats denote posterior mean point estimates. In our empirical work we find that for some marginal models , the values are close to uniformly distributed, suggesting a mixture of three normals provides a good fit to the regression disturbances. However, for other regressions they deviate meaningfully from a uniform distribution, so that to use these values as copula data would provide poor estimates of the copula parameters. Therefore, we construct the empirical distribution function from the values , and compute the copula data as . For each supply region , provides a nonparametric estimate of the marginal distribution function in Equation (2.4). The approach of employing nonparametric estimates for marginal distributions, followed by a parametric copula, is popular in multivariate modeling; for example, see Shih & Louis (1995).
One advantage of using a Gaussian copula for is that, conditional on the copula data, estimation of the copula parameter matrix by maximum likelihood is straightforward using the VAR representation. Let and , then the MLE of the autoregressive parameters of a stationary VAR() for the time series are readily computed; for example, see Mauricio (1995).^{5}^{5}5The stable VAR models for the series are estimated in Matlab using the routine ‘vgxvar’. From these, point estimates of the autocovariance matrices can be computed as in Lütkepohl (2006; pp.2831), along with estimates of the corresponding autocorrelation matrices , where the diagonal matrix . The matrices are the component blocks of , and the resulting estimate is its MLE; for example, see Song (2000). We note here that Smith & Vahey (2016) propose an alternative Bayesian estimator for the copula parameters of this model based on a drawable vine representation of the Gaussian copula density (Bedford and Cooke 2002) when the latent Gaussian process has Markov order 4. However, in our empirical analysis of halfhourly electricity prices, a long Markov order of is used, so that a vine decomposition involves far too many paircopula terms to be evaluated in practice. In comparison, computation of the MLE of this copula time series model using the VAR representation is computationally feasible using existing software.
Another advantage of using the Gaussian copula is that it is straightforward to obtain inference on the meancorrected time series . This includes measures of dependence and the predictive distributions, as we now outline. If is the th element of , then the pairwise dependence between and , can be measured using Kendall’s tau, which for a Gaussian copula is . These values can be arranged into matrices , which have th elements . Then, measures crosssectional dependence of the vector , and measures serial dependence at lags . We note that is symmetric for , but asymmetric for . We call these matrices ‘autodependence’ matrices, as they are direct generalizations of the autocorrelation matrices in linear multivariate time series, and we compute some examples in our later empirical work. Similar matrices can also be constructed from other pairwise measures of dependence, such as Spearman correlations.
The step ahead predictive distribution with conditional density
(4.2) 
can be evaluated in a Monte Carlo fashion as follows. Iterates of the latent time series steps ahead are generated from the fitted VAR() model,^{6}^{6}6This is undertaken in Matlab using the routine ‘vgxsim’. conditioning on the previous values (which are computed from as discussed above). The elements of each iterate are then transformed as , and further transformed as . These successive transformations produce iterates of that are distributed with the density given above in Equation (4.2).
5 Density Forecasting
Predictive distributions of prices are constructed from the statistical model in two ways. The first is conditional on observed values of supply and interconnector flows. However, while AMEO is likely to have accurate forecasts for these values, they will not be known to market participants in general. Therefore, we also outline a second approach where supply and interconnector flows are also forecast. In either case, because there is a separate statistical model for each supply region , there are separate forecast distributions. We then combine these distributions to produce an ensemble forecast distribution.
5.1 Conditional on supply and interconnector flows
We first consider the case when supply and interconnector flows are assumed known during the forecast window. Let be a vector of the flows along all interconnectors at time . Then for supply region , the predictive density of the price vector steps ahead from forecast origin is
(5.1) 
This forecast is conditional on information at time (ie. the filtration ) and also the future supply and interconnector flows . In Equation (5.1) the value can be computed from the fitted marginal models by plugging in and , while the predictive distribution of can be computed via simulation from the fitted copula time series model as discussed in Section 4.2.
Because there is a separate predictive density for each supply region , we combine these into an ensemble predictive distribution with density
(5.2) 
The weights should not be confused with the weights of the mixture distribution at Equation (2.4). They are assumed equallyvalued (ie. ) in our empirical work for simplicity, although other approaches for determining these can also be employed (Jore et al. 2010). Each component in the ensemble has predictive mean
which can be computed from the fitted model. The weighted sum of these predictive means is used as a point forecast.
5.2 Joint with supply and interconnector flows
To predict supply, we first predict load in each region at an intraday resolution. Forecasting intraday electricity load is a wellstudied problem, and there are a number of effective solutions which can be employed here; see the overview by Weron & Misiorek (2008). For region , at time we denote load as , total imports as and total exports as . Then, from Equation (3.1), supply in region at time is .^{7}^{7}7For simplicity, in our later validation study we use actual demand and allocated transmission losses, but note that forecasts for both are usually readily available to system operators and can be used instead.
Therefore, given load forecast , the problem becomes one of forecasting the future transmission flows from which and can be computed. To undertake this we solve an optimization problem that minimizes differences between expected prices in each region, where the expectations are obtained from our statistical model. For the Australian NEM, this approach mirrors the objectives of interregional transmission scheduling by the network management organization AEMO. There is a separate optimization problem for each component of the ensemble distribution, with an objective function , where
(5.3) 
The values are weights, which we set as proportional to the demand in regions and at time , so that . This increases the contribution of price differences between regions with higher demand, and downweights that for regions with lower demand. The objective function is minimized with respect to , given system constraints on interconnector flows.
For our NEM data we incorporate two sets of system constraints. The first set of constraints is the upper bound on the capacity of the interconnectors. Table 2 outlines nominal upper capacity on the interconnectors in the NEM. However, actual upper capacity differs at any given point in time depending upon a number of factors, so that the observed maximum flows provide more realistic estimates of upper bounds, and we use these as the first set of constraints. The second set of constraints derive from our construction of the directed flow variables from net bidirectional flows along each interconnector. For example, and are directed flows along the NSWVIC Interconnect, with only when (and vice versa). This corresponds to the equality constraint ; similarly for the other directed flow pairings depicted in Figure 1.
The optimization is repeated for all periods in the forecast horizon and each supply region , producing forecast values of interconnector flows, from which supply values are also computed. Forecasting of prices then proceeds as outlined in Section 5.1.
6 Empirical Analysis
6.1 Monotonic Regressions
We employ the logarithm of spot prices as the price data . Because prices are occasionally negative, for each price series the logarithm is computed after subtracting the minimum pool price of $1000, and adding $1. We note that even after a logarithmic transformation, prices are right skewed and have a heavy right tail; which is typical of wholesale electricity prices generally. The Bayesian method is employed to estimate each of the 25 regressions, and posterior means of the monotonic functions and model parameters are computed from the Monte Carlo samples.
For example, Table 4 provides the parameter estimates of the mixture distribution components in each of the five VIC supply regressions. We focus on these regression results because VIC plays a key role as the central region of the interconnected NEM; for example, see the recent discussion by Han, Kordzakhia & Trück (2017). Component 1 has low price variation and occurs 75% of the time for NSW prices, 86% of the time for QLD prices, 90% of the time for SA prices, 90% of the time for TAS prices, and 93% of the time for VIC prices. For price formation in NSW, QLD, TAS and VIC, Component 2 captures prices with the same mean, but with a standard deviation between 5 and 7 times larger than Component 1. For the same four regions, Component 3 captures the price spikes, with substantial increases in average prices, a standard deviation between 126 and 209 times larger than Component 1, and infrequent occurrence between 0.4% and 2.7% of the time. For price formation in SA, Component 3 still has a higher average price, but variation that is only 4 times greater than that of prices in Component 1. Instead, Component 2 captures prices with extreme variation and a low incidence of 0.8%. The results suggest that price formation in SA differs from that in other regions. Certainly, the SA price distribution differs in Table 3, with more negative prices and higher average prices than other regions. While not reported here, interpretation of the mixture components are similar for the regressions of prices based on the other four supply regions.
There are five monotonically increasing estimates of each (inverse) supply function , which we combine into a single ensemble estimate , that is also monotonically increasing. Figure 4 plots the ensemble estimates of each supply function, and a number of observations can be made. The left hand side of each function is flat, as supply is made up of low cost baseline capacity, and a kink and rapidly increasing righthand side corresponds to capacity with rapidly increasing marginal costs. Table 1 provides a summary of registered generation capacity in each region at 2011. Coal is the lowest cost, followed by hydroelectric, gas and lastly ‘other capacity’ (the latter of which includes renewable sources). For example, in VIC there is 8.8MW of coal and hydroelectric capacity, which is approximately where the kink in the supply function occurs. In TAS, prices rise when supply goes beyond around 2MW, close to the hydroelectric capacity of 2.3MW.
Figure 5 plots the posterior mean estimates of the cost functions for the major interregional trade flows in Table 2. NSW is the major importer of electricity in the NEM, and panels (a), (b) and (c) show that imports in excess of 150 MWh on the Terranora, 1000 MWh on the QNI, and 1400 MWh on the NSWVIC interconnect, correspond to an increase in NSW prices. SA is the second major importer of electricity in the NEM, with panel (d) showing that increased flows from VIC correspond to higher SA prices. The Basslink interconnector is the only third party commercially operated interconnector in the NEM, with all others being owned and operated directly by AEMO. Imports into TAS produce only a limited impact on local prices, but exports to VIC in excess of around 520MWh from TAS correspond to a sharper increase in VIC prices. Nevertheless, increased flows on the interconnectors are unrelated to substantial regional price variation, which is consistent with a similar observation by Higgs et al. (2015).
6.2 Copula Models
The copula data are computed from the fitted marginals, and the five copula multivariate time series models estimated using MLE for each supply region . We use the Bayesian information criterion (BIC) to identify a lag structure from all thirtyfive combinations of lags between 1 and 5 halfhours, and also lags at the same time of the day between 1 to 7 days previously (i.e. lags at 48, 96, 144, … , 336 halfhours). An almost identical Markov structure was identified for all five models, with lags at the same halfhour of the day 1,2,3 and 7 days previously, and also at the three (TAS, VIC supply regions) and four (NSW, QLD, SA supply regions) halfhours immediately previous. In all cases, there is substantial serial and crosssectional dependence in the copula data.
To summarize the overall dependence structure, Figure 6 plots the autodependence matrices for lags , 1, 2, 3, 48 and 335, arising from the copula model with the VIC supply regressions as margins. These are computed from the autocovariance matrices , which are evaluated using sparse matrix algebra applied to the VAR(1) representation of a VAR() process; see Lütkepohl (2006; pp. 2631). Sparse calculations are important here as this involves algebra employing dimensional matrices with and .^{8}^{8}8This is implemented in Matlab in routines var2auto.m and plotautocorr.m found in the Supplementary Material. These measure dependence in prices, corrected for the impact of the VIC supply regression marginals, and a number of observations can be made. First, there is positive serial dependence in prices throughout. This is consistent with previous research that finds that there is significant time series dependence in electricity prices, even when corrected for a number of demand and supplyside variables (González et al. 2012). Second, this serial dependence declines as the lag increases, which is due to the assumption of stationarity in the copula time series model. The decline is fastest for the prices in QLD, with , and slowest for prices in VIC, with . Third, QLD prices are least dependent with those in other regions. We note also that QLD is the region with the lowest level of average imports and also lowest prices. Last, VIC, SA and TAS exhibit the strongest crosssectional dependence in prices at all lags, suggesting a higher level of market integration between these regions.
6.3 Event Studies
A key feature of the proposed econometric model is that it can be used to measure the response in regional prices to supply shocks or price impulses in any one region. Supplyside shocks are transmitted to prices through the supply function estimates, while price impulses are transmitted to future prices through the copula multivariate time series model.
6.3.1 Supply Shock
Consider the impact of a hypothetical supply shock in region at time by an amount , equivalent to a major baseline generator being taken offline. The marginal impact on prices can be evaluated using the th supply region regressions, but where the supply curve is shifted to the left by , producing a new supply curve . The expected (logarithm of) price in region from this model is then , where . ^{9}^{9}9Because we are modeling the logarithm of prices, the expectation of actual price can be computed as the sample mean of exponentiated iterates simulated from the marginal regression model in Equation (2.3).
To illustrate, we consider a supply shock of 560 MWh in VIC, which is equivalent to Steam Turbine Number 1 at the Loy Yang A Power Station going offline. The Loy Yang A Power Station is the largest single registered baseline generator in the VIC region, and the loss of Steam Turbine Number 1 is equivalent to the loss of approximately 4.6% of total 2011 registered capacity in that region. Outages of this size or larger occur from timetotime in the NEM due to plant or line failure. Figure 7 plots the posterior means and for each of the five price series . Using these estimates, along with the other components of the fitted model, we compute the marginal expected prices with and without the shock on 30 June 2010 at 08:00. Supply during this halfhour was MWh, which is the 99.3th percentile recorded in the data for this halfhour, and the 97.95th percentile over all halfhours, so that VIC supply is at a high level before the shock. Therefore, the shock results in a substantial increase in prices, with the vector of increases in marginal expected prices being (194.70, 1.28, 32.74, 1.85, 5.07) $/MWh for (NSW, QLD, SA, TAS, VIC). The increase is greatest in NSW, which Figure 7(a) shows has an estimated supply function with a kink at the lowest VIC supply value (around 8GWh). Prices in QLD, SA and VIC are less responsive to VIC supply shocks, with kinks occurring at higher values of VIC supply.
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.^{10}^{10}10For 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 halfhours 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 halfhours 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 halfhours 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 nonsmooth nature of the presented densities is not due to either Monte Carlo sample error,^{11}^{11}11We 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 nonsmooth.^{12}^{12}12Alternatively, 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 halfhours 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.^{13}^{13}13Hourly, rather than halfhourly, data are used to reduce the computational burden of the validation study. However, the models fitted to hourly and halfhourly 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.^{14}^{14}14In 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 crosssectionally 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 crosssectional 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 hourbyhour 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 demandweighted (log) price across all five regions,
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.^{15}^{15}15The 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 modelbased 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 highdimensional Gaussian copula to capture additional serial and crosssectional 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 meancorrected 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 twostage 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 interregional 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 5dimensional copula to capture the dependence between prices in the NEM, while Manner, Türk & Eichler (2016) employ a 4dimensional copula to capture the dependence between price spike incidence in four regions in the NEM. However, these lowdimensional copulas capture only crosssectional dependence, whereas we employ a parsimonious dimensional copula to capture both crosssectional 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 lognormal. 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 crosssectional and serial dependence in prices overandabove that induced by supplyside relationships in the marginal regressions. We illustrate the usefulness of our model in two event studies. The first is a supplyside 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 supplyside 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 interregional 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 timevarying variance, skewness, kurtosis and tail probabilities. The extent to which our proposed model captures interregional 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 reparameterize 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
where
is the tth row of , and
is the prior distribution defined in Section 4.1. Letting
and following Shively et al. (2011), we introduce a scalar latent variable such that
(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