A Regime-Switching Model for Electricity Spot Prices Gero Schindlmayr EnBW Trading GmbH [email protected] May 31, 2005

A Regime-Switching Model for Electricity Spot Prices

Abstract Electricity markets exhibit a number of typical features that are not found in most financial markets, such as price spikes and complex seasonality patterns. This paper proposes a stochastic model for electricity spot prices that is based on a regime-switching approach applied to average daily prices . Two different regimes represent a “normal” and a “spike” regime, the latter characterized by high volatility and strong mean-reversion. The model is calibrated via a maximum-likelihood optimization in connection with a Hamilton filter for the unobservable regime-switching process. Given the daily prices, the hourly price profiles are modelled using a principal component analysis for the 24-hour price vectors and afterwards setting up a time-series model for the factor loads. Example results are shown for spot price data from the European Energy Exchange EEX.

Keywords: electricity prices, regime-switching, PCA, Hamilton filter, seasonality

1

Introduction

In the course of the liberalization of energy markets it has become very important for utilities and energy trading companies to develop stochastic price models to be used for option pricing and risk management. A particular challenge is modelling electricity prices, since their stochastic behaviour differs substantially from other commodities or financial products. This section lists the most important properties of electricity spot prices and gives an overview of the model proposed.

1.1

Behavior of electricity spot prices

Electricity spot markets exhibit a number of typical features that are not found in most financial markets. The most important of those features are Seasonality: Electricity spot prices show a seasonal pattern on different time scales (yearly, weekly, daily) that reflect the typical electricity demand patterns. Spikes: Spot prices may exhibit extreme price spikes at times of high electricity demand (e. g. cold weather) and limited production capacity (e. g. plant outages). No cash-and-carry arbitrage: There is no cash-and-carry arbitrage relation between spot and futures prices, since electricity is not efficiently storable. Thus, the forward curve may have a very complicated seasonal pattern. Figure 1 shows the hourly spot prices of the European Energy Exchange (EEX) from year 2001 to May 2005.

1.2

Model overview

Most published work on electricity prices avoid the additional complexity of hourly spot prices by dealing with a daily average price process. Regime-switching models for daily electricity prices were studied in [4] for a continuous-time setup and in [3] and [1] for a discrete time setup. In [3] and [1] the processes for the stable regime and the spike regime were considered to be independent, which simplifies the analytical treatment of the model. in [6] a regime switch is triggered by crossing a certain price threshold. Modelling spot prices on an hourly basis, two approaches could be taken 1. Specify a scalar price process on an hourly granularity 2. Specify a 24-dimensional vector process on a daily granularity The first approach is taken in e. g. [2]. To account for the 24h-seasonality, seasonal ARMAprocesses with a 24h time lag were used. If instead of ARMA processes more complicated models such as regime-switching models are to be used, the 24h seasonality can be difficult to take into account. The second approach uses a vector process instead of the 24h seasonality. Following this approach, [7] apply a principal component analysis (PCA) to a load process, that is considered as a fundamental driver for the spot price. Here, the PCA is applied directly to the hourly spot price profiles. This PCA approach serves two purposes: First, we can choose 1

EEX spot hourly 1000

EUR/MWh

800

600

400

200

05 20

04 20

03 20

02 20

20

01

0

Date

Figure 1: EEX hourly spot price independent stochastic processes on a daily basis for the factor loads without the hourly granularity. Second, taking only the most significant principal components into account, the model can be reduced to less dimensions which speeds up computations and reduces memory usage. Regarding the spot process on a daily basis may also seem more natural with respect to actual trading, where spot prices are usually the result of an auction for the day (or weekend) ahead. The following notation is used throughout the paper: • St : Daily spot price on day t • St = (St1 , . . . , St24 ): Vector of hourly spot prices on day t • st = log St : vector of logarithmic hourly spot prices • st =

1 24 24 ∑i=1 st :

mean logarithmic price

We model the hourly spot prices on a logarithmic scale and separate the mean level and the hourly profile: 24

st = st + ht ,

∑ ht = 0

(1)

i=1

For the scalar process st and the vector process ht we use the following time series models: 1. The non-seasonal component of st is a AR(1) model with regime-switching 2

2. The stochastic component of ht is decomposed via a principal component analysis (PCA) into factor loads uti (i = 1, . . . , 24) which are then modelled as independent ARMA processes. Business days and non-business days are modelled as independent processes. This is an approximation done mainly for practical purposes. Non-business days have a completely different behavior regarding volatility and spikes. Thus, fitting a single process to both business and non-business days generally will give not very good results. Obviously, in reality prices on business days and non-business days are not independent but correlated. However, this effect is neglected here and could be incorporated in future.

2

The Daily Price Process

This section is devoted to the (logarithmic) daily price process st . This process has a yearly seasonality and is sensitive to weekday and holiday. To account for price spikes in the daily prices, a regime-switching approach is followed.

2.1

Seasonality

As a first step we take out the seasonal component of st . For this purpose we set up a regression model: Nd

st =

∑ 1Jd (t)βdA + 1JdDT (t) cos(2πt/365)βdB + 1JdDT (t) sin(2πt/365)βdC + 1JdDT (t)tβdD

d=1

+ 1J DT (t)1JV P (t)βdE + yt ,

(2)

d

where yt is the residual and 1J (t) the indicator function  1 t ∈J 1J (t) = . 0 otherwise The sets JdDT , d = 1, . . . , Nd define a partition into day types (e. g. Mo, Tu-Th, Fr, Sa, So, Holidays) and JV P defines vacation periods. The regression coefficients βdA , . . . , βdE have the following meaning coefficient βdA βdB , βdC βdD βdE

description mean level Amplitudes for yearly seasonality deterministic drift price effect of vacation period

The regression model is solved using a least squares algorithm. An example of the seasonal component is given in figure 2.

3

5 log price daily seasonal component

4

log price

3

2

1

05 20

04 20

03 20

02 20

20

01

0

Date

Figure 2: Seasonality of the daily logarithmic price from EEX

2.2

Stochastic model

As a next step, a regime-switching ARMA process is calibrated to the yt data. At this stage, the process yt is divided into business days and non-business days: yt = 1J B (t)ytB + 1J H (t)ytH , where J B denote the business and J H the non-business days. In the following we describe the time series model and calibration for either of the processes ytB or ytH . To simplify the notation we leave out the superscript and work with a time series yk , k = 1, . . . , N observed at times tk . To account for price spikes, we choose an AR(1) model with regime switching, where the regime change is modelled via a discrete Markov chain. Thus, the model equation is given by
yk − µrk = φrk yk−1 − µrk−1 + σrk εk , (3) where rk ∈ {1, 2} denotes the current regime at time tk and εk ∼ N(0, 1). The Markov chain is characterized by the transition probability matrix   p11 p21 P= p12 p22 4

with pi j = P(ri = j | rk−1 = i). Since the regime state is not observable we use the Hamilton filter described in [5]. Since the Hamilton filter requires the transition probabilities to depend only on the current state, we have to transform the two-state model above into a four-state model (see [5], p. 691) with states r˜k ∈ {1, 2, 3, 4}, such that r˜k = 1 r˜k = 2 r˜k = 3 r˜k = 4

if rk = 1 and rk−1 = 1 if rk = 2 and rk−1 = 1 if rk = 1 and rk−1 = 2 if rk = 2 and rk−1 = 2.

The transition probability matrix becomes  p11 0 p11 0  p12 0 p12 0 P˜ =   0 p21 0 p21 0 p22 0 p22

   

In this way, the model is of the general form yk = µ˜ r˜k + φ˜r˜k yk−1 + σ˜ r˜k εk .

(4)

i for the probability P(r = i | y , y The Hamilton filter produces estimates ξˆk|k k k k−1 , . . .) that the observation was generated by regime i based on observations up to time k. As a byproduct of the Hamilton filter, the log maximum likelihood function is calculated as ! Nk 4 L= log ξˆ i η i (5)

∑

∑

k|k k

i=1

k=1

Example calibration results for EEX spot price data from Jan 01 to May 05 are the following:

regime 1 regime 2

µ −0.004 −0.02

φ σ 0.74 0.11 0.68 0.30

 P=

0.94 0.28 0.06 0.72



Since we have only estimated probabilities about the regime during the calibration period, we can only calculate estimates about the residuals: yk − µ˜i − φ˜i yk−1 εˆk = E[εk ] = ∑ P(˜rk = i | yk , yk−1 , . . .) σ˜ i i=1 4





A sample autocorrelation function of νˆ k is shown in figure 3. Figure 4 shows a QQ plot for the innovations both for the regime-switching process with two regimes and a non regimeswitching AR(1) process. Using a simple AR(1) process instead of the regime-switching process, the innovations are clearly non-Gaussian. The innovations εˆk of the regime-switching process provide a much better fit to a Gaussian distribution. 5

Sample Autocorrelation Function (ACF)

0.8

Sample Autocorrelation

0.6

0.4

0.2

0

−0.2

0

2

4

6

8

10 Lag

12

14

16

18

20

Figure 3: Autocorrelation of EEX daily price residuals (Jan 01 – Apr 05) QQ Plot of Sample Data versus Standard Normal 8 no regime−switching regime−switching 6

Quantiles of Input Sample

4

2

0

−2

−4

−6 −4

−3

−2

−1

0

1

2

3

4

Standard Normal Quantiles

Figure 4: QQ Plot of innovations for daily price process (EEX Data Jan 01 – Apr 05)

3

The Hourly Profile Process

From equation 1 the profile process is defined as the logarithmic spot process normalized to mean zero: ht = st − st The steps for modelling the profile dynamics are the following: 1. de-seasonalize the profile for each hour 2. use PCA to decompose the vector process into factor loads 3. model each factor load as an ARMA process 6

3.1

Seasonality and PCA decomposition

Daily spot price profiles have a pronounced yearly and weekly seasonality (see figure 5) due to different temperature and light conditions. To de-seasonalize the profiles a similar technique as described in section 2.1 is used for each single hour. In this way the profiles can be written as ht = hˆ t + ∆ht , where hˆ t is the deterministic seasonal component and ∆ht is the stochastic residual from the regression.

summer winter 2

Value

1.5

1

0.5

0

5

10

15

20

Hour

Figure 5: Mean summer and winter hourly profiles The residuals ∆ht form a vector process with a covariance structure reflecting the different volatilities of each single hour and the correlations between the hourly profile values. As before, the process ∆ht is divided into business days and non-business days: ∆ht = 1J B (t)∆htB + 1J H (t)∆htH . The two processes ∆htB and ∆htH are treated separately. For notational convenience we regard only one of these processes and denote it by ∆hk , k = 1, . . . , N. By means of a PCA the vector process ∆hk is now decomposed into independent factor loads. Let p1 , . . . , p24 be the principal components of ∆hk (after column-wise normalization to standard deviation 1), i. e. the eigenvectors of the matrix (∆h)T ∆h/N, where ∆h is a N × 24-matrix. It is assumed that the corresponding eigenvalues are in descending order λ1 > . . . > λ2 4. The 24×24-matrix (p1 , . . . , p24 ) containing the principal components as columns is denoted by P. The N × 24matrix of factor loads is then given by Z = (∆h) P. An example of the principal components is shown in figure 6.

7

The i-th eigenvalue λi represent the variance contributed by the i-th principal component. Figure 7 shows the percentage of the total variance that is explained by a given number of principal components. In the case shown there, a large number of principal components is needed explain the variance.

0.4

0.3

0.2

Value

0.1

0

-0.1

-0.2

-0.3

pc 1 pc 2 pc 3

-0.4

-0.5

0

5

10

15

20

Hour

Figure 6: First three principal components (EEX spot data, business days)

100

explained variance in %

90

80

70

60

50

40

business days non business days 0

5

10

15

20

number of pc

Figure 7: Explained variance with respect to the number of principal components (EEX data)

8

3.2

Stochastic model

Let zik = Zk,i , i = i, . . . , 24, k = 1, . . . , N be the time series of the i-th factor load, which is the i-th column of the matrix Z. The series zik are modelled as independent ARMA process. From the simulated factorloads Z the process ∆ht can be calculated as ∆ht = ZPT . Figure 8 shows the autocorrelation function of the innovations calculated for the first factor load with a AR(1)-model and an ARMA(1,1)-model respectively. The results are qualitatively similar for the factor loads of the next higher principal components. ARMA(1,1): Sample Autocorrelation Function (ACF)

0.8

0.8

0.6

0.6 Sample Autocorrelation

Sample Autocorrelation

AR(1): Sample Autocorrelation Function (ACF)

0.4

0.4

0.2

0.2

0

0

−0.2

0

2

4

6

8

10 Lag

12

14

16

18

−0.2

20

0

2

4

6

8

10 Lag

12

14

16

18

20

Figure 8: Autocorrelation plot for ARMA-innovations of first factor load A more careful study of the profile dynamics reveals that on days with an extremely high daily price the profiles behave differently, since those high daily prices often are the cause of spikes at single hours. Thus, to get a more realistic price behavior, on days where the regime is in state 2 (volatile regime) the model randomly draws a historical profile according to the given day type, yearly season and regime state.

4

Simulation Results

While each of the previous sections focused on a separate part of the model, this section shows simulation results of the complete model. In figure 9 ten sample simulation paths for one year are shown. This demonstrates that the model indeed produces a spike behavior very similar to the behavior observed in historical EEX spot price data. A more explicit analysis is given in figure 10 where the histograms of the historical and simulated spot prices over the time period Jan. 2001 to May 2005 are compared for the price range 100 to 250 EUR/MWh. Again, the qualitative behavior is very similar in this range of high spot prices.

9

path path path path path path path path path path

1600

1400

Hourly Price EUR/MWh

1200

1 2 3 4 5 6 7 8 9 10

1000

800

600

400

200

07 n Ja

ec

06

06 D

ov

06 N

O ct

Se

p

06

06 g

l0

6 Au

Ju

06 n Ju

ay

06

06 M

r

06 Ap

ar

b

06 M

Fe

Ja

n

06

0

Date

Figure 9: Ten sample hourly EEX spot simulation paths for year 2006

−3

x 10

historical data simulated data

3.5

3

frequency

2.5

2

1.5

1

0.5

0 100

150

200

250

spot price in EUR/MWh

Figure 10: Histogram of historical and simulated spot prices in the range 100 − 250 EUR/MWh

10

References [1] M. Bierbrauer, S. Trück, and R. Weron Modelling electricity prices with regime switching models Lecture Notes in Computer Sciences, 07/2004 [2] M. Burger, B. Klar, A. Müller, G. Schindlmayr A spot market model for pricing derivatives in electricity markets. Quantitative Finance, Vol. 4, pp. 109-122, 2004 [3] C. De Jong, R. Huisman, Option formulas for mean-reverting power prices with spikes, Working Paper, Erasmus University Rotterdam, 2002. [4] S. Deng, Stochastic Models of Energy Commodity Prices and Their Applications: Meanreversion with Jumps and Spikes, Working paper, Georgia Institute of Technology, 2000. [5] J. D. Hamilton. Time Series Analysis. Princeton University Press, 1994 [6] A. Roncoroni, H. Geman, A class of marked point processes for modelling electricity prices, Working Paper, ESSEC, 2003. [7] P. Skantze, A. Gubina, M. Ilic, Bid-based stochastic model for electricity prices: The impact of fundamental drivers on market dynamics Report, MIT Energy Laboratory (2000).

11