Filtering and parameter estimation for electricity markets by Alberto Molina-Escobar B.Sc., Universidad Nacional Autónoma de Mexico, 1996 M. S., Universidad Nacional Autónoma de Mexico, 2000 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Mathematics) The University Of British Columbia October, 2009 © Alberto Molina-Escobar 2009 Abstract The growing complexity of energy markets requires the introduction of in creasingly sophisticated tools for the analysis of market structures and for the modeling of the dynamics of spot market and forward prices. In order for market participants to use these markets in an efficient way, it is important to employ good mathematical models of these markets. This has proved to be particularly difficult for electricity, where markets are complex, and ex hibit a number of unique features, mainly due to the problems involved in storing electricity. In this thesis we propose three models for electricity prices. All are multi- factor models, that is, as well as an observable spot price they assume the existence of an unobservable long term mean’ process. The introduction of such additional processes helps to explain the relation between spot and futures prices. In the first part of the thesis we introduce a two factor Gaus sian model for prices. Using the Kalman filter, and based on both spot and forward prices, we successfully estimate parameters for simulated data. We then estimate parameters for the German EEX market, and compare our fitted model with the observed prices. We find that this model does capture some features of the EEX market, but it fails to exhibit the price spikes which are a prominent feature of true spot prices. We therefore introduce a second model, which includes jumps. The inclusion of jumps has the potential to give a better explanation of the behavior of electricity prices, but it creates difficulties in the estimation of parameters. This is because as the model noise is non-Gaussian the Kalman filter cannot be applied satisfactorily. We implement the particle filter adopting the Liu & West approach for the jump model. This method allows us to identify the hidden process in the model, and to estimate a small number of parameters. The third model is a new model for electricity prices based on the inverse Box-Cox transformation. This model is non-linear with Gaussian noise, and can generate price spikes using fewer parameters than a multi-factor jump-diffusion model. In this context, we successfully applied the Unscented Kalman filter to estimate the parameters. 11 Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements viii Dedication ix 1 Introduction 1 1.1 Commodity markets 1 1.2 Electricity markets 2 1.3 The relationship between spot and futures prices 6 1.4 Previous work 8 2 Filtering 25 2.1 State space formulation . 25 2.2 The Kalman filter 27 2.3 The unscented Kalman filter 29 2.4 Particle filter 32 2.5 Parameter estimation via maximum likelihood 38 2.6 Parameter estimation via Bayesian methods 43 3 MROU model 44 3.1 Double mean-reversion model 44 3.2 Radon-Nikodym theorem for Ornstein-Uhlenbeck processes . . 45 3.3 Future price 49 3.4 Formulation in Kalman filter terms 51 in Table of Contents 3.5 Empirical results 56 3.5.1 Simulated data 56 3.5.2 The German electricity market . 58 4 MROU with jumps 66 4.1 Description of the model 66 4.2 Valuation of electricity futures 67 4.3 Particle filter setup 70 4.4 Simulated data with known parameters 76 4.5 Likelihood function estimation 77 4.6 Sequential parameters 78 4.7 Empirical results 80 5 NLMROU model 89 5.1 The model 89 5.2 Future price 90 5.3 Unscented Kalman filter setup and estimation procedure . 93 5.4 Simulation results 95 5.5 Parameter estimation based on historical data 95 6 Conclusions 101 6.1 Future work 102 Bibliography 105 Appendix A 113 Appendix B 116 iv List of Tables 1.1 Models for electricity prices 22 3.1 The data are taken from Hikspoors Jaimungal (columns 1 and 4), and Nomikos & Soldatos (columns 2 and 3) 47 3.2 Five different maximization runs, on the same set of simulated data 57 3.3 Estimation using one futures contract (average of 50 simula tions) 57 3.4 Estimation using two futures contracts (nrr300) 58 3.5 Estimated values for f(t) by least-squares fitting 60 3.6 Estimated values for the EEX market using St and F(t, T1) 63 3.7 Estimated values for the EEX market using S,, F(t, T1) and F(t,T2) 63 3.8 The table shows the first four moments of the logarithmic deseasonalized price returns of observed data and the average of 50 simulated trajectories 65 4.1 Sample of 8 estimated values for )‘x 84 4.2 Individual estimates for parameters in MROU with jumps model. 84 5.1 Five different maximization runs with 800 observations 96 5.2 Estimation using one futures contract (n = 800) 97 5.3 Estimated values for the EEX market using S, and F(t, T1) . 98 5.4 The table shows the first four moments of the logarithmic deseasonalized price returns of observed data and the average of 50 simulated trajectories 100 V List of Figures 1.1 Classification of commodity markets 1 1.2 World net electricity consumption 2004-2030 2 1.3 If the total load is low, the plants with the lowest variable production costs are used (nuclear, hydro); if the total load is high, gas or oil fired plants with high fuel cost are running additionally, producing a huge effect on the price 3 1,4 Seasonal patterns by hours and by week for the German market. 4 1.5 Average daily spot price in German market for years 2002-2007 4 1.6 The factors exerting a major influence on electricity wholesale price 5 2.1 A graphical representation of the particle filter with impor tance sampling and resampling 37 3.1 Sampling distributions with (column 4) where we add 10% deviation from Hikspoors & Jaimungal’s parameters (column 1) 48 3.2 The upper graph shows the simulated spot price St and the future price F(t, T1) with maturity of one month. The lower graph is the long-term mean process L 55 3.3 Electricity spot and nearby monthly futures price in German market 60 3.4 The upper graph shows the log spot-price of the EEX mar ket and the seasonal component h(t) and the lower graph the deseasonal series Xt = logS — h(t) 61 3.5 The upper shows the spot price St with exp{h(t)} and and the lower graph the deseasonal spot price S, = exp{Xt} 62 3.6 Simulation of spot and future prices (upper graph) and long term mean process (lower graph) using estimate values for part 1 64 vi List of Figures 4.1 Plot of the true state L and estimate of the particle filter. 76 4.2 The log-likelihood for different Ax values 77 4.3 The graph shows a simulated trajectory of the MROU model with jumps 82 4.4 The graph shows the upward and downward jumps for the simulated trajectory of the MROU model with jumps 83 4.5 Sample particle filter trajectories for the estimate of Ax. . 84 4.6 Sample particle filter trajectories for the estimate of T1.u and ‘id using Liu and West approach 87 4.7 Sample particle filter trajectories for the estimate of cxx and AL using Liu and West approach 88 5.1 The graph shows a simulated trajectory of the NLMROU model. 96 5.2 Plot of the true and estimated processes L and X of the NLMROU model (first 80 observations) 97 5.3 Simulation of spot and future prices price using estimated val ues for the whole data set 100 vii Acknowledgements First and foremost I must acknowledge my thesis supervisor, Martin Barlow. Throughout my time as a graduate student Martin has provided me with unending guidance in research, constant encouragement and many years of financial support. Without his undying support this thesis would not have been possible. I am also grateful for the helpful comments of Rachel Kuske, Ulrich I-Iorst, Joel Friedman and Kevin Murphy. Many thanks also to Arnaud Doucet for his feedback and recommendations that help me improve my project. I would also like to thank Gabriel Mititica for his help and encouragement. Thank you also to CONACyT for its financial support during the first years of my program. viii To my wife and daughter. ix Chapter 1 Introduction In this thesis we are interested in the commodity futures markets, and in particular in electricity futures markets. 1.1 Commodity markets In this section we are going to describe some of the unique characteristics present only in commodities markets. Figure 1.1 shows the three fundamen tals groups of commodities. Commodities Agricultural Metals Energy I Vegetable Animal Industrial Precious Upstream DownstreamGoods Goods • Corn • Live Cattle • Caper • Gold • Crude Oil • Coffe • Pork Bellies • Aluminium • Silver • Heating Oil • Electricity • Cotton • Lean Hogs • Lead • Platinum • Natural Gas Figure 1.1: Classification of commodity markets. Commodities markets exhibit some characteristics that are not present in financial markets due to physical constraints and also due to the variation of demand due to changes in consumption. The commodity spot price is defined by the intersection of supply and demand curves. Thus the spot price can be affected by changes in consumption, production or inventory. Unlike financial assets, which are traded for investment purposes, commodities are traded in order to be consumed or used in an industrial process, with the partial 1 Chapter 1. Introduction 35.000 30000 20,000 19,554 :::: 10,000 3004 2010 Figure 1.2: World net electricity consumption 2004-2030. exception of some precious metals. This close link with the real economy causes commodities prices to have seasonal behavior and also mean-reversion [31]. This is one reason why many of the standard financial theories may not be applicable to commodities markets. 1.2 Electricity markets Among the various commodities the energy market is the most recent market to be transformed. Since the early 1990s, electricity markets have been and continue to be developed as a result of the deregulation of electricity markets worldwide. In many regions the market structure has moved from a monop olistic to a competitive one. Traditionally, there was only one company or government agency that produced, moved, distributed, and sold electricity power and services. This transformation has been already taken place in the Americas (parts of Canada and US, Argentina, Chile, Peru, Paraguay and Colombia) in Europe (Norway, Finland, Denmark, Germany, France, Netherlands, Spain, Poland, and Romania) and in the Asia/Pacific region (Australia, New Zealand and Japan). In theory, deregulating the electricity market should increase the effi ciency of the industry by producing electricity at lower costs and passing those cost savings on to customers [39]. Electricity is a growing market. In 1973 electricity consumption accounted for 11% of the total world energy Year 2 Chapter 1. Introduction 0 Figure 1.3: If the total load is low, the plants with the lowest variable pro duction costs are used (nuclear, hydro); if the total load is high, gas or oil fired plants with high fuel cost are running additionally, producing a huge effect on the price. demand and has grown to 18% today. The absolute growth rate of electricity consumption in the future is estimated at an average of 2.4% per year. The projected growth in electricity consumption is shown in Figure 1.2. Electricity is considered a secondary energy source, which means it is created from the conversion of other sources of energy, such as coal, oil, natural gas, nuclear power, or hydropower, all of which are referred to as primary energy sources. To understand the behavior of electricity prices we have to note that electricity possesses a unique feature; it is very difficult and expensive to store and quite difficult to transmit from one region to another. As a result of this, the spot price of electricity is set by the short-term supply-demand equilibrium, and supply and demand must be in balance at each time. Figure 1.3 displays a schematic supply-demand curve. The sup ply and demand are affected by many factors that influence the seasonality and volatility of prices. For example, supply may be affected for transition constraints (breakdowns) or fluctuation of fuel prices (oil, gas). Demand exhibits seasonal fluctuations, which are due to climate conditions. In ad 1 (ETA) System for the Analysis of Global Energy Markets (2007). Demand Capacity 3 0000 Chapter 1. Introduction 0111111 M T s (a) Hours average price from Jan-Dec 2002. (b) Daily average log-price throughout the week from 12/2002 to 05/2005. Figure 1.4: Seasonal patterns by hours and by week for the German market. 350 300 250 200 150 100 50 07/01/02 01/01/03 07/01/03 01/01/00 07/01/04 01/01/05 07/01/05 00/01/06 07/01/06 01/01/07 Figure 1.5: Average daily spot price in German market for years 2002-2007 4 Chapter 1. Introduction dition, electricity demand is also not uniform through the week. It peaks during weekday working hours and is low during nights, holidays and week ends due to lower industrial activity, see Figure 1.4. Also unexpected weather conditions can cause abrupt and dramatic disruptions, producing jumps and spikes in the spot price. Finally the constraints on transmission mean that power markets are geographically distinct. In some markets (such as Alberta or Norway), demand is higher in the winter months due to the use of power for heating. In other markets, such as California power, usage peaks in the summer due to use of electricity for air conditioners. Figure 1.6 shows the factors that influence the determination of electricity prices.2 Long-term factors • Economic cycle • Politcal decisions • Capacity expansionfcloaures ThameS power plants plonta Electricity Power prices I Lighting end consumer behavior 4)i Vacations PnblIe hotidaysCloudiness Time a? day l Factors of a.ipply Factors of demand Factors of supply and demand Figure 1.6: The factors exerting a major influence on electricity wholesale price. The most unusual feature of electricity spot prices is the presence of “price spikes”; a phenomenon which does not have any parallel in other commodity markets. See for example Figure 1.5 which shows that in some days in July 2Source: (RWE AG) Shares of Primary Energy Sources in Total Electricity Generation in Europe (2008). 5 Chapter 1. Introduction 2006 the spot price in Germany reached 300 €/MWh, compared with a normal daily price of 30-50 €/MWh. If such an event occurred for a conventional commodity, such as say cop per, holders of the material would be able to made substantial profits by selling the commodity during the spike, and then repurchasing it at a nor mal price a few days later. But, because, it cannot easily be stored, this is not possible with electricity. 1.3 The relationship between spot and futures prices The relationship between the spot price and the futures (forward) prices is important for risk management and option pricing theory [30]. Across all these commodities ranging from agricultural products to pure financial assets certain common principles of futures valuation and futures price behavior apply. Let (, F, IP) be a complete probability space endowed with the natural filtration {}>o. In the financial world, the relation between spot and future prices, under the risk-neutral measure Q, is given by the formula Q/ I T ‘1F(t,T)=E ISTexpj’ rdnj} (1.1) where r is the risk-free interest rate. The proof is based on the no-arbitrage argument (see [9]), which proceeds by comparing returns on a portfolio con sisting of the future contract with one consisting of cash and the commodity. However, unlike financial assets, storage of commodities is costly. Con sequently, physical ownership of the commodity carries an associated flow of services. On the one hand, the owner enjoys the benefit of direct access, which is important if the commodity is to be consumed. On the other hand, postponing consumption and storing the commodity means that storage ex penses have to be paid. The net flow of these services per unit of time is called the convenience yield Ct. Since the convenience yield is the result of subtracting the cost of storage from the benefit attached to the physi cal commodity it can be both positive or negative at different times. (A positive convenience yield implies an instantaneous benefit from holding the commodity, a negative one an instantaneous cost.) 6 Chapter 1. Introduction Again, by a no-arbitrage argument, the relationship between the spot and forward price is given by: F(t, T) = (ST exp { ftT(ru — cu)du}), (1.2) where c is the instantaneous forward convenience yield [31]. Note that the convenience yield plays the same role as dividends play for stocks. Some authors have argued that as a consequence of the non-storability of electricity the notion of convenience yield is irrelevant in power markets. Therefore the relation between spot and futures (forward) prices cannot be established through the no-arbitrage argument (see [14, 32, 33]). For exam ple, Geman and Roncoroni comment in [32]: “Our view is that a convenience yield does not really make sense in the context of electricity: since there is no available technique to store power (outside of hydro), there cannot be a benefit from holding the commodity, nor a storage cost. Hence, the spot price process should contain by itself most of the fundamental properties of power.” Other financial theories view the futures (forward) prices F(t, T) and the expected future spot price Et(ST) as related but not identical. The difference is the risk premium, i.e. F(t, T) Et(ST) + ir(t, T). (1.3) The full specification is not straightforward to establish. The theory of a positive risk premium is termed normal backwardation. The opposite situa tion where the futures prices is set above the expected future spot price (a negative risk premium) is called contango. An alternative approach is the actuarial one, which values a forward con tract as its discounted expected real world payoff, see [44]. This is the ap proach we will adopt in this thesis: we will assume that the relation between spot and futures prices is given by (1.2), and that the risk free interest rate r,, and convenience yield c are constant. 7 Chapter 1. Introduction 1.4 Previous work The main motivation for the development of models for electricity prices is the need for such models by market participants. For example, a power company has the choice of selling its power either on the spot or forward market, and would wish to make the optimum choice. In addition there is the need to price derivatives such as forwards, options and swaps. Hence the model should be sufficiently sophisticated for realistic modelling but sufficiently simple for pricing of derivatives. This issue is very important for computing risk measures, testing hedging strategies and evaluating investment policies. Various approaches have been developed to describe the stochastic price process in energy markets. There are significant parallels between corrimodity markets and interest rate markets. For commodity markets, the traded assets are both the spot and various forward or future contracts. For interest rates, the main traded assets are futures (represented by different types of bonds), while the spot or instantaneous rate of interest plays a more minor role. Given these parallels, it is natural to use the interest rate theory as a base for electricity price models. In general, interest rate models can be separated into two categories: short-rate model and forward-rate models. The short-rate models describe the evolution of the instantaneous interest rate as stochastic process, and the forward-rate models capture the dynamics of the whole forward curve (Heath-Jarrow-Morton models). These interest rate models are then applied to arrive to arbitrage-free pricing of bonds or other derivative products. The same division of models arises for power prices, where the models may be broadly divided into two groups:3 statistical models (spot price based models) and fundamental models (forward based models). For the forward based models, the futures prices are the main objects of study, and the dynamics of the whole futures prices curve is modeled using the Heath-Jarrow-Morton [42] theory for interest rates. See for exam ple, Clewlow and Strickland (1999) [19] and Manoliu and Tompaidis (2002) [61], and for more recent papers see Borovkova (2006) [10], Koekebakker and Ollmar (2005) [55]. 3There is another approach based on econometric time series model that we will not consider in this work (see [56, 63]) 8 Chapter 1. Introduction A general discussion of HJM-type models in the context of power futures is given in Benth and Koekebakker (2008) [8]. They dedicate a large part of their analysis to the relation of spot, forward and swap-price dynamics and derive no-arbitrage conditions in power future markets and conduct a statis— tical study comparing a one-factor model with several volatility specifications using data from the Nord Pool market. The disadvantage of such approaches is that futures prices do not reveal information about price behavior on a daily timescale and provide a poor approximation to the complex observed spot behavior in power markets. In this thesis, following most of the literature, and the philosophy outlined by Geman and Roncoroni above, we will consider spot based models. In principle these models should provide a reliable description of the evolution of electricity prices. Moreover, these models are versatile in the sense that it is relatively simple to aggregate characteristics to an existing family or class of models by for example adding a seasonality function. Securities (stocks) are usually modeled by Geometric Brownian Motion with drift S, = Soexp{at+aW}, as in the famous Black-Scholes model. This model is not found suitable for commodities, since ‘mean reversion’ is typical feature of these markets [30, 31, 66]. The simplest stochastic process with mean-reverting behavior is the Ornstein Uhlenbeck process [661. Here the process X, is a diffusion process satisfying the stochastic differential equation dX = —A(X — a)dt + crdWt (1.4) where W is a standard Brownian motion, a the volatility of the process, and A the velocity with which the process reverts to its long term mean a. Many electricity price models use this process or variants as a basic building block. For example, Lucia and Schwartz (2002) [59], give models of the form St = h(t) + X (1.5) or S = exp {h(t) + X} (1.6) 9 Chapter 1. Introduction where S is the spot price, X, is an Ornstein-Uhlenbeck process, and h(t) is a deterministic component, intended to account for seasonal and weekly effects. Benth et al. (2008) [6] called models like (1.5) ‘arithmetic models’ and (1.6) ‘geometric models’, i.e. geometric models represent the logarithmic prices by a sum of processes. The incorporation of a deterministic component of this kind is an im portant feature of nearly all spot price models. Spot prices are higher on weekdays than on weekends, due to higher demand, so a correction h(t) which compensates for this is essential. See for example Figure 1.4b. Spot price models can be divided into ‘single’ factor or multi-factor mod els. For single factor models the spot price is itself a Markov process, while in multi-factor models the spot price is a function S = g(X’, ..., X) of a. mul tidimensional Markov process. Here g : IRk R+, and as g is not one-to-one these models have unknown or hidden components. As well as the model of Lucia and Schwartz mentioned above, other single factor models are in Cartea and Figueroa (2005) [16], Barlow (2002) [3], Kanamura and Ohashi (2007) [49], and Geman and Roncoroni (2006) [32]. Many of these models, unlike that of Lucia and Schwartz, include mecha nisms to take account of price spikes. One of the simplest of these is in Cartea and Figueroa [16], which adds a jump term to the Ornstein-Uhlenbeck pro cess: log St = h(t) + Y, d = —cYdt + crdWt + JdN, (1.7) where W, is a Brownian motion, h(t) is assumed to capture the seasonal patterns of the spot price, and the third term JdN enables the process to have discrete random spikes: these are a combination of a Poisson process, which determines the jump frequency, and a jump-size distribution, which gives the jump magnitude conditional on a jump occurring. In (1.7) the process dN is approximated by a Bernoulli process with parameter ldt and J is log-Normal, i.e. log J N(—u2/2,u2). Cartea and Figueroa apply this one-factor mean-reverting jump diffusion model for the electricity spot price, adjusted to incorporate seasonality effects and derive the corresponding forward in closed-form to the England and Wales market. 10 Chapter 1. Introduction However, the rather short period for which electricity prices were available and the small number of spikes caused difficulties with parameter estimation. Such models require a high speed of mean reversion in order to reduce the spot price following a large positive jump, and this has the effect of removing too much variability in the series over the non-jump time-periods. Barlow [3] introduces a nonlinear Ornstein-Uhlenbeck model for spot power prices. The price is obtained by matching the demand level with a deterministic supply function which must be nonlinear to account for price spikes. He proposes the inverse function of the Box-Cox transformation. I fa(Kt)’ 1+tX > E0St = , 1 + dX —..\(Xt—a)dt+udWt, where f(x) (1 + x)’, 0 fo(x) e. When o = 0, an exponential Ornstein-Uhlenbeck process is retrieved for St. The case — 1 yields a regular Ornstein-Uhlenbeck process. The model has been estimated by maximum likelihood on the Alberta and California markets. Another paper sharing the same theoretical idea is found in Kanamura and Ohashi [49]. Instead of using the inverse function of the Box-Cox trans formation they assume that the supply curve has a ‘hockey stick’ shape. Setting X = D — D, D describes the seasonal component and s — f (a1 +b1D), D D0(a2 +b2D), D > D0 dX = (—\X)dt+udW. This model captures the observed mean-reverting behavior of electricity mar kets and it accounts very well for the observed price spikes, allowing for a better fit to market data. But the assumption of a deterministic supply 11 Chapter 1. Introduction function is probably too restrictive since this implies that spikes can only be produced by surges in demand. Geman and Roncoroni [32] built up a jump-reversion model for electricity spot prices. The model assumes that the natural logarithm of power price dynamics is described by a stochastic differential equation dE(t) = [h(t) + 8(t(t) — E(t))]dt + crdW(t) + f(E(t))dJ(t), (1.8) where h(t) is a deterministic seasonality function, 8 is the mean reversion speed, and o is a constant instantaneous volatility, i(t) is the mean rever sion level. The process reverts to a deterministic mean level rather than the stochastic pre-spike value. The last term in equation (1.8) represents the discontinuous part of the model featuring price spikes. This effect is charac terized by three quantities defining occurrence, direction, and size of jumps. f is a function which is ±1 depending on the level of the spot price. E1t — f +1, if E(t) <r(t) ‘ “ “ — —1, if E(t) r(t) The process J(t) is a time-inhomogeneous compound Poisson process with intensity function A(t) = (1 + I sin[n(t - 7)/61I _i) where the expected maximum number of jumps per year is represented by Ic. Jump sizes are modeled by a sequence of independent and identically distributed truncated exponential variables. This model generates trajectories similar to those observed in the elec tricity market, and also it gives a good fit of the empirical moments of order 1, 2 and 4, i.e. mean, variance and kurtosis.4 Neither of the last two models includes the convenience yield as a factor, nor considers the valuation of futures contracts or any other kind of deriva tive. The single factor models are quite tractable and their parameters are 4The kurtosis of a random variable X with mean m and variance a2 is defined by: E((X-m)) When i is much greater than 3, it means that the density in the tail is higher than that which prevails for a Gaussian distribution. 12 Chapter 1. Introduction relatively easy to estimate. However they have a serious limitation: they do a poor job explaining the relation between spot and futures prices, see [3] and [12]. This limitation can be avoided if changes in spot prices are allowed to depend on more than one factor. The copper mine example of Brennan and Schwartz (1985) [11] assumed that the spot price followed a geometric Brownian motion arid incorporated a convenience yield to their model, assuming it was proportional to the spot price. dS = ,LtSdt+uSdz, C(S,t) = cS. The idea of a constant convenience yield only holds under restrictive assump tions, since the theory of storage is rooted in an inverse relationship between the convenience yield and the level of inventories. Gibson and Schwartz (1990) [35] take an important step to a more realistic model of the econ omy by introducing a stochastic convenience yield rate. The spot price S of the commodity is described by a geometrical Brownian motion and the convenience yield rate cY is described by an Ornstein-Uhlenbeck process with equilibrium level a and rate of mean-reversion f: dS = (—öt)Stdt+ciiStdzi, dc5t = ic(cr — S)dt +u2dz, dz12 = pdt. Significant contributions have been made by Schwartz (1997) [73]. He reviewed one and two factor models and developed a three factor model under stochastic convenience yield and interest rates. Including the interest 13 Chapter 1. Introduction as a third factor makes forward and futures prices different. dS = (rt — 6)Sdt + uiSdzi, d5 = ii(a’—c5t)dt+u2dz, dr = a(m — rt)dt + cr3dz, dz12 pidt, dzdz3 = p2dt, dz1 = p3dt. This model was originally developed for two commercial commodities (copper and oil). He used the Kalman filter algorithm to estimate the parameters in the models. In the paper by Lucia and Schwartz [59], where they analyze the Nordic Power market, the spot price is modeled by = dX = —.AXdt+uxdWx, dY = ,udt+ciydWy, dWxdWy = pdt. The function h(t) is deterministic, and it is intended to capture the pre dictable component in the spot price, i.e. seasonal effects. This function distinguishes between weekdays, and includes a monthly seasonal compo nent employing dummy variables. The idea of this model is to have a non- stationary process for the long-term equilibrium price level Y and short-term mean-reverting component X. They estimated all the parameters simulta neously by nonlinear least squares methods. The multi-factor models described so far do not capture one of the most characteristic feature of power prices, jumps or spikes. Several authors, Deng (2000) [24], Villaplana (2003) [80], and Xiong (2004) [82] extend such models to even more factors with both diffusion and jumps. In the work of Villaplana power prices are modeled according to non-observable state variables that account for the short-term movements and long-term trends in electricity 14 Chapter 1. Introduction prices. inSt = h(t)+X+Y dX ‘cxXtdt + uxdWi + JdN()) — JddN(\d) dY = —ky(i—)dt+uydW2 dW12 = pdt. The jump components are characterized by N(\), and N(Ad), i.e. Poisson processes with intensities )4, and ‘\d respectively, and by random jumps of size J,, and Jd with some specified distribution (Gaussian/Exponential). Deng (2000) and Villaplana (2003) set their models in the affine jump diffusion (AJD) framework which enables them to use transform results of Duffie et al. (2000) [27] to derive tractable closed-form solutions for a variety of contracts. Deng proposes more sophisticated mean-reverting jump dif fusion models with deterministic/stochastic volatility and regime switching, which may be a good way of addressing the dramatic changes in spot prices. However the trajectories produced by the model are fairly different from the ones observed in the market. Cartea and Villaplana (2008) [17] build a model for wholesale power prices explained by two state variables (demand and capacity) and calculate the forward premium. Writing D, C for the demand and capacity, they model D and C by = fD(t)+XP, C = f0(t)+X, where fD, fc are deterministic functions, and X, X° are independent Ornstein-Uhlenbeck processes. They then take the spot price as given by St exp{D +‘yG}. They perform empirical research embracing PJM, England and Wales, and Nord Pool markets. They find that, depending on the market and the period under study, the volatility of capacity and the market price of capacity risk could either put upward or downward pressure on forward prices. They also 15 Chapter 1. Introduction find that the forward premium follows a seasonal pattern, being positive in the months of high volatility of demand and close to zero or even negative in the months of low volatility of demand. Inspired by Cartea and Villaplana (2008), Lyle and Elliott (2009) [60] present an hybrid model that uses a supply-demand approach for price elec tricity derivatives. They assume that the system demand D(t) is given by D(t) = f(t) + b(t) where .b(t) is an Ornstein-Uhlenbeck process and f(t) a deterministic func tion. For the supply side, they model the curve S(t, P) which gives the supply at time t if the price is P. They consider curves of the form S(t, P) = aSb(t) + blog(cP + ) where Sb(t) is the base portion of the system supply and a, b, c, and are positive constants. They consider two different models for Sb(t): a mean- reversing model, and a Markov chain model. The equilibrium price is given by P(t) = (exp {_aSb(t) - D(t) }) Using these equations Lyle and Elliott are able to obtain closed-form solutions for European options. They test the model on Alberta prices data calculating the first four empirical moments. The model gives a good fit for the mean and standard deviation but not for the skewness and kurtosis. Benth et al. (2007) [7] propose a non-Gaussian Ornstein-Uhlenbeck pro cess which takes into account seasonality and price spikes. Their model is S(t) = h(t) + X(t) where S(t) is the spot price, h(t) is a deterministic periodic function and X(t) is a sum of independent Levy-driven Ornstein-Uhlenbeck components. X(t) d}4(t) = —cY(t)dt + u(t)dL(t), Y(0) yj. 16 Chapter 1. Introduction Here the processes L(t) are independent, possibly time inhomogeneous pure jump Levy processes with E(L(1)) < oc. L(t) can be written in terms of their jump measures N2(dt, dz), i = 1, ..., n. ft foo Li(t)=J J zN(ds,dz).00 The deterministic predictable compensator of N(ds, dz) which is the jump measure of L(t) is of the form: v(dt,dz) = p(t)dtv(t)dz. Here p(t) is a deterministic function that contrOls seasonal variation of the jump intensity, u(t) controls the seasonal variation of the jump sizes, cx is the level of mean reversion, and L(t) controls the variation of price such as the daily volatile variation and price spikes. Benth et al. (2007) provide closed- form and semi-closed form solutions for forwards and options on forwards. This model, coupled with a good description of price seasonality, provides a precise characterization of electricity spot price behavior. In addition, due to its arithmetic structure, it is analytically tractable when it comes to futures and other derivatives pricing. Although the model seems to capture the stylized facts of spot prices market such as mean reversion, seasonality and price spikes, the authors did not make a precise statistical analysis of the quality of the model. However, they suggested the particle filter as a possible method to estimate the parameters in the model. Parameter estimation for this model would appear to be a significant challenge. Hikspoors and Jaimungal (2007) [44] consider two models for oil prices. The first is a two factor version of the model of Schwartz: = exp{h(t) + X}, dX = ).x(Y — X)dt + crxdWt, (1.9) d = )y(b — Y)dt + ydZ, dWdZ pdt. (In fact, in order to value spread options, they consider two different com modities satisfying (1.9)). The second is a modification of (1.9) with an 17 Chapter 1. Introduction additional jump component to handle price spikes: = exp{h(t) + Xt + J}, (1.10) where X, is given by (1.9), and dJ = —,cJt-dt+dQt (1.11) with Q a compound Poisson process, and t denotes the instant immediately before time t. Through measure changes induced by a pseudo-numeraire, they obtained, for both models, exchange option and futures prices in closed form both under real-world and risk-neutral measures. Also they consider the problem of model calibration. For the jump model (1.10) they suggest a modification of the procedure of [20]. This is to identify the price jump by searching for days with a price change more than 3 times the standard deviation of the daily price change. (This procedure is then run several times to produce a ‘despiked’ series). They do not apply this method to electricity prices, but do estimate parameters for their first model for oil, from a period of about 3 years of data. They consider parameters built under the real world probability IP, and a risk neutral probability Q. Let us write a bar on the model parameters to denote parameters with respect to Q. If P Q on the filtration generated by both X and , then one has 0XJX, OyJy, pp. The parameter estimation proceeded in a number of steps: 1. Using least squares, the risk neutral parameters were estimated from the forward data, as is the unknown process Y. 2. Given (Xe, Y), the remaining real world parameters Ax, Ay, were estimated by regression. One surprising feature of the estimates is that they obtain p = —0.96, so that the long-term process adds little extra randomness. The advantage of two or more factor models is that they allow for a good mathematical description of the problem. Furthermore they have a better fit to historical data and provide a better relation between spot and futures prices, see [73]. 18 Chapter 1. Introduction Many of the models in the literature incorporate ‘regime switching’ com ponents. For example, the basic model of Nomikos and Soldatos (2008) [64] for Nord Pool prices, is similar to that of Hikspoors and Jaimungal: S =h(t)+exp{Xt+}, (1.12) where X is an Ornstein-Uhlenbeck type process, and Y, is a jump process intended to take account of price spikes. However, both X and incorporate regime dependent terms. X, is given by dX = — X)dt + uxdW, (1.13) where R is the water reservoir level, and is assumed to follow a two-state Markov chain with state space {wet, dry}. (Much of the Nord Pool electricity production is by hydro.) is driven by a ‘jump process’, and satisfies dY = ,c2Ytdt + dLi. (1.14) Here L is a jump process, with rates and jump distribution dependent on the season. Davison et al. (2002) [22] they propose a hybrid model based on the ratio a(t) between demand and capacity. At each time step t the spot price S is drawn from a distribution P(t) which is a mixture of Gaussian distributions given by P(t) = (1- t)))PL(t) + f(c(t))PH(t). Here PH(t) is the price-spike distribution, PL(t) is the low-price distribution, and is is a function of c that represents a relative demand-capacity ratio. The distribution e plays the role of a switching variable that determines whether the price is to be drawn from PL(t) or PH(t), i.e. the probability of a spike. They assumed (t) = tanh(20 * (a(t) — 0.85)) + where the constants are determined by historical PJM spot prices. The distributions PL and PH are taken to be Gaussian and the function (t) is deterministic. It appears that each time step t independent samples are taken from the distribution PL, PH. The choice of independent samples from PL would 19 Chapter 1. Introduction lead to highly oscillation prices, and in Anderson and Davison (2008) [1] they replace FL by a Brownian motion. To test the model, they simulated trajectories and compare statistical moments. Applying the Kolmogorov Smirnov test, they concluded that their model is able to simulate data that are from a similar distribution to the observed prices. Eriwein et al. (2008) [29] develop and analyze an exponential Ornstein Uhienbeck process with an added jump process based on hidden Markov model (HMM) setting.The jump component is a Poisson process where the mean and variance are controlled by a discrete time HMM. That is, the spot price, that is partially observed (the underlying economic state is hidden) is given by S = D exp{Xt} where D is a deterministic function and dX = a(z)(/3(zt) — X)dt + a(zt)dW + JdN. Here Zt is a Markov chain with 2 or 3 states, and the jump sizes J are conditional Gaussian distributed, i.e. JIZt N((zt),u2(zt) . They apply the EM algorithm to estimate the parameters in the model using data from the Nord Pool market. The model captures some of the spikes presents in the real data for the 2 and 3-state Markov chain. A puzzling feature of this paper is that the estimates for the transition matrix of zt suggest that Zt are close to i.i.d. random variables. In the same theoretical framework a continuous-time process is derived by Kholodnyi (2001) [52], where self-reversing non-Markovian spikes are added to a Markovian regular price process. One sees in the literature the need to balance two competing demands. Simple models, particularly simple one factor models, have relatively few pa rameters, and these parameters may be relatively easy to estimate. However, these models generally fail to capture one or more features of real markets of which the most difficult are the existence of price spikes, and the relation between spot and future prices. 20 Chapter 1. Introduction The need for a better fit with data leads to more complicated models, usually with hidden variables, and sometimes with multiple regimes. How ever introducing more factors requires introducing more parameters into the model. Parameter estimation then becomes a significant challenge. Benth et al. (2008) [6] remark: “The question of estimating such models on data is not an easy one... For multi-factor models this may be an even more challenging problem, involving highly sophisticated estimation techniques”. Such parameter estimation is, of course, an essential preliminary to the valuation of options or derivatives based on the commodity. Table 1.1 sum marize some of the statistical models. Although the standard statistical procedure for estimation of a partly unobserved process involves filtering methods, few of the papers in the liter ature use those techniques. An exception is Culot et al. (2006) [21]. They consider a model of the form log St = h(t) + + yTX, where h(t) is deterministic, and X, and X are spikes and long-term fac tors. The spike process X is an m state Markov regime switching pro cess, while X is a 3-dimensional Ornstein-Uhlenbeck type process given by = (X1),X2),X3)), where dX = —kdt + udW, and W are independent Brownian motions. After estimating the jump term X and subtracting this from the series, the authors combine Kalman filtering techniques with maximum likelihood estimation, using both spot and forward prices, to estimate the parameters for X. Another exception is Kellerhals (2001) [51]. He developed a model for short-term electricity forwards. The suggested stochastic volatility model uses the non-tradeable spot price St of electricity and its variance rate Vt as state variables. The stochastic specification of the processes is given by dS = Sdt + S dvt = udt + a/dZ, dWdZ = pdt. 21 Chapter 1. Introduction [ SPOT PRICE BASED MODELS Authors Model Specification Gibson / dS = (i — S)Sdt + aiStdzi Schwartz (1990) = — ö)dt + cx2dz dz12 = pdt Schwartz (1997) dS = (rt — 6)Sdt + uiStdzi dö = — ö)dt + cr2dz drt a(m — rt)dt +j3dz dz12 = pidt,__dz2dz3=p2dt,__dz1z3= p3dt Lucia/ S=h(t)+X+Y Schwartz (2002) dX = —,\Xdt + cxx dWx dY, = pdt + crydWy dWxdWy = pdt 1 (1 + X)’1, 1 + crX >Barlow (2002) St c 1/crL. o 1+aXteo dX = —)(Xt — a) dt + crdWt Villaplana (2003) ln St h(t) + X, + dX = —icxXtdt + cxxdWi + JdN\) — JddN2d) dY = —Icy(p — Yt)dt + aydW2 dW12= pdt Geman / S(t) = exp{E(t)} Roncoroni (2006) dE(t) = [h(t) + O((t) — E(tj)]dt + crdW(t) + f(E(tj)dJ(t) Hikspoors / St = exp{h(t) + Xt + J} Jaimungal (2007) dX = ).x(Yt — X)dt + cxxdWt dY = \y(q$ — Y)dt + crydZ dJ_=_—,cJt-dt_+_dQt Benth / S(t) = h(t) + X(t) et al. (2007) X(t) (t) dY(t)_=_—)ç’Yj(t)dt_+_dL2(t),__Y(0) = Table 1.1: Models for electricity prices. Using maximum likelihood estimation based on Kalman filtering he reports empirical results on electricity data from the California market. Given observations which derive from hidden state process u1, ‘u2, ..., one may distinguish ‘online’ methods from ‘batch’ methods. The online methods provide an estimate ii which can easily be updated given an addi tional observation vt+i; while the batch methods estimate the whole series 22 Chapter 1. Introduction (ui, ...,Ut) from (vi, ...,Vt). Following most of the literature in this thesis, we have used online meth ods (such as the Kalman filter and particle filter), rather than batch methods such as the EM algorithm. One motivation for doing this is the need of mar ket participants to update their models in real time. As is clear from the survey above, many authors have attempted to de sign models which capture the typical properties of electricity prices, namely seasonalities, spikes, and stochastic mean-reversion, but none of the models proposed so far has commanded wide assent. In this thesis we consider some relatively simple multi-factor models, with a relatively small parameters. A second emphasis is the use of both spot and forward prices for parameter estimation, and a third is the use of filtering techniques to estimate the hidden processes, and hence the model parameters. We will introduce three different spot price models from which we can also extract the futures prices. All of these models are capable of capturing some of the features of the spot price dynamics and imply certain dynamics for futures prices. The first model (MROU model) is a Gaussian two-factor model where the spot prices is an mean reversion process which reverts to a stochastic mean, also fluctuating as Ornstein-Uhlenbeck process. Since the presence of spikes is a fundamental feature of electricity prices, and any relevant spot price model should take this feature into account, we introduce a second model, an extension of the MROU model with a jump component (MROU with jumps). However, the inclusion of the jump component introduces two kinds of problems in parameter estimation. The first is that the inclusion of jumps adds several new parameters, to describe the jump frequency and distribu tion. The second is that the jump models are non-Gaussian, and the best known filtering technique for these models, the particle filter, is not easily adapted to handle parameter estimation. To avoid this problem, we introduce the third model (NLMROU model) based on the MROU model that produce spikes introducing only one more parameter. 23 Chapter 1. Introduction In general terms, a statistical model is good if it successfully captures the main features of the observed data. Various statistical tests (such as the Kolmogorov-Smirnov goodness of fit) can be used to test and compare statistical models. From this viewpoint, the theory of electricity prices is rather undeveloped. No systematic comparison of the various models in the literature has been made. One reason for this it that in many cases authors have proposed models, but have not yet developed techniques for parameter estimation. In this thesis, we have followed other workers in this area in using fairly simple tests for our models. We have compared moments and sample paths of simulated and real data. Even these simple tests indicate that our models do not capture all the features of real prices. 24 Chapter 2 Filtering 2.1 State space formulation A state space model is a representation of the joint dynamic evolution of an observable random vector Vt and a unobservable state vector Ut. It is based on two important sets of system equations: the measurement equation and the transition equation. The transition equation describes the evolution of the state vector and the measurement equation reflects how the state interacts with the vector of observations. The evolution of the state is assumed to be autonomous, that is it does not depend on the measurement equation. We consider the non-linear and non-Gaussian state-space model, which is represented in the following general form: For t 1, 2, ..., ftQut_i, qt—i) (transition equation), (2.1) Vt = h(u, rt) (measurement equation), (2.2) where ft : Rn” x Thq H-> IR’ and h : R x R i—* are vector functions, which are assumed to be known and possibly non-linear. The process and measurement noises, qt and rt respectively, are independent with known but arbitrary densities. In addition, we assume that the initial distribution of ‘u0 is available, that is pQuo) := p(UOIVO). Associated with a state-space model is the problem of estimating the unobservable state using a set of observations. To do so, from a Bayesian perspective we need to estimate the filtering density p(utIVit), where Vi,t {V1,V2, ..., Vt} is the past history of the observed process up to time t. If possible, we wish to do this recursively, so that, p(utlvi:t) can be calculated by updating the estimate p(Utvit_i) with the new observation Vt. The esti 25 Chapter 2. Filtering mate of the filtering density can be obtained in two stages (prediction and updating) as follows. Applying the Chapman-Kolrnogorov equation we can write the time up date iteration as: p(utlvi.t_i) Vi:t_i)P(Ut_iIVi.t_i)dUt_i = fP(uut_i)P(ut_lIVl:t_l)dut_i (2.3) by using in the last equation the Markov property. Now, after the observation Vt 1S available we use the Bayes rule to have p(VitIU)p(ut) p(utIVi;t) = / PPi:t — p(Vt, V1:t_i Iu)p(u) — p(vt,V:t_i) = p(Vt IV1:t_i, ut)p(Vi:t_i Iut)p(ut) P(VtIVI:t_i)p(Vi:t_i) = p(VtI Vi:t_i, Ut)P(UtIV1:t_i)P(Vi:t_i)P(Ut) p(VtIV1t_1)p(V1:t_1)p(ut) — p(VtInt)p(utIVit_i) 2 4) — p(VtIVIt_i) where the denominator could be written as = J p(vtlut)p(utlvit_i)dut. (2.5) Unfortunately, in general there do not exist closed-form expressions for equations (2.3) and (2.4). The main exception to this is where (2.1) and (2.2) are linear and the noise processes qt, rt are Gaussian, and in this case the solution is given by the Kalman Filter [48]. 26 Chapter 2. Filtering 2.2 The Kalman filter Equations (2.1) and (2.2) reduce to the following special case where a linear Gaussian state-space model is considered. To include a more general case we have included two additive components C, and A to the transition and measurement equation respectively. We have Ut = C + D Uti + q_ , (2.6) (nx1) (nxn,) (nx1) Vt = A + B u + rt , (2.7) (nx1) (n,xn,,) (mx1) where the process and measurement noises are normally distributed (:) N(() (q )). (2.8) The initial condition becomes no ‘—‘ NQZZo,o) (2.9) and the matrices G, D, A, B, q, and r are assumed to be known. N(p ) denotes a Gaussian density with mean p. and covariance , that is: N(p., ) := 2I_1/2 exp{—(x — p.)’’(x — p.)}. (2.10) Here . denotes the determinant. For the model above (2.6)-(2.9), it follows that the transition density p(ut+iut) and the measurement density p(VtI’ut) are normal. It can be shown that this implies that also the prediction and filtering densities are normal, see [77] for details. We have p(ut Ivi:t_i) = N(u1_, (2.11) p(’ut Vt) = ]V(’ut1t, (2.12) p(vtIVit_1) = ]VQOtit.....i,F1_), (2.13) 27 Chapter 2. Filtering where the conditional means tIti, ut(t, tit_i and conditional covariances tit_, >tit and F1_ are computed by the following pseudo code of the Kalman filter. Algorithm 1 (Kalman filter) • Step 1, Initialization Set ‘U010 = , = z0, and set t = 1. • Step 2, Prediction Compute = D D + • Step 3, Innovation Define et = Vt — tt_i (2.14) with = A + B u_i (2.15) and compute F1_ = B B + (2.16) • Step 4, Updating Compute K = t_i B F’1, (2.17) = u-i + Dtt_i B F’1 Ct, (2.18) = — B F’ B (2.19) • Step 5, Looping if t <n, set t = t + 1 and go to Step 2; else stop. 28 Chapter 2. Filtering In the linear Gaussian case, the Kalman filter has strong optimality prop erties. (2.11)-(2.13) give the maximum likelihood estimator of u given Vit, and this is also the minimum mean square error estimator. That is, in the mean square sense no other algorithm can perform better than the Kalman filter in the Gaussian environment; see [40] for details. In cases where the measurement or the transition equation are nonlinear, sub-optimal solutions such as the extended Kalman filter (EKF) and the un scented Kalman filter (EKF) are commonly used to solve the problem. The extended Kalman filter simply linearize all nonlinear transformations and substitutes a Jacobian matrix for the linear transformations in the Kalman filter equations [41]. Although it is easy to implement it has a number of limitations especially if the system nonlinearities are severe or the true dis tribution is multimodal or highly skewed - see [75]. The unscented Kalman filter gives a more accurate performance for nonlinear equations, but does rely on the noise being Gaussian. 2.3 The unscented Kalman filter An alternative filter with performance superior to the extended Kalinan filter is the unscented Kalman filter. Unlike the extended Kalman filter it does not approximate the nonlin ear function of the process and the observation, it uses the true nonlinear models to approximate the distribution of the state variable ut by applying an ‘unscented transformation to it. The unscented transformation uses the so-called sigma-points that capture the mean and covariance of the original distributions and, when propagated through the true nonlinear system, cap ture the posterior mean and covariance accurately to third order. For more details see [47, 75]. Unlike the particle filter considered in the next section, which requires a large number of points, the unscented transform only re quires 2n + 1 points to capture the mean and covariance of a probability distribution in ll’. Let us consider a simplified version of the UKF formulation, where we as 29 Chapter 2. Filtering sume that both the transition and measurement noises are additive Gaussian, that is, Ut = f(ut_i)+qt_i, (2.20) Vt h(n) + rt, (2.21) where qt—i N(0,q) and rt ‘- N(O,). Here Ut e R, Vt E R. The algorithm can be described as follows: Algorithm 2 (Unscented Kalman filter) • Step 1, Initialization Set = E(uo) and o E[(uo — i7o)(uo —z)’]. Set t = 1. • Step 2, Unscented transformation Compute the sigma points and weights. i = 0 Xt_i(O) = (m) ,\Wo = i = 1,...,n Xt_1(i) flt_i + (V(n+)t_1), (m) (c) 1 — — 2(n+)’ i = n + 1,..., 2n Xt-i(i) = flt_i — (Vn + (m) (c)_ 1 — — 2(n+A) Here the subscripts i and i—n correspond to the jth and j_flth columns of the square-root matrix (Cholesky factorization). ..\ = a2(n + ic) — n is a scaling parameter. a determines the spread of the sigma points around and is usually set to a small positive value. frC is a secondary 30 Chapter 2. Filtering scaling parameter which is usually set to 0, and is used to incorporate prior knowledge of the distribution of Step 3, Time update Fori=O,1,...,2n, XtIt_i(i) — f(t_i(i)) Vtt_i (i) = h(x1_(i)) 2n — (m) UtIt_1 W Xtjt—i (%) = Wc(xtt_, (i) — (XtIt-i (i) — tIt-1)’ + q tt-i = Wm)VI_1(j) (2.22) w ( VI_1 (i) — i-i) (Vtit_i (z) — Vtit-i)’ + (2.23) w (xtit-1(i) — (Vt1t_(1) —-1)’. • Step 4, Measurement update Calculate K = (2.24) Ut = tt_i + K(v (2.25) Zt = — (2.26) • Step 5, Looping if t < n, set t = t + 1, update Ut and > and go to Step 2; else stop. 31 Chapter 2. Filtering 2.4 Particle filter A different approach to filtering has recently become popular [26, 79]. In this approach, we use Monte Carlo simulations instead of Gaussian approxima tions for p(utlvt), as in the Kalman filter. This method allows for a complete representation of the filtering distribution, so that any statistical estimate can be easily calculated. This filter has the advantage that it allows one to deal with fundamentally non-Gaussian situations. The idea is based on the importance sampling technique [70, 76]. Recall that UOt = {uo,ui, ...,ut} is the (unknown) true state of the system, and {vi, v2, ..., Vt} are the observations. Suppose we wish to calculate the conditional expectation E(f(uo.t) IUi:t) = f f(ü:t)p(0: Vi:t)dU0: . (2.27) Inmost non-Gaussian or non-linear situations, the true distribution p(uotvit) will be impossible to calculate. Importance sampling works by instead sam pling from a proposal distribution q(uotlvit), which can be easily sampled from. The support of q(uo:tIvit) is assumed to cover that of p(uotvit). We can write E(f(uot) Ivi:t) = f f(UO:t) 0:tl q(uotvmt)duot, f p(vi.tI’uo.t)pQuo.t) = i f(ot) q(uo:tvi:t)duo:t, J pVi:t)qtUO;t V1.t) = f f(uo:t) W(UO:t) q(uotIvit)duot (2.28)p(vi;t) where W(UO:t) = p(vi.tluo.t)p(uo.t) (2.29)q(notlvit) is defined as the filtering non-normalized weight at step t. Now 32 Chapter 2. Filtering E(f(uo:t)Ivi.t) = f f(uot)w(uo:t)q(uo;tvit)duo:tp(vit) — f f(uo:t)w(uot)q(uotIvi:t)duot q(uo.tIvi.,)I p(vi;t uo:t)p(uo:t) q(uo.t Ivi.i) du0, — f f(uot)w(uot)q(uo:tvit)duot — J W (U:t) q(uo;t,ji. ) duo — Eq(f(uo:t)w(uo:t)) — Eq(w(uo:t)) = Eq(f(uo:t)’th(uo:t)), where = Eq(w(uo:t)) (2.30) is defined to be the filtering normalized weight at step t. Now let u, i = 1, 2, ..., n7, be a Monte-Carlo sample from q(uotvit). Then by the law of large numbers E(f(uo.t) Ivi:t) f(())) (2.31) where now i) = . (2.32) w(u) Note that (2.31) is correctly normalized: if f 1 then both sides of (2.31) equal 1. 33 Chapter 2. Filtering Thus, provided we can sample from the proposal distribution q(UotIvi;t), and calculate the weights W(UO:t) given by (2.29), then we can estimate E(f(uo:t)Ivi:t). Note also, that if we can only calculate w(.) to a multi plicative constant c, which can depend on Vi:t, then this constant drops out when we calculate z7 in (2.32), and so the procedure will still work. Suppose that the importance density is chosen to factorize, so that q(uo:tIvi:t) = q(utuot_ivit)q(uot_iIvit_i). (2.33) Then (i) (i) / (i)\ — p(Vit1U0:t)P(.t WUo.t) — (i) q(u0. Vit) (i) (i) (i) (i) (i) = p(vt, V1:t_1UO.t_i, u )P(ut Uo:t_i)P(Uo.t_i) q(r4 Iz4_i, vi:t)q(U_1IVi:t_i) The measurement equation (2.2) implies that Vt depends on O:t only through t and that v1,v2, ..., ‘Vt are conditionally independent given ‘UO:t. Therefore (i) (i) (i) (i) p(Vt, Vi:t_1 IUO:t_i,Ut ) = p(VtJU )p(Vi:t_i 1U0:t_i) Further, (2.1) implies that Ut is lVlarkov, so p(UIU1)= p(uU1). Combining these eqilations, (j) = p(Vi;t1 Iui)P(u i)p(VtIu)p(r4IUi) q(u_1v1;_)q(4tIu_i, Vi:t) — (i) p(vt )p(u Iu1) —w_1 (i) (i) ( . q(U UO:t_i, Vit) 34 Chapter 2. Filtering Thus, the importance weights defined in (2.34) can be updated in a simple way at each time step. (2.31) implies that the filtering density p(UtIVit) can be approximated by p(utlvit) (ui) (2.35) where S() denote the point mass at u. The success of this operation de pends on how close the proposal distribution is to the posterior and whether the resulting point-mass approximation is an adequate representation of the distribution of interest. Although sequential importance sampling poses only one restriction on the importance density, equation (2.33), with the number of choices being unlimited, the design of the appropriate proposal function is, in fact, one of the most critical issues in importance sampling algorithms [79]. Poor choice leads to poor approximation in (2.35), and to poor algo rithm performance in general. See [26, 69] for more details and variants of the particle filter. One major problem with this algorithm is that the variance of the weights increases steadily over time. If one starts with a fixed number n,., of particles, then in practice after a while nearly all the mass of the distribution in (2.35) is concentrated at one particle. Not surprisingly, this leads to poor algorithm performance. In order to solve this, we resample the points to create copies of particles with large importance weights and to remove those with negligible importance weights. This ensures that there are sufficient particles exploring regions of high probability in the next time step [37]. Various methods have been suggested for this [2, 13, 58]. In the particle filter literature four basic resampling algorithms can be identified: 1. Multinomial resampling Generate n, ordered uniform random numbers a = aj+i, = with — U[O, 1) and use them to select the new particles (1) according to the multino 35 Chapter 2. Filtering mial distribution. That is, — (i) with i s.t. a E s) s)) where F’ denotes the generalized inverse of the cumulative probability distribution of the normalized particle weights. 2. Stratified resampling Generate n ordered random numbers a= (i—i)+ witha1U[O,1) rip and use them to select () according to the multinomial distribution. 3. Systematic resampling Generate ni,, ordered numbers (j—1)+a with a U[O, 1) np and use them to select according to the multinomial distribution. 4. Residual resampling Allocate n = copies of particle to the new distribution. Additionally, resample rn = n, — n particles from {u) } by mak ing n,’ copies of particle where the probability for selecting n is (i) /proportional to w = rlw — n2 using one of the resampling schemes mentioned earlier. A illustration of generic particle filter is shown in Figure 2.1. The whole particle filter algorithm can be implemented in the following way: 36 Chapter 2. Filtering i1 rplO particles 0 . (u 1} * np A / u, . • • • • (u .w Figure 2.1: A graphical representation of the particle filter with importance sampling and resampling. Algorithm 3 (Generic Particle Filter) Step 1, Initialization For time step t = 0, choose u0 and for each i between 1 and n (number of particles), take p(uo) where p(. ) is the initial distribution. Also take (importance weights). np While 1 < t < n (number of observations) • Step 2, Prediction For each index i sample q(utlu2i,vit). 37 Chapter 2. Filtering Step 3, Importance sampling Calculate the probabilities p(vt u) (likelihood distribution), Iu2) (prior distribution), (2 36) q u21,Vit) (proposal distribution), and the associated weights for each i (i) (i) p(vt Iu )p(u u1) wt = ‘wt_i (i) (i)q( IUt_i,Vi:t) • Step 4, Normalizing Normalize the weights (i) _____ wt Z,;=lWt • Step 5, Resampling Resample the points z4 and reset w = = 1/np. • Step 6, Looping Increment t and go back to Step 2. Stop at the end of the While ioop. 2.5 Parameter estimation via maximum likelihood Up to this point we have assumed that the functions f and h, with the distribution of the noise qt, Tt are fully known. We have discussed the filtering problem - that is how to estimate Uot from observations V1t. But in many applications, and in particular nearly all financial applications f, h, q,, and rt, will depend on unknown parameters. In this case the structure of the nonlinear state-space model becomes 38 Chapter 2. Filtering ut = ft(ut_i, qt_i, 61) (transition equation), (2.37) Vt = ht(ut, Tt, 61) (measurement equation), (2.38) where 0 e e C Rn0 denotes the parameters in the model. Given the above structure, in this section we addresses the problem of estimating the parameters 0 from the observed data. The parameter estima tion problem for state-space models has generated a lot of interest over the past few years and many techniques have been proposed to solve it. These methods could be broadly classified as Maximum Likelihood or Bayesian. Using a Maximum Likelihood formulation the estimate of 0 is the maxi mizing argument of the likelihood of the observed data, i.e. =argmaxp9(Vi,v2,...,Vt) (2.39) OEO where po(’ui, V2, ..., ‘Vt) denotes the joint density of the observations up to time t. In a more convenient form we can rewrite (2.39) as 0= argmax J(Vi:t), Co(Vi.) = log po(Vi:t) (2.40) oee where V1t = {v, ‘v2, ..., ‘Vt}. The joint density can be written as the product of the conditional densities: pe(Vi:t) = fJP9(VkIV1:k_1), (2.41) where pe(viJVo) = po(Vi). Thus the log-likelihood function is = logpo(VkIV1k_1). (2.42) The material presented up to here has dealt with the state-space model using a quite general formulation. The algorithm described above is in prin ciple the full answer to the problem of parameter estimation. Now, from 39 Chapter 2. Filtering a classical approach we can use some numerical optimization search proce dure (Newton’s method, Nelder-Mead) on (2.42) in order to calculate the maximum likelihood estimate 8. Maximum likelihood estimation (MLE) of 6 is particularly simple in the linear Gaussian state-space model (2.6)-(2.7), since the density func tion Po(VkIVl:k_1) is the normal distribution with mean VkIIc-1 and covariance matrix FkIk_1 given by equations (2.15) and (2.16) respectively. Thus, po(vkIvl:k_1) = I(2)IFkIk_1II” exp{—(vt — kIk_1)’F_l(vt — VkIk_1)} and the log-likelihood function becomes £O(Vi:t) — [log IFkIk_1I + (Vk VkIk_1)Fl(Vk — i’)] (2.43) where FkIk_1I denotes the determinant of FkIk_1. Thus finding the MLE is quite straightforward for the Kalman filter. In the general nonlinear case this approach is non-trivial, since the distri bution Po (v1,v2, ..., Vt) 1S generally not available in closed form. However, if we are using the UKF algorithm then, one can approximate the true (non Gaus sian) distribution Po(VkIV1.k_1) by a Gaussian distribution NQikIk_1, VkVk), where vkIkl and Zvkvk are given by equations (2.22) and (2.23). Hence one obtains — [log IVkVkI + (Vk — VkIk_1)’vk(Vk — Since the UKF algorithm is based on the equations of a Kalman filter, the maximization of the log-likelihood can be done exactly as in the Kalman filter. For the general nonlinear case, with non Gaussian noise, we have seen that the particle filter provides a technique for filtering with known parameters. In 40 Chapter 2. Filtering general, we can approximate (2.42) using equation (2.5). Given the likelihood at step k, po(VkV1:k_1) = fPo(vkIuk)Po(ukIvl:k_1)duk this could be written as I Po(ukIVl:k_1) PO(VkI’01:k_1) = i e(v’I’) qe(ukIuk_1,v1:k)duk,qoukuk_1,v1:k) and given that by construction the u’s are distributed according to q, we can write the Monte Carlo approximation po(vkIvlk_1) ± (2.44) Thus, we can estimate £(&) by slog (± Ew). (2.45) k=1 ‘ j=1 While this does give an approximation to the log-likelihood, this approach has several problems if we try to use it to obtain the MLE. For a fixed 9 equation (2.45) gives a random variable, where the randomness comes from the particle filter. See Chapter 4 below for an account of the difficulties this causes. An alternative approach to maximize the likelihood is employed the Ex pectation Maximization (EM) algorithm [74]. The objective of the algorithm is to maximize the likelihood of the observed data (2.42) in the presence of the hidden variables (Uot = {uo,u1,.. . , Ut}). The basic idea is that if we could observe UOt, in addition to the observations Vi,t then we would con sider {UOt, Vi:t} as the complete data with the joint density po(uo:t,vi:t) po(uo)llpo(ukIuk_1)llpo(vkuk) (2.46) 41 Chapter 2. Filtering and seek the maximum log-likelihood estimate of & via 6= argmax £O(Vi:t, Uo:t), Q(Vi:t, uot) logpo(vi:t, flo:t) (2.47) oEe The EM algorithm for maximizing L(Vit, u) is a two step procedure. 1. (E-Step) Computes the expected value of £0(Vi:t, flit) over the hidden (missing) data i:t based on the current value of the parameters 6C) and the observations V1:t Q(616(i)) = E(L0(Vi:t, UO:t)IV1:t, e()) = flogp (UOt, v1:t)Pe(,) (UO:t IVi:t) 2. (M-Step) Update the parameter estimate maximizing Q(&j6()) with respect to 6, argmax Q(&1& )) (2.48) 0 and we repeat this two-step process until a fixed point of Q is obtained. Unfortunately, there are very few situations where an exact and tractable solution exists for these two steps. One exception is the linear Gaussian state- space. The algorithm is described in [71]. Work applying the EM algoi’ithin to nonlinear dynamical systems in the form of (2.20) and (2.21) is reported in [23, 34, 41]. To use the EM algorithm for a general nonlinear state-space the sequential Monte Carlo (particle filter) methods are employed to approximate Q(&1& )) numerically. For the maximization step there is no standard method to solve this problem, and so it is necessary to proceed on case-by-case basis. In one general approach we calculate gradients with respect to 6, and use a gradient-based search procedure to find the maximum. Recent works that use this technique are in [36, 65, 81]. None of these methods is simple. 42 Chapter 2. Filtering 2.6 Parameter estimation via Bayesian methods Because of these difficulties one would like an alternative to the MLE approx imation. One Bayesian approach, described in ([43, 53, 57]) is to consider a Bayesian estimation by concatenating the state vector u with the unknown parameter 0, and introduce an artificial dynamic on the parameter. That is, we replace 0 by Ot and define a new state vector Yt (Ut) (2.49) where Ot - p(Ot—i Vit). Then one applies the particle filter to this augmented state-space. However this method has a number of difficulties and problems [50]. At a theoretical level it is not altogether satisfactory to replace a fixed parameter 0 by a random evolution process Ot. In addition, there are various practical problems associated with the choice of the artificial dynamic of 0. Thirdly, there are problems with the performance of the algorithm. These difficulties have been discussed in the literature, for example, in [18, 50]. Under the Bayesian framework, more sophisticated and new methods are being proposal. See [50] for an overview of particle filters methods for param eter estimation considering a nonlinear non-Gaussian state-space models. 43 Chapter 3 MROU model In this chapter we present a Gaussian two-factor model known as Mean- Reverting to Ornstein-Uhlenbeck model (MROU) for the spot price and the convenience yield that captures some of the characteristics that we had de scribed above of the power market and the dynamics of the futures prices. 3.1 Double mean-reversion model There are many parallels between interest-rate models and modeling com modity prices, so many models originally developed for stock and interest rate markets have been applied to the energy market. We implement, for the valuation of electricity futures contracts, a two-factor mean-reverting model originally proposed in [5], and considered previously in [44] for oil commodity prices. Let (p2, F, {}>o Q) be a complete filtered probability space where Q is the risk-neutral measure. If S is the spot price then = exp{Xt + h(t)}, (3.1) dX —)‘x(X — L) dt + x dW’, (3.2) dL = —AL(Lt — L) dt + L dW2. (3.3) Here X, is the observed deseasonal log spot price, and L is a non-observed long-term mean process. We assume that both processes are given under the risk-neutral measure and the two Brownian motions, W’ and TT/2 satisfy d(T’V1,W?) = pdt. h(t) is an arbitrary deterministic function that accounts for seasonality. The difference between this model and the Gibson and Schwartz model 44 Chapter 3. MROU model is that here, both the log spot price and the convenience yield follow an Ornstein-Uhlenbeck type process. The seasonal component h(t) combines trend and seasonality. Usually it consists of sum of sinusoidal functions which incorporate predictable daily and annual periodicity and dummy variables which incorporate predictable workday/weekend and holiday effects. In this thesis, we consider a sum of two cosine functions with distinct periods with a linear trend, that is h(t) = + 0t + cos (Ti +2nt) (3.4) where m represent the seasonality period and the parameter = needs to be estimated. The seasonalities have been discussed extensively in the literature, see for example Lucia and Schwartz [59], Cartea and Figueroa [16], and Benth et al. [7]. Although there are several ways of deseasonalising the data, usually it is estimated by means of non-linear regression methods. That is, if t = t1, t2, ..., t we estimate the seasonality function by fitting h(t) to the log- prices using least squares estimation 2 = argmin (i(tj) — log(Sti)) ,2,rl ,T2) Write (t) for the estimate. The deseasonalized log-spot price is given by: = logS .-(t). 3.2 Radon-Nikodym theorem for Ornstein-Uhlenbeck processes Before we continue with parameter estimation for the MROU model, we consider a simple case, that is the problem of parameter estimation for the Ornstein-Uhlenbeck process OU(A, a, a) defined by dX = —A(X — a)dt + crdWt. (3.5) 45 Chapter 3. MROU model If one observes the whole process x = {X, 0 t T} the parameter o can be estimated exactly, using the quadratic variation of X, but ) and a cannot. The simplest statistical problem is when one has two alternatives: H0: a=a0,) A vs. H1: a=a1,) ) Denoting Po, P for the probability measures for (3.5) associated with H0 and H1, the Radon-Nikodym theorem gives the optimal test for this in terms of the likelihood ratio. If — L(xIHi) — dP1 T L(xlHo) dIP0 then the test takes the form Reject H0 if ZT c(a), where c(cl) is given by IP0(reject H0) = IPo(ZT c(c)) = c. The power of the test is given by p = IP1(accept H1) = IP1(ZT c(a)). Using the result in [44] (Thm. 3.2) we obtain log(ZT) = = +c1X)dW3— fco + ciX5)2ds. The constants c0, and c1 satisfy — Ao — c0= , c1= . (3.6) ci ci We used Monte Carlo simulation to find the power of the test for these sample tests - see Table 3.1. 46 Chapter 3. MROU model H-J 1 N-S 1 N-S 2 H-J 2 A0 0.73 5.78 5.78 0.73 a0 4.21 4.83 4.83 4.21 A 0.15 5.78 5.78 0.66 a1 3.27 5.43 5.43 4.60 0.63 1.03 1.03 0.63 T 3 1 0.5 25 0.10 0.05 0.10 0.10 p 100% 95% 86% 90% Table 3.1: The data are taken from Hikspoors & Jaimungal (columns 1 and 4), and Nomikos Soldatos (columns 2 and 3). Given A0, A1, a0, a1, a, and T we calculated c0, and c1 from (3.6). We performed n = 2000 simulations of a standard OU(A, a, a) under IP0 and IP to estimate the power of the test. Using these simulations, we can estimate how much data is needed to obtain reliable parameter estimation for the OU process. We began by looking at some parameter values for OU processes found in the literature. In [44] Hikspoors & Jaimungal study the MROU model (1.9). Having used futures data to estimate the parameters for X, and }‘ in the risk neu tral measure, and also estimate the process i’ itself, they then estimate parameters for }‘ in the real world measure. The column H-J 1 of Table 3.1 gives the values of those parameters. Considering these values, we have that the distribution of H0 and H1 are sig nificantly different taking 3 years of observations, i.e. we distinguish (A0, ao) from (A1, ai) almost perfectly. Nomikos Soldatos [64] consider a regime switching model for the Nord Pool market - see (1.12)-(1.14). The parameters a, i = 0, 1 differ according to whether the weather is wet’ or dry’. We can ask how long a period of observations is necessary to distinguish reliably between a a0 and a a1, in this situation. Simulation results shows that if T = 1 then this is possible at the 95% confidence level but if T = 0.5 then the power of the test is only 47 Chapter 3. MROU model 86% at = 0.1. If we imagine data to be split into fixed periods in which the regime is constant and the task is to distinguish which regime holds during the period, then one can do so reliably if the period is 1 year, but shrinking the period to 6 months will give rise to an error probability of about 10%. In the fourth column we consider two sets of parameters which differ by about 10%: the first set being similar to the values of ao) in H-J 1. In this case we see that even 25 years of data is not enough to reliably distinguish between the two parameter sets. Figure 3.1 plots the distribution of log ZT under the two hypotheses: a substantial overlap in the distribution is apparent. Null hypothesis distribution 0.1 i I I I I I I 0.16 0.14 - 0.12 0.1 - 0.08 - 0.06 0.04 0.02 - —8 —8 —2 2 Alternative hypothesis distribution 0.Z . I I I I I I I I 0.18 0.16 0.14 0.12 0.1 0.06 0.06 0.04 0.02 I I — - 10 —8 —6 4 —2 Figure 3.1: Sampling distributions with (column 4) where we add 10% devi ation from Hikspoors & Jaimungal’s parameters (column 1). ---I —10 I I I I 4 8 8 10—4 48 Chapter 3. MROU model It would be interesting to perform an analysis using the Radon-Nikodym theorem for hypothesis testing with the MROU model. Calculating the like lihood ZT with respect to the u-field u(X3,L3, 0 s T) is straightforward - see [44] (Thm. 3.2). However, we need to calculate E(ZTIXS, 0 s T), and this would require calculating conditional expectations of the form E (exp {JT + ciL5)dW — fT +c1L3)2ds} x8, o s T). This does not seen to be an easy problem. However, from the results for the fourth column, it seems clear that many decades of data would be needed for an accurate estimate of the parameters for the hidden process L. Financial time series have one significant property which makes them unique in a statistical sense. As well as the data itself (e.g. the spot prices), we also have available prices of various derivatives of the product. These produce a substantial amount of extra information. In the context of the IvIROU model, we will therefore use futures as well as spot prices in our parameter estimation. 3.3 Future price This model is an special case of Affine Jump-Diffusion model (AJD), so to obtain closed form formula for the price of futures contracts we use the results by Duffie et al. (2000) [27], see Appendix A. The model can be rewritten as d — 0 dt —A ‘X dt L — LL + 0 ‘L L + — p2x pux dW (3 7)0 dW? Thus, U [Xe, Li]’, and we have dU = (K0 +K1U)dt + dW, 49 Chapter 3. MROU model where dW is a Brownian motion with covariance H0 and — 0 — —)‘ )‘x — PxLK0 — , K1 — , H0 — 2 (3.8)?LL 0 —AL PxL L The functions H1 and lo given in the Appendix A satisfy H1 = 0 and 1 = 0. Duffie et al. gave expressions for various functionals of U, and in particular for ‘(u, t, T, U) := E[exp{u. UT}IUt]. (3.9) By setting u = (1,0)’ in (3.9) one can obtain the future price of U (t, T) := I1((1, 0)’, t, T, U) exp{h(T)} (3.10) = E[exp{XT}Utjexp{h(T)}. (3.11) By equation (A-3) in Appendix A, the future price is of the form E[exp{XT + h(T)}IUj = exp{h(T)} exp{M((1, 0)’, t, T) + N1(1, t, T)X + N2(0, t, T)L}, (3.12) where M((1,0)’,t,T), N1(1,t,T) and N2(0,t,T) satisfy the following equa tions: = AxN1, (3.13) AxN1 + ALN2, (3.14) = LLN2— (uN? + N) —N12pUXJL, (3.15) with the boundary conditions M((1,0)’,T,T) =0, N1(1,T,T) = 1, and N2(0,T,T) =0. 50 Chapter 3. MROU model Solving the initial value problem, we have N1(t,T) = et_T), (3.16) N2(t, T) = e(tT)m — e)t_T)m, (3.17) M(t, T) = (e(t_T) — 1)mi + (et_T) — 1)m2 + (e2_T) — 1)m3 + (et_T) — 1)m4 + (e2t_T) — 1)m5, (3.18) where _______ m2cx rnpcrxcrL m = — rn4 + 7722 Lrn )\X+L )‘X+,\L )\LLm / u om mpuxuLN m2cr m1=— m3=—j—————+ + J 1725=— \4-\X 4x 2.Ax .j 4’\L Thus, we have that the price of the futures contract is given by F(t,T) = exp{M(t,T) +N1(t,T)X +N2(t,T)L}exp{h(T)}, (3.19) where the functions Ni(t,T), N2(t,T), and M(t,T) are as (3.16)-(3.18). The deseasonalized log-future price is given by log F(t, T) F(t, T) exp{—h(T)}. (3.20) 3.4 Formulation in Kalman filter terms We assume that data are available in the form of the spot price S and various futures or forward prices F(t, Ti), i 1,2, ..., m. Most data sets are available in the form of daily prices. Spot prices are traded for all hours of the week, but final2cial markets, which trade future contracts, are only open Monday- Friday. Since in any case, due to low consumption at the weekend (see Figure 5 Appendix B for the calculations) 51 Chapter 3. MROU model 1.4b), spot prices on Saturdays and Sundays exhibit a different behavior to the rest of the week, we will disregard weekends and also holidays, and just consider weekday data. Given this data, we wish to estimate the parameters X, L, X, L, L and p. The Kalman filter method has been applied previously to electricity mod els [14, 51]. To use the Kalman filter we need a discrete time set of equations, so we replace (4.2)-(4.4) with the forward Euler approximation: St e, (3.21) = X_1 — )x(Xt_i — L_1)At + — p2uxAWtl + pcxxAW/, (3.22) L = — AL(Lt_l — L)At + JLAW. (3.23) Here we have made a slight abuse of notation, in writing X,, X_1 for the successive values of X. More precisely, in (3.22) we should write X — —L1)At + — p2uxAW + puxAW, but to avoid two levels of subscripts we have used the form (3.22)-(3.23). Here At 1/250 (the number of trading days in a year), and AW1, AW? are independent Gaussian random variables with mean 0 and variance At. To apply the Kalman filter, the model must be expressed in its state space form. Taking the state variable as u = (Xe, Li)’ a discretization of the time t = t, t2, . . . , and At = (t — t_1) the transition equation becomes: ut=Ct+Dtut_i +qt_i, (3.24) where ( 0 ‘ (1—At AAtCt = LLAt j’ = 0 1— 52 Chapter 3. MROU model and the process noise covariance matrix is / 2A °X’- PxL tZq cov(qt_i) = . (3.25) The measurement equation is given by the functions M, N1, and N2 calculated according to equations (3.16)- (3.18): Vt = ( ) = A+B( ) +r. (3.26) We write X for the observed deseasonalized log-spot price, and Zt,T for the observed deseasonalized log price at time t of a future contract with maturity T > t. We will assume that there is some noise in the measurement of X, and Zt,T, so that X=X+, (3.27) = log F(t, T) + (3.28) Here are iid N(O, o) random variables and are iid N(O, cr) random variables. We have two reasons for making this assumption about non-zero noise. First, at a fundamental level, it is reasonable to allow for some pricing errors due to large bid-ask spreads. (This may be particularly relevant in the futures markets, which are not always heavily traded.) Secondly, the Kalman filter involves matrix inversion - see (2.17), where has to be computed. If the model has degeneracy, then severe numerical problems can arise. Adding the noise in (3.27) and (3.28) avoids this difficulty. Using (3.19) we can write (3.28) as = M(t,T) +XN1(t,T) +LN2(t,I)+. (3.29) The measurement equation is therefore given by xt* ZtT Vt = 1 = A + B L, ) + rt, (3.30) Zt,Trn 53 Chapter 3. MROU model where 0 1 0 = M(t,T1) B = N1(t,T) N2(t,T1) (3.31) M(t,Tm) Ni(t,Tm) N2(t,) and the measurement noise rt has covariance matrix J 0 ... 0 o cj ... 0 = cov(rt) = : o (3.32) o •.. 0 Zm Simulated trajectories of the MROU model using equations (3.24) and (3.26) with parameters ‘x = 130, )\L = 3, ux = 5, = 0.5, L = 3.5, p = 0.3 and Zt 1/250 can be seen in Figure 3.2. The long-term mean process L reverts towards the mean L, and as we expected, the spot price St and the future price F(t, T1) mimic the long-term mean process but with different volatility. The observation and state equation matrices C,, D, q, A, B, and depend on the unknown parameters of the model. Based on this state-space formulation the parameters that we need to estimate are: & = {x,L,ux,JL,L,p,Js,uz}. Note that if we use two different maturity contracts then rn 2, and we will have and crz2 in the parameter vector 6. The log-likelihood function Jü(Vit) for the linear Gaussian space-state is given by equation (2.43). This function can be maximized with respect to 0 using an appropriate numerical optimization procedure. O(Vi.t) only depends on the prediction error Ct and its covariance matrix F1_. Both in turn are outputs of the Kalman filter, equations (2.14) and (2.16). Thus the maximum likelihood estimate of & can be obtained as follows: 54 Chapter 3. MROU model 0 Spot S and Future F(t,T1)processes Figure 3.2: The upper graph shows the simulated spot price S and the future price F(t, T1) with maturity of one month. The lower graph is the long-term mean process L. Algorithm 4 (Kalman Filter optimization) • Step 1 Choose a initial value for 6, say 6o. • Step 2 Run the Kalman Filter (Algorithm 1) and use the sequences et and F1_ to compute the log-likelihood £(Vit) by (2.43). • Step 3 Employ an optimization procedure that repeats Steps 1-2 until a max imizer 6 of (2.43) has been found. Days Long—term process L Days 55 Chapter 3. MROU model Some practical problems arise with the optimization procedure, and so the performance and accuracy of the Kalman Filter are affected since the problem may be poorly scaled. An optimization problem is poorly scaled if changes in the decision variables produce large changes in the objective function for some components and not for others [25]. We solved this problem by rescaling the variables. There are several numerical search algorithms available to maximize the log-likelihood (Step 3). Barlow et al. [4] used the Nelder-Mead method6 to minimize £Q(Vi.t). Here we decided to apply a quasi-Newton algorithm, the so-called BFGS method7 (which is similar to the method used by [51]) to get the initial point, and then used the Nelder-Mead method to find the optimum. 3.5 Empirical results In this section we report some empirical results based on simulated and real data to examine the Kalman filter method applied to the model described in the previous section. 3.5.1 Simulated data We first ran the algorithm on deseasonalized simulated data. To analyze the algorithm performance, first we simulated series with 100 and 800 ob servations respectively using equations (3.24) and (3.26), considering only the nearest monthly futures contracts in the log-future price, i.e. Zt,T1 with T1 = 30 days. In this case we have a 8 dimensional parameter space. We started the maximization procedure with a different initial values each time. Examples of some runs are given in Table 3.2. Since we are searching for the maximum likelihood, we did 25 runs and we took the one which gave the largest value for £(8). Repeating the same 6The Nelder-Mead method is used to minimize a function of multiple variables without derivatives [see Nelder and Mead (1965)]. 7Details and derivation of the BroydenFletcherGoldfarbShanno method in the context of filtering can be found in [28]. 56 Chapter 3. MROU model Ti = 300 True Value Run 1 Run 2 Run 3 Run 4 Run 5 Ax 130 129.032 127.923 126.335 138.586 130.232 AL 3 3.478 3.301 2,348 2.788 2.985 5 4.917 4.845 4.658 4.795 5.123 0.5 0.522 0.489 0.453 0.522 0.509 L 3.5 3.495 3.509 3.503 3.499 3.497 p 0.3 0.318 0.313 0.247 0.382 0.311 0 0.004 0 0.005 0.000 0 crz1 0 0.009 0.002 0 0.000 0 -620.952 -629.702 -643.223 -631.359 -601.72 CPU time 163.141 192.954 152.183 155.216 123.2167 Table 3.2: Five different maximization runs, on the same set of simulated data. n=1000 m=300 True value Estimator Std. Estimator Std. Ax 130 129.369 1.055 129.042 4.221 AL 3 2.954 0.237 2.968 0.643 x 5 4.998 0.098 4.862 0.127 L 0.5 0.493 0.004 0.489 0.023 L 3.5 3.507 0.001 3.502 0.009 p 0.3 0.308 0.018 0.301 0.077 °s 0 0.0003 0.0007 0.001 0.002 z1 0 0 0.003 0 0.001 Table 3.3: Estimation using one futures contract (average of 50 simulations). procedure with 50 different series we obtained the following results, which are given in Table 3.3. We can see that the estimation results recover the true values reasonably well in both cases. As expected, the standard deviation increases for all estimators when we reduce the number of days. Note that the estimator is close to zero for variables s and z1• This is not surprising since we are estimating the true model that generated the data and the noise in the model only comes from two sources. Next, we simulated 25 new series with 300 data, but this time we included 57 Chapter 3. MROU model log-futures prices with two different maturities, T1 = 30 and T2 = 60. In this case we have an extra parameter to estimate, az2. In Table 3.4 we summarize the results, repeating the same procedure as before. TL 300 True value Estimator Std. 130 129.801 0.436 3 3.000 0.000 ax 5 4.845 0.178 aL 0.5 0.499 0.001 L 3.5 3.505 0.017 p 0.3 0.294 0.102 as 0.010 0.009 az1 0.001 0 az2 0 0 Table 3.4: Estimation using two futures contracts (n=300). Thus we obtained good approximations for all the parameters using only a short data series. This could be useful due to the scarcity of data in the electricity markets. One example is [16] where the authors comment that there is too little data for parameter estimation in the UK market. We remark that if we tried to estimate the model parameters just using the spot price, then it will require many decades of data to make accurate estimates. Notice that we used an alternative Kalman filter formulation (U-D fil tering) since numerical problems arise when we include more parameters to estimated. During one of the recursions of the filter the covariance matrix failed to be positive semi-definite and consequently the estimated parameter differed from the true values. This problem arose because the matrices . and Zq were ill-conditioned. Since we are minimizing a fixed function, it is legitimate to discard runs which fail in this fashion. For further details about the U-D filter refer to [15, 38, 75]. 3.5.2 The German electricity market We now wish to apply these techniques to real data. While there are many markets in which electricity is traded, the data we require (that is, both spot and futures prices) in many cases are not available. For example, the Alberta 58 Chapter 3. MROU model Power Pool makes spot prices available, but forward prices are known only to market participants. One market for which both sets of data are available is the German EEX market. The European Energy Exchange (EEX) is Germany’s energy exchange. It is one of the biggest power markets in Europe. EEX emerged in 2002 from the merger of EEX Leipzig Power Exchange and EEX European En ergy Exchange Frankfurt. Both exchanges initially started spot trading for physical contracts in 2000. In 2001 EEX Frankfurt also introduced trading of standardized financial contracts. Commonly traded products in the power markets are baseload, peakload and hourly contracts. At the German market the times for peakload are defined as weekdays between 8:00am and 8:00pm. In the futures market contracts on both baseload and peakload are traded. The usual delivery periods are one month, one quarter and one year. In the spot market of the EEX baseload, peakload and hourly contracts up to the next weekday are traded. The estimates are based on historical daily average spot price and monthly baseload futures price covering the period from July 2002, when the EEX and LPX markets merged, until the end of June 2007, almost five years of historical data. This data contains prices for 1267 days. Figure 3.3 depicts the price trajectories of the spot and the nearby monthly futures prices between July 1, 2002 and June 29, 2007 for the EEX market. From the graph we note that there is a strong mean reversion and the spot prices show extreme spikes as well as high volatility which changes rapidly over short time periods. Moreover there is an linear drift over the years in the spot and futures prices. Figure 3.3 shows much greater volatility, and more price spikes, for the period Jan 2005-June 2007 than in the earlier period, Jul 2002-Dec 2004. We therefore split the data into two parts (7/1/02-12/31/04, 1/1/05-6/29/07), and repeated the runs on each of these. Following [8], we removed the seasonality by representing it as a linear combination of cosines including a trend, a weekly, and a yearly cycle of the form (r1+2i t’\ 1r2+2i t\h(t)=r+/3ot+/3cosç 250 ) +j32cos ). (3.33) 59 350.00 300.00 250.00 200.00 150.00 100.00 50.00 0.00 Chapter 3. MROU model Figure 3.3: market. Electricity spot and nearby monthly futures price in German Assuming 250 trading days in a year, we estimate the seasonality by fitting the h(t) function to the log-price series by ordinary least squares. The results can be seen in Table 3.5 and Figures 3.4, and 3.5. We see less strong season ality in the EEX market that in the more hydro-dependent Nordic electricity market that depends more on hydro. The seasonal component is highest in winter. Parameter ñ i Est. value 3.2965 0.0005 -0.0695 23.5235 0.0224 0.8218 Table 3.5: Estimated values for f(t) by least-squares fitting. Based on equations (3.24) and (3.26) we estimated the model using the spot price and one futures contract with one month maturity, for both peri ods. Also we estimated using two futures contracts with one and two month —Fl —SpoL Price 7/1/2002 1/1/2003 7/1/2003 1/1/2004 7/1/2004 1/1/2005 7/1/2005 1/1/2006 7/1/2006 1/1/2007 60 Chapter 3. MROU model .c) 5.5 4.5 3.5 2.5 Log—spot price Iog(S) = X + h(t) Figure 3.4: The upper graph shows the log spot-price of the EEX market and the seasonal component h(t) and the lower graph the deseasonal series = logSt — h(t). maturities. The results of the parameter estimation ale shown in Tables 3.6 and 3.7. For all the periods the estimates of p are quite small - the largest value being 0.115 for Part 1 when estimated using one future contract. This may be compared with the estimate —0.96 obtained in [44] in the context of the oil market. Unlike the case of simulated data, the parameter estimates using two futures prices differ somewhat from those obtained with just one future. The explanation is presumably that the model does not perfectly describe the real data. Iog(S,) h(t) fitted 200 400 600 800 1000 1200 Days Deseasonalized log—spot price X = log (Se) — h(t) 600 Days 61 Chapter 3. MROU model Spot price S and exp(h(t)) Figure 3.5: The upper shows the spot price S with exp{I(t)} and and the lower graph the deseasonal spot price St = exp{Xt}. In the absence of noise, the model (4.2)-(4.4) gives exponential decay toward the long run mean L. We can therefore interpret as the ‘half lives’ of the processes X and L (measured in years). For the X process the estimates in Tables 3.6 and 3.7 give a half life 1.5 — 4 days, while for the L process the half life estimates vary from about 6—12 months. Thus the estimated process X and L do play a satisfactory role in separating out short and long-term fluctuation in the spot price. The long run standard deviation of the OrnsteinUhlenbeck process L is (aL/2,\L)1/2.Using the numerical values for the whole period from Table 3.7, q) Days Deseasonalized spot price S q) 600 Days log 2 TX = Ax log 2 TL = AL 62 Chapter 3. MROU model Whole Part 1 Part 2 Estimator Estimator Estimator Ax 65.030 42.439 77.509 AL 0.544 1.319 1.060 1.719 1.172 2.013 0.517 0.410 0.586 L 3.108 3.105 3.156 p 0.101 0.115 0.065 0.177 0.193 0.167 Q.z1 0.000 0.000 0.000 CPU time 613.88 485.99 312.33 Table 3.6: Estimated values for the EEX market using S, and F(t, T1). Whole Part 1 Part 2 Estimator Estimator Estimator Ax 98.507 116.807 98.227 AL 0.441 1.231 0.579 3.161 3.191 2.965 cxL 0.537 0.337 0.507 L 3.169 3.168 3.222 p 0.021 0.019 0.114 us 0.000 0.001 0.047 0.089 0.074 0.065 __________ 0.000 0.006 0.024 CPU time 1440.58 1208.56 703.68 Table 3.7: Estimated values for the EEX market using S, F(t, T1) and F(t, T2). this gives an standard deviation of 0.57. The corresponding quantity for the X, process is 0.23, so both components contribute significantly to the long run variance of the log-spot price. In order to investigate if the estimated parameters make sense, we sim ulated a path of the spot price, future prices and long-term mean process that describe the model using the estimated values from Table 3.6 Part 1, see Figure 3.6. We now compare those simulations with Figure 3.5, which gives the real 63 Chapter 3. MROU model Spot S and Future F(t,T1)processes 3.u I I I I I L, L 3. 0 100 200 300 400 500 600 700 Figure 3.6: Simulation of spot and future prices (upper graph) and long-term mean process (lower graph) using estimate values for part 1. deseasonalized data. The generated trajectories differ somewhat from those observed in the EEX market. The most notable difference is the absence of price spikes: the simulated data is all in the range €10 - €401 while the real deseasonalized data has about 6 spikes with prices above €70. The absence of such spikes is not surprising since the model contains no mechanism for generating them. The second feature is that, as with the real data, the spot and futures price do tend to follow each other. This is not surprising, since both prices do relate to the same commodity. Since the long run process L is available for the simulated data, we have also plotted this - see lower graph in Figure 3.6. Comparing this with the future price in upper graph we noted that for those parameters, the dominant effect on the future price is from the oscillations a •0 Days Long—term process L 3.6 3.4 28 2.6 Days 64 Chapter 3. MROU model of the long-term mean. The graphs in Figure 3.6 and Table 3.7 together show that while this model does capture some features of the real data, it does have significant defects, in that it does a poor job of capturing the extreme events, spikes or jumps, which appear in the real market. Another test for the appropriateness of the model is to compare empirical moments for the real data with those from simulated data. If we compare the first four empirical central moments of the log-return8 of the sequence of the simulated and real spot prices we can see that there is a good fit for the mean value and for the standard deviation, see Table 3.8. The empirical distribution has fatter tails than the normal distribution (kurtosis > 3), indicating a higher occurrence of extreme events, i.e. jumps. Real data (1) Sim. data (1) Real data (2) Sim. data (2) Mean -0.0019 -0.0002 0.0000 0.0002 Std Dev. 0.2214 0.1770 0.2066 0.20402 Skewness -1.3125 -0.0037 -0.0800 -0.0160 Kurtosis 29.8250 2.9855 8.1288 2.9313 n data 634 634 631 631 Table 3.8: The table shows the first four moments of the logarithmic de seasonalized price returns of observed data and the average of 50 simulated trajectories. 8Log-return for a sequence of prices S are defined as i 1n(Sj+j./S) 65 Chapter 4 MROU with jumps In this chapter we consider a jump-diffusion model. This model is similar to the ones considered by Hikspoors and Jaimungal [44] and Nomikos and Soldatos [64]. [44] considers a model of the form = exp{h(t) + Xt + J}, (4.1) where X is the first component of a pair (Xe, Y) satisfying (1.9), and J is an independent jump process, see equation (1.11). In one respect our model represents a simplification of the model in [44], in that there is only one Gaussian factor. However, unlike the model given by (4.1), the Gaussian and jump component in our model are not easily separated into independent processes. There are several reasons for considering a jump-diffusion model. First, actual spot prices do exhibit spikes - see Figure 3.3, and adding jumps to the process is one way of modeling this. Second, data for spot prices show that the fourth central moment (kurtosis) of the log-returns is much bigger than 3 (see Table 3.8). Diffusion models such as the MROU model in Chapter 3 tend to give the kurtosis close to 3. 4.1 Description of the model Let (Q, F, Q) be a complete filtered probability space. The dynam ics of the state variables are given by the following stochastic differential equations: = (4.2) dX = —)x(X — L) dt + x dW’ + JdN — (4.3) dL = —AL(Lt — L) dt + L dW2. (4.4) 66 Chapter 4. MROU with jumps As before S is the spot price. Here the Brownian motions W’ and W? are independent. The jump behavior of X is governed by two types of jumps: upward jumps and downward jumps. The upward jumps J are exponentially distributed with positive mean i/i, and the downward jumps J are also exponentially distributed with mean 1/nd. In this model, N and N are two independent Poisson processes with arrival rates A and respectively. The function h(t) denotes a deterministic seasonality function. 4.2 Valuation of electricity futures This model also belongs to the class of (AJD) process. We can rewrite equations (4.3) and (4.4) according to equation (A-i) in Appendix A to obtain: d — dt ‘X ‘X X dt x 0 dW L — LL + 0 L L + 0 L dW? JdN1 — J?dN 0 Again defining U = (Xi, Li)’, the expressions for K0, K1, H0, and, H1 remain the same as for the MROU, see equation (3.8). However we now have 10 = >u + ‘d. (4.5) To obtain the formula for the future contracts F(t, T) we need to calculate the jump transform function in order to include them in the ODE for M(t, T). The other two equations N1(t, T) and N2(t, T) are the same: see equations (3.16) and (3.17). The density of the distribution of jumps of X is given by 67 Chapter 4. MROU with jumps x>O vX(X) = d x < 0 The transform for the jumps is given by (A-4), = f eN1z1+N2z2d(zi,z2), (4.6) where (•) is the jump distribution on R2, and N1, N2 are such that the integral (4.6) converges. However, since L does not have jumps, is con centrated on the subspace z2 0 and we have (N1,2)= feN1zd(z), (4.7) wherever this integral converges, that is, wherever —‘rid < N1 < . Then = f:exp{Nlz}v(z)dz 00 0 = u f e_(_N1)zdz + d f de1)dZ0 u+d - — ( 7u d ( 7?d — u+Adu—Nl) +Add+N1 Therefore = -ALLN2(JN+u ) A ( —1—A ( dU U—Nl(tT) d d+Nl(tT) )‘ with the boundary condition M((1, 0)’, T, T) = 0. 68 Chapter 4. MROU with jumps The equations for N(t, T) are, as before, \xNi, and = —\XN1+; L2. We can solve the system with the corresponding boundary conditions to obtain M(t, T) = mi(et_T) — 1) +m2(e(t_T) — 1) +m3(e2t_T) — 1) + — 1) +m5(e2t_T) — 1) (4.8) + in (?7u_— e__(t_T1 + in (‘1d + eX(t_T)’ ?]—1 / \x ??d+l / where the constants m1,m2, ..., m5 and the solution for N1(t, T) and N2(t, T) are given by equations (3.16) and (3.17). See Appendix B for the calculation of M(t, T). Thus, the expression for the future prices is given by E(t,T) = exp{h(T)}exp{M(t,T) +N1(t,T)X +N2(t,T)L}, (4.9) where the functions M(t, T), N1 (t, T) and N2(t, T) are given by equations (3.16), (3.17), and (4.8). Note that (4.8) requires that < 1. This restriction is to be expected, since the future price is given by (t,T) = E(STI) = E(e)TT)IT). (4.10) Since XT contains, in general, terms involving jumps with an exponential distribution, the expectation in (4.10) will diverge if the upper tail of the jump distribution is sufficiently large. See [44]. 69 Chapter 4. MROU with jumps 4.3 Particle filter setup Since this model is non-Gaussian, due to the jump process JdN — JdN, we cannot employ the IKalman filter to estimate the parameters )‘, ‘d, 71d We therefore wish to employ the particle filter, which is in principle able to handle quite general distributions. As for the MROU model, we assume the data is given by X, Z1, t= 1,2,...,n. (For simplicity we just considered one future price). Here X7 is the deseason alized log-spot price, and Zt,T1 the deseasonalized log-future price F(t, T1). The first step is discretise the model. Let to, t1, ..., t, be the times, and = t, — t_1. As before for simplicity we abuse notation slightly, and write X for X1. Then we have, using the forward Euler approximation: Xt = — x(Xt_i — L_i)L\t + uxAW + JN — jdNd = L_1 — L(Lt_l — L)t + JLW. Here J are exponential random variables with parameters , ‘rid respec tively and N(O,tu), N(O,zto). zNu, /N1 are Poisson random variables with parameters )At, )\dAt re spectively. Since P(/N 2) = 1 —et— /te_>t the probability of two or more jumps in one day is small. We therefore approximate AN by a Bernoulli random variable with parameter that is we take P(AN = 1) = P(z\N = 0) = 1 — 70 chapter 4. MROU with jumps with similar approximation for AN. The transition equation is therefore (Xt N 1 0 N (1— AxAt AxAt NI X_1 Ut= L ) LLAt 0 1—ALAt ) L_1 411 + As before, we assume the measurements are subject to noise, so the mea surement equation is Ix*N / 0 N / 1 0 N/X Vt = M(t,T) ) + N(t,T) N2(t,T)) L (us 0 NI 412) + ü (. where °, , and are standard normal random variables. For the particle filter we initially assume the parameter 6 is known, and set it up by writing in the following way. We write i4(k) for the kth component (k = 1,2) of the ith particle (i = 1,2, ..., ni,) at time t. A key part of the implementation of the particle filter is the choice of the proposal density q(utut_i,vt). One choice would be to simply take q() to be the transition density p(utIuti) arising from (4.11). However, this choice is not likely to be optimal. Most of the time, we will have AN AN = 0, so most particles will not make a jump at time t. If the observed process X, does make a jump at time t, then most particles will be left behind by this jump. It is therefore better to choose q(ututi, Vi:t) to exploit the information available in Vt. (We remark that the optimal choice of q(ututi, Vit) would be the distribution p(utIuti, Vit), but it is not feasible to calculate this distribution). We therefore choose q(.) so that the particles are propagated according to the following equation: 71 Chapter 4. MROU with jumps u(1) X + u3t, (4.13) u(2) = —(u1(1) — L)At + JL/tt. (4.14) Note this means that the second component of u just follows the evo lution equation of L, but we use the new data available in X to move the first coordinate of the particles close to X. Here are independent N(O, 1) random variables (i = 1, 2). We now calculate the likelihood and prior densities. We have p(vtU) = p(X, Zt,TIU(1), r4(2)) = p(XIU(1))p(Zt,TIU(2)), (4.15) where xIu(1) N(U(1),u) and Zt,TIU(2) N(mz,u) (4.16) with mz = M(t, T) + N1(t, T)U(1) + N1(t, T)U(2). Now, for the prior density or transition density, (i) (i) (i) (i) (i) (i) it Iu_1 = P(Ut (1), Ut 2)Iu_1(1),u_12)) (4.17) = Here U(2)IU1( ) N(mL,SL) (4.18) with mL = L(1 — e_t) +U1(2)e_t and SL = (u/2L)(1 — e_2t), since L is an Ornstein-Uhlenbeck process. 72 Ohapter 4. MROU with jumps Using the Bernoulli approximation for and AN, and neglecting the O((At)2)probability of both and upward and downward jump in the same period, p(u(1)Iu1(1)) = (1- )At - dAt)fO(U(1)IU2l(1)) + + AUtfd(ut(1)Iutl(1)). (4.19) Here fo denotes the density corresponding to the jump-free case and f and fd correspond to the case of a single upward and downward jump respectively. Set /x u1(1) — (11) — Lx)t, and s2 = then fo(u(1)Iui(1)) = (2ns2)h/2exp{(_(u(1) - x)2/2s2)}. (4.20) f and fd are obtained by the convolution of fo with the distributions of JU and jd So f(u (1) u1(1)) = (4.21) Now 73 Chapter 4. MROU with jumps (i)(i)(u (1) — y — x)2/2s + uY = (2s)’ (( — (U (1) — x))2 + 2s2y) = y — (u (1) — x))2 + 2s2y)(2s)’ (( (t) (i) = (2b)’ (( — (u (1) — — s2))2 —(s4 — 2(u(1) — x)s2)) / (%) — (y — (1) — — — 2s 2 2 (i) —s T/ + (u (1) — x)u. (i) 2Hence writing A(x) = (u (1) — — s (i) (i)f(u (1) Iu_1(1)) = ue322e_1)_ P00 I0 P00 J (2)_1/2et/ dt A(x)/s = e82/2e_+(1)_) I (2)_h/2e_t/ dt -00 = where J(.) is the normal cumulative distribution function. So (i) (i)f(u (1)Iut_l)) = ue822e_1)_ (4.22) — — (z)Let u” = —Ut (1). To get the distribution fd we use the calculations for f, 74 Chapter 4. MROU with jumps (z) fd(ut (1)I21(1)) = f(2ns2)_h/2e(_1)_2/2s2dedy 0 T= (27rs)_1/2e(_(u’_Y_Px)/ s)e_??d dy0 = So (i) (i) fd(ut (l)u_i(l)) = (4.23) — —s2d)/(ux(At)’12)). (i)(i)Combining (4.19), (4.20), (4.22), and (4.23) gives p(u(l)Iu_, ). Finally, the proposal density is z) (i)q(u Iu_,,v) = p(u(1),u(2)Iu,(1),u(2),v) = p(u(1)Iu(1),vt) p(u(2)Iu,(2)), (4.24) (i) with P(t (2)Iu1)) given as (4.18) and p(u(1)Iu,(1), Vt) = (2nu1l2 exp{(-(u(1) -X)2/2u)}. (4.25) Combining (4.15), (4.17), and (4.24) we calculate the associated weights (i) (i) (i)(i) p(VtIu ) P(’Ut u._,) wt (i) (i)q(ut Iu_i,Vt) (i) (i) = (2)u,(2)) p(u(1)Iu1(1),Vt) p(u(2)u,(2)) (i) = p(Zt,Ttu(2)) P(Ut (1)1(1)). 75 Chapter 4. MROU with jumps Here the first term is given by (4.16), and the second by combining (4.19), (4.20), (4.22), and (4.23). 4.4 Simulated data with known parameters We tested the implementation of the particle filter presented in Section (4.3) for the MROU jump model to estimate the long-term mean process L. We simulated a series of 100 time data points according to (4.3), (4.4), and (4.9), taking and = 1. Figure 4.1 shows the particle filter estimates of the state using 100 particles. Figure 4.1: Plot of the true state L and estimate of the particle filter. As is clear from Figure 4.1 the particle filter does a good job of estimating the transition state for this example even using just a few particles, provided the parameters are known. 76 Chapter 4. MROU with jumps 4.5 Likelihood function estimation Using the particle filter for a fixed parameter 0, one can obtain an estimate of the likelihood function Jü(Vit) by (2.45). However, as we already mentioned above, severe difficulties arise when one tries to optimize this function since, for each value of 0 one is using a different randomization. To investigate the problem, we fixed all the parameters except A at their correct value (Ax = 110), and estimated £x(Vit) for simulated data. Taking 1000 points and 2000 particles, we estimated £ (Vi:t) for Ax in integer steps between 102 and 120. As is clear from Figure 4.2, there is little prospect of satisfactory use of a hill-climbing algorithm given the level of noise. We repeated this estimation with 15000 particles and obtained a similar curve, but with oscillations roughly 10 times smaller. These oscillations however were still large enough so that a hill-climbing algorithm would have difficulties in locating the optimum. It seems likely that increasing the number of particles by a factor of k will increase accuracy of the estimates of the log-likelihood by If so, one Figure 4.2: The log-likelihood for different Ax values. 77 Chapter 4. MROU with jumps might need 1O61O8 particles in order to use a hill-climbing approach even in one dimension. 4.6 Sequential parameters Our goal is to estimate the long-term mean process together with the 11- dimensional parameter vector 0 given by according to the information available at a given time t, i.e. Vit. As we explained in Section 2.6, from a Bayesian point of view we concatenate the state vector and the parameters and apply the filter to this augmented state. Then we define Yt (Ut) (4.26) where Ot is the particle approximation to 0. It is useful to add some noise to the transition for 0 19t Ot—1 + ct_i (4.27) where {t—1}t1 is a small artificial noise with a decreasing variance A with t. A wide variety of choices of and A are possible. It is clear that some care has to be taken in the choice of A. If A is too large, then particles & will oscillate too much to give a satisfactory estimate of 0 (the larger the covariance, the more quickly older data are discarded). If A is too small, then unless the initial value is nearly correct, 0 will not move enough to reach the correct value. The literature contains a number of suggestions on how A should be chosen. In [43] the author suggested to be white noise. Others proposed 91f not, we oniy sample particles in 0 space at time 1 and never modify their locations, then after a few time steps p(9IVit) is approximated by a single particle. 78 Chapter 4. MROU with jumps to set A as a diagonal matrix annealing to zero - see [41, 53, 68] for more options. First we take a fairly simple approach. We will take A to be diagonal, and write t9, for the j-th component of &, and for the j-th element in the diagonal of A. We set t,j =b (4.28) where b and cj are constants. Note that var(t,j) = 00. The ‘annealing’ given by (4.28) is not fast enough to make Oj converge almost surely. Later we use the Liu and West [57] approach. They suggest approximat ing the distribution p(OIvit) by using a mixture of Gaussian distributions, that is p(OIV1:t) n8m,h2At. (4.29) The quantity m = + (1 — a) is the kernel location for the i—th component of the mixture where = (4.30) and the matrix A is an estimate of the posterior variance covariance-matrix, i.e. = — )(& — t)’. (4.31) The constants h and a, that measure the extent of the shrinkage and the degree of over-dispersion of the mixture, are given by h2 = 1— ((2y— 1)/2)2, a = — h2, where the discount factor ranges between 0.95 - 0.99. In this case the artificial dynamic of the parameter is given by (i) 1/2 = m + hA Et+1, ‘ n(0,I), (4.32) 79 Chapter 4. MROU with jumps where n(xIp u2) denotes the density of the Gaussian distribution. With these specifications the likelihood, prior and proposal densities for z4 now depend on 8. So we have p(vtIy) = n(XIu(1),u)x n(Zt,TIMt(6) + Ni,)u(1) +N2t(O)’u(2), ui), (4.33) (i) (i) (1) (i) (i) (i) P(Yt I-) =P(ut (1)Jut_i(1),ut_i(2),Ot_i) x p(u(2)Ju1(2 ,6) fln(IO, , (4.34) (i) (i) (i)(1) * 2 (i) (i) (i) q(yt y_1,vt) = n(u IX ,us)p(Ut (2)I_1, x fln(OIO). (4.35) Substituting (4.33)-(4.35) into the weights formula gives (i) — p(vt iu) p(uJu1) ‘wt — (i) (i)q(u u_1,vt) = n(Zt,TIMt(O) + Ni,t(O)u (1) +N2,t(O)u (2), x p(u(1)Iu21(1 ,u2(2), O2). (4.36) The last density in (4.36) is given by (4.17). 4.7 Empirical results In this section we report some empirical results based on simulated data to show the performance of the method. The true parameters were fixed to be: 80 Chapter 4. MROU with jumps )x=l’O x3 \i=s L=3.2 L’ 1?1.5 A=5 7d=3 d=’ For simplicity we took us = 0.1 and z = 0.1 to be fixed and assumed to be known. Figures 4.3 and 4.4 show a simulation of a MROU with jumps process using the parameter values as above.The upper graph in Figures 4.3 shows the simulated spot price and the future price with maturity of one month. The lower graph is the long-term mean process L. Figure 4.4 shows the upward and downward jump processes. For this data set the skewness and kurtosis is -1.8895, and 21.069 respectively. Considering the first approach (4.28) and using 1000 particles we ran the particle filter to optimize all the parameters, but the filter failed to obtain good results on most of the parameters. Difficulties of this kind in a high-dimensional situation are not surprising. Even in the case of the Kalman filter, where we optimized a deterministic function, the hill-climbing algorithm sometimes failed to find a point close to the true maximum. The particle filter replaces the hill-climbing point with a cloud of particles, where the weights w are higher when the particles are close to the true value. This should mean that the particle cloud drifts towards the true value of the parameters, but it seems intuitively unlikely that it will perform as well as a hill-climbing algorithm. A known weakness of optimization algorithms is the following. The higher the number of parameters, the worse the performance of the algorithm. This means that a one-parameter optimization should perform best. To test this, we allowed in turn each of the parameters to vary. Thus in each run we fixed all but one parameter at its correct value, and ran the particle filter with just one 0t,j varying. We took the number of particles to be 1000, and the number of time steps to be 1200. We chose as in (4.28), with cj = 100 The Figure 4.5 shows the dynamic under the artificial noise (4.28) of the parameter := for t = 1, 2, ..., 1200 using 1000 particles and its optimal value (110). To test the particle filter, we started the parameter some way away from its correct value. The noise is decreasing quite slowly, 81 Chapter 4. MROU with jumps 450 400 350 300 250 200 150 100 50 SpotS I dXt=—?.x(Xt_Lt)dt+axdWl +JdNJddNd Figure 4.3: The graph shows a simulated trajectory of the MROU model with jumps. and hence the variance of the process remains large. Figure 4.5 suggests that there is little “push” in the particle filter towards the correct value, and this is confirmed by Table 4.1, which shows a fairly wide dispersion of estimates around the true value, for different runs of the filter on the same data set. A similar pattern arises for the other parameters. We ran the particle filter 25 times for each parameter. Let O(k) be the value of 8t,j in the kth run, where k = 1, 2, ..., 25. For each run k, we estimated by taking its average over the last 400 observations, that is, Sexp(X1)’ F(t,Ti)=Et(exp(Xri)) 200 400 Long—term mean Lt I dLt = — 2L ( Lt — L ) dt + cYL dW2 600 800 1000 12 Days Days 82 I I I I 0 200 400 600 800 1000 1200 Days Figure 4.4: The graph shows the upward and downward jumps for the sim ulated trajectory of the MROU model with jumps. 1200 = ZOt,j(k). We then calculated the mean and standard deviation of (k), k = 1, 2, ..., 25. Table 4.2 summarizes the outcome of the whole procedure. The first column is the true value of the parameters. The second column is the initial value of , that is b in (4.31). The third and fourth columns give the mean and standard deviation of the estimates (O(k), k = 1, 2, ..., 25). As we can seen from the Table 4.1, the performance results are rather Chapter 4. MROU with jumps Upward juwps Ju a Jd=ExponentiaI0d) 600 Days Downward jowpn Jd 025 0.2 0.15 0) 0.1 0.05 83 Chapter 4. MROU with jumps Figure 4.5: Sample particle filter trajectories for the estimate of Ax. True ri r2 r3 r4 r5 r6 r7 r8 Ax 110 126.13 161.49 102.76 87.72 76.41 105.29 124.15 83.71 Table 4.1: Sample of 8 estimated values for Ax. True value b mean value std. Ax 110 30 99.29 19.64 x 3 1 3.26 0.048 AL 5 2 8.147 0.729 L 1 1 1.16 0.017 L 3.2 1 3.37 0.037 , 1.5 0.5 2.444 0.754 A 5 0.5 8.301 2.983 rid 3 1 14.87 6.917 Ad 1 0.5 4.218 1.817 Table 4.2: Individual estimates for parameters in MROU with jumps model. limo step 84 Chapter 4. MROU with jumps mixed. In some cases the particle filter was able to obtain a reasonable estimate for the parameter, but in others, particularly for the parameters associated with the jumps, the estimates are far from the true value. The particle filter literature suggests that even with a very large number of par ticles one may not be able to obtain accurate results. Next we tested the Liu and West approach. We used 4000 particles to approximate the distribution of interest. A problem that we noticed is that the estimated posterior variance-covariance matrix A collapses to zero after a few hundred iterations. We solved this problem choosing an efficient resampling scheme that kept low the variance in the particle filter algorithm. We found that the residual resampling kept the covariance matrix positive. We ran two examples. In the first one we fixed all but one parameter at its correct value, and ran the algorithm choosing a reasonable initial point for the free parameters. We obtained good results even for the jump parameters. (For simplicity we consider again a5 = 0,1 and z = 0.1). The results are displayed in Figure 4.6. It is interesting to note that the algorithm detects precisely the parameters Xi,, A, r, 71d associated with extreme events (spikes). For the second example, we took the vector parameter to estimate as 0 = {Ax,AL,ax,uL,L} while the rest of the parameters were fixed to their optimal values. Again, we started the algorithm choosing a reasonable initial point for the free pa rameters 0 as well a small initial variance. Figure 4.7 presents part of the results. From this example, we noticed that the algorithm implemented gave ac curate estimates for the parameters L, crL, and More difficulties arose when estimating the speed reversion for long and short-term process Ax, and AL. There are slightly under and over estimated respectively. Overall, the algorithm provided more precise estimated values but diffi culties arose for the parameters related to the jump process. Further, the algorithm required a significant amount of tuning, i.e. choosing the initial 85 Chapter 4. MROU with jumps value and variance of the artificial noise. However the graphs in Figures 4.6 - 4.7 indicate a significantly better performance than that obtained by using (4.28). Given these difficulties, we did not feel confident that the algorithm would perform satisfactorily on real data. 86 Chapter 4. MROU with jumps Figure 4.6: Sample particle filter trajectories for the estimate of , and T)d using Liu and West approach. 2. 100 200 000 4 000 600 700 000 900 tine step It It One step 87 Chapter 4. MROU with jumps 100 200 0 400 505 600 000 500 900 1000 Figure 4.7: Sample particle filter trajectories for the estimate of ux and L using Liu and West approach. 88 Chapter 5 NLMROU model Our difficulties with the particle filter led us to look for more tractable models which have the potential for explaining price spikes. While the incorporation of jumps is the most natural way to account for price spikes in the spot price, other explanations have been offered. Barlow [3] proposed a non-linear diffusion model, which can produce price spikes similar to those observed in real data. This single factor model is unlikely to provide a good explanation for the observed relation between spot and future prices. Here we present a two-factor model of the same kind. The model is estimated using data from the European Energy Exchange. 5.1 The model The nonaffine term structure two-factor model for futures prices is known as the Non-Linear Mean-Reversion Ornstein-Uhlenbeck model (NLMROU). It uses the inverse of the Box-Cox transformation to generate price spikes that fit the observed data observed in the power market. Let (Q, F, {}t>o, Q) be a filtered probability space. The dynamic of the spot price under the risk-neutral measure Q is the following: = f(X)h(t) where X is an Ornstein-Uhlenbeck process which reverts to a stochastic mean L, also fluctuating as Ornstein-Uhlenbeck processes and fa is the inverse of the Box-Cox transformation. This transformation was introduced in the context of electricity markets in [3]. The deterministic component h(t) incorporates the seasonality effects in the model. More precisely, 89 Chapter 5. NLMROU model I (1+aX)’/, 1+aXt > 6oS = h(t) x 1/a 1 (5.1) 1+aA.tEo dX = —Ax(X — L) dt + x dW, (5.2) dL = \L(Lt — L) dt + 01 dW. (5.3) where the two Brownian motions W’, and W2 have correlation p and a 0. If a = 0 then S, is the MROU model, see Chapter 3, with a cutoff at e0. If a < 0, the function (1 + aXt)1/a increases more rapidly that an exponential function. An important advantage of this approach compared to other methods for produce spikes in the spot price process is the inclusion of just one more parameter a in the model to be estimated. The deterministic seasonality function is the same as is described by equation (3.33) and its components are estimated by least-square fitting exactly as in Chapter 3. We denote the deseasonalized spot price by St, that is = S/h(t). where as before h(t) is the estimated seasonal correction. 5.2 Future price Based on the risk neutralized process (5.1)-(5.3) we calculate the Future price. Assuming a deterministic interest rate the Future price is the expected future spot price under the risk neutral measure Q, i.e. fr(t,T) h(T)E(fa(XT)IXt = x,L = 1) = h(T)f faWy’2()ex{_ (w_(:x1))2}d (5.4) where i(s, x, 1), and a(s) are the mean and variance of XT respectively. 90 Chapter 5. NLMROU model Taking y = (w — t(s, x, l))/u(s) (t, T) = h(T) f(s, x, 1) + u(s)y)e_Y2/2dy, (5.5) = 1d + a((s, x, 1) u(s)y))1e_Y2/2dy + A0 f°° *e_Y2/2dy, (5.6) d1V = 1d + ((s, x, 1) + u(s)y))1e_Y2/2dy +A0I(—d1), (5.7) with A i/ — — 1 (s, x, 1)j-j0—C , ui— — ______ u(s)ci u(s) and (.) denotes the cumulative normal distribution function. Using the fact that E(e’) = eJ2/2 with Y N(t, cr2) we are able to calculate the mean (s, x, 1), and variance u(s) of XT in (5.4). Let and s=T—t. — From previous calculations, see equation (3.9) and Appendix (6.1) E(exp{(u1,u2) (XT, LT)’}I(Xt, Li)’) = exp{M(s) +N1(s)X +N2(s)Lt} where Ni(s) = uOx(s), N2(s) = and 91 Chapter 5. NLMROU model M(s) = — 1) + (rLuio_ LLU2) (8() — 1) +( — — 2PcrxJLuao) (&2X(S) — 1) + (_J — ou + 2ouiuco ) 0( — 1) + (‘ — JLU1U2a0 + PJxJLU?ao — pJxJLUlU2’\ X + L )(8x+(s)—1). Taking u2 = 0 E(exp{u1Xr}J(X = x, L = 1)’) = exp{u1(M1s)+ Nii(s)x + N21(s) ) + Mi2(s)} (5.8) where LLO &X(8) — 1) ALL0O+ (&L(s) — 1),Mu(s) = ( ,XL Na(s) = N21(s) = cio(Ox(s) and M12(s) — — —2PXLO (&2 — 1) — 2 ___ — 1)+ 2u +2PuxuLao(0 () 1).(02L(S) (Ax + AL) Therefore the mean and variance of XT are respectively x, 1) = Mn(s) + Nn(s)x + N21(s)l 92 Chapter 5. NLMROU model and u(s) = M12(s). Unfortunately the integral in (5.7) does not admit a closed-form solution. However, this is not a significant obstacle, since these integrals can be evalu ated quickly by numerical methods. The Futures prices have been calculated for this model, though in a rather non-explicit fashion - see [78]. 5.3 Unscented Kalman filter setup and estimation procedure This model is non-linear but has Gaussian noise, so an appropriate technique is to use the Extended Kalman filter (EKF) or the Unscented Kalman filter (UKF), see Section 2.3. The starting point of our inferential procedure involves an Euler dis cretization. The model is thus evaluated at a set of discrete times {t : i = 0, 1, ..., n} such that At = t — t_1. Writing X for X as before, the Euler scheme for equations (5.1)-(5.3) can be written as: = f(Xt), (5.9) X = X_1 — x(X_1— L_1)At + —p2uXAT’/ + puxAW2, (5.10) = L_1 — AL(Lt_i — L)At + uLAW2. (5.11) In order to apply the Unscented Kalman filter, we use the state-space repre sentation for the NLMROU model. We defined the state equation as ( X N 1 0 N I 1—AxAt )xAt NI X_ L ) = ALLAt 0 1LAt ) L_1 + I ‘ —P2x”1+ puxA (5 12) uLfW? . The observations we have available, after seasonal correction, are the spot and filture prices (S F(t, T)). However, the UKF allows us to take as 93 Chapter 5. NLMROU model our measurement any transformation of these observations. To reduce non linearity in the model we therefore took our measurement equation as — ( logS N (/AW’ NVt — logF(t,T) ) + uzW2) ( log f(Xt) = log{ f(i + (s, x, 1) + u(s)y))1/ae_Y2/2dy+A0(—d1)} I usvW1 N+ (5.13) where St denotes the deseasonalized spot price, and F(t, T) is the deseasonal ized future price at time t with maturity T. In this case the noise covariances are given by — ( /1 —p2ax/M4’ç’ + puxzW N — I u%At PJXULAt N q — coy ) — PJXJL/t p and 1us/M’V’\ (cT. 0 = coy Aw) = 0 4 We set up the specific characteristics of the state-space model for the spot and future prices using the transition and measurement equations (5.12) and (5.13). The UKF parameters were set to c = 0.001, = 2, and ‘ 0. Based on this state-space formulation we are able to run the unscented Kalman filter algorithm in order to estimate the parameter set &={)\x,)L,ux,uL,L,p,a,us,uz} by means of maximum likelihood according to Section (2.5). We run Algo rithm 4 in Section 3.4, but using the UKF instead of the Kalman filter in (Step 2) to obtain 0. 94 Chapter 5. NLMROU model The log-likelihood function is calculated as £O(VLn) — [log IvkvkI + (Uk — VkIk_1) — VkIk_1)] (5.14) where EVkVk, and tikIk1 are given by (2.23) and (2.22). 5.4 Simulation results Assuming the parameters are known, and using the Euler scheme discretiza tion in (5.9)-(5.11) we simulated a path of the deseasonalized spot and future price with n 800 observations of the NLMROU model using A = 100, AL = 3, ux 1.7, L = 0.4, L = 1.9, p = 0, c = —0.4 and = 800a with = 1/250. Figure 5.1 shows the trajectories. As we can see in Figure 5.2 the UKF is able to recover the ‘hidden’ states. For simplicity we fixed the parameters p = 0, s = 0.1, and og = 0.01. The optimization method was repeated 25 times with random re-initialization for each run to obtain: = argmax£9(i). (5.15) 8e0 and we proceeded in the same way with 50 different trajectories to obtain the following results. See Tables 5.1, and 5.2: this shows that the algorithm is able to obtain quite good estimates of the parameters for simulated data. 5.5 Parameter estimation based on historical data For our empirical analysis we again use the data from the European Energy Exchange (EEX) in Leipzig, Germany. In our analysis we considered the spot price of the EEX baseload index and monthly baseload futures prices. The spot price is an equally weighted average of all 24 hourly spot prices for that particular day. Holidays and weekends have been removed from the data set. 95 Chapter 5. NLMROU model SpotS=f(X) I 101 0 100 200 300 500 600 700 8 Figure 5.1: The graph shows a simulated trajectory of the NLMROU model. n = 800 True Value Run 1 Run 2 Run 3 Run 4 Run 5 Ax 100 102.562 105.627 101.162 104.041 106.12 AL 3 3.523 2.871 3.466 3.540 3.618 clx 1.7 1.582 1.530 1.845 1.728 1.652 L 0.4 0.419 0.457 0.420 0.380 0.379 L 1.9 1.951 1.875 1.952 1.803 1.974 c -0.4 -0.402 -0.382 -0.461 -0.389 -0.406 —(6) -186.276 -179.261 -197.417 -186.248 -191.87 Table 5.1: Five different maximization runs with 800 observations. Our data comprise almost five years baseload day prices from July 1, 2002 to June 29, 2007, totaling 1267 observations. The dynamics of the spot prices for the considered period are shown in Figure 3.3. Days Future F(t,T)E(STSt) I S =f (X)tco t 400 Days MROU process X & L Days 96 Chapter 5. NLMROU model Figure 5.2: Plot of the true and estimated processes L and X of the NLM ROU model (first 80 observations). n = 800 ‘Iue value Estimator Std. Ax 100 102.4069 4.340 AL 3 2.964 0.795 x 1.7 1.786 0.216 L 0.4 0.3954 0.0921 L 1.9 1.8968 0.0216 a -0.4 -0.4009 0.038 Table 5.2: Estimation using one futures contract (n = 800). To estimate the parameters of the deterministic part IT1 + 2irtN T2 + 2irth(t)=r+/3ot+/3icos 250 )+/3cos( 5 ), tine Lep 97 Chapter 5. NLMROU model we ran the least squares method as we did in Chapter 3. See Table 3.5 for the estimated values of , , 2, r1, and r2. In the following we will analyze the three time series of the EEX market for the periods July 1, 2002 - December 31, 2004, January 1, 2005 - June 29, 2007, and the whole series. The seasonalities have been removed. For simplicity and to reduce the number of parameters to estimate, we fixed the noise parameters to be constants. s = 0.1 and crz = 0.01. In view of the low correlation found between the processes X and L in Chapter 3, we took the covariance p = 0. The results on the parameter estimates are shown in the Table 5.3. Part 1 Part 2 Whole Estimator Estimator Estimator Ax 127.599 116.036 121.557 AL 0.11178 0.753711 0.203681 ox 1.25847 0.756485 1.15815 0L 0.0834349 0.0730882 0.0901108 L 1.86049 1.77443 1.79335 c -0.438183 -0.489661 -0.414042 Table 5.3: Estimated values for the EEX market using S, and F(t, T1). The parameters A and AL relate, as in Chapter 3, to the ‘half-life’ of the mean reversion of the short-term and long-term processes. We obtained similar results for TX log 2/Ax, showing a half life of 1 — 2 days. The estimated half-lives for L differ significantly for the two periods, being about 6 years for the first period, and about 8 months for the second. While the MROU and NLMROU models give similar estimates for the speed of mean reversion of the short-term component X, the NLMROU model gives significantly slower mean reversion for the long-run process L. The estimates of the nonlinearity parameter c are quite similar for the two periods. As in the case of the MROU process, comparison of the variances VL = o/2AL, Vx = 4/2Ax show that the dynamics of X, and L contribute significantly to the long run variation of S. As before the long run process 98 Chapter 5. NLMROU model gives a somewhat greater contribution. We obtain the following results: w Period 1 0.18 0.08 Period 2 0.06 0.05 Note that because we are applying a different function (f() rather than exp{.}) to X to obtain the spot price, it does not make sense to directly compare the values of x and crL with those obtained in Chapter 3. In order to investigate if the estimated parameters make sense, we simu lated a path of the spot price, future prices and long and short term mean process that describe the model using the estimated values from Table 5.3 considering the whole data, see Figure 5.3. Empirical moments of the EEX spot price versus simulated moments (averaged over 50 simulation paths) are shown in Table 5.4. Comparing the real spot price with the price produced by the model, we see that the model is able to produce significant price spikes with values similar to those for the real data. As with the real data, the spikes tend to bunch together. As well, the model exhibits periods of high variance - compare for example the period 700 — 950 in Figure 3.5, and 750 — 950 in Figure 5.3. The estimated parameters for the long-run mean process give quite slow mean reversion, and Figure 5.3 shows a process of this kind. This model will tend to generate more spikes when the L process is large, and one sees this at the end of the simulation, when the simulated spot has many small spikes in the range from €150 on. Similar features are seen in the real data - in the period 800 to 1100 in Figure 3.5. The model does however have some defects. The first is that the skewness of the log returns is close to zero - see Table 5.4 (this is likely to be a feature of any model without jumps). Then this model is not able to capture the skewness which is present in the real data. It also appears that the model underestimates the kurtosis, compared with the real data. In spite of these problems, this model appears to offer a significant im provement over the MROU model at a quite moderate ‘cost’ - the cost being in terms of additional parameters and complications of parameter estimation. Although the model is far from perfect, its performance suggests that it is 99 40 Chapter 5. NLMROU model SpotS=f(X) I Figure 5.3: Simulation of spot and future prices price using estimated values for the whole data set. Real data (1) Sim. data (1) Real data (2) Sim. data (2) Mean 0.0019 0.0000 0.0000 -0.00037 Std Dev. 0.2214 0.60029 0.2066 0.46088 Skewness -1.3125 -0.02066 -0.0800 0.04927 Kurtosis 29.8250 4.77539 8.1288 3.99024 n data 634 634 631 631 Table 5.4: The table shows the first four moments of the logarithmic de seasonalized price returns of observed data and the average of 50 simulated trajectories. well worth considering more elaborate models of this type in the search for a good description of electricity prices. 50 30 20 Days Future F(t,T)E(STISt) I Sf(X) 100 200 300 400 500 600 Days MROU process X & L I I I 700 800 900 11 I0 500 Days 100 Chapter 6 Conclusions In this thesis, based on the specific properties of electricity we have proposed three process models that incorporate various features of power prices. We calibrate these models, using both spot and futures prices, to artificial and historical data applying three filtering methods. The first model is a two-factor linear Gaussian model. This models the log-spot price as mean-reverting process, where the mean reversion is to a sec ond (unobservable) “long run mean” process. This second process is modeled by an Ornstein-Uhlenbeck process. In this thesis (unlike the work reported in [4]) we used both the spot and future prices to estimate parameters. The space-state formulation of this model is suitable to the application of Kalman filtering techniques, and we used a maximum likelihood estimator based on the Kalman filter. This worked well for simulated data, and we then applied it to estimate parameters for the German EEX market. Simulations suggest that this model, with the fitted parameters, does fit some features of the real data. However, it definitely fails to exhibit some features of the real data, such as the jumps or price spikes seen in Figure 3.3. This defect in the first model led us to consider a second model, which incorporates jumps. We kept the same basic form of a log-spot price, and a long-term mean process, but added jumps to the log-spot price. For sim plicity we took the distribution of both the upward and downward jumps to be exponential, with possibly different rates and parameters. The standard Kalman filter cannot be applied satisfactory in this case, since the model is non-Gaussian. One alternative, which requires rather weak assumptions on the distribution involved, is the particle filter. We developed code to use the particle filter for the second (jump diffusion) model based on the kernel ap proximation of the posteriori suggested by Liu & West (2001). An empirical application on simulated data is presented to study the performance of the implemented algorithm. In general we observed that while the particle filter can work satisfactorily in estimating the unknown parameters, it is very sen 101 Chapter 6. Conclusions sitive to the particular form of dynamics for the artificial parameters used in the parameter estimation. In view of this it is hard to apply to real data. The third model we presented is an extension of the nonlinear Ornstein Uhienbeck model (NLOU) proposed in Barlow [3]. While Barlow used only one factor to describe the dynamics of the spot price in this thesis we consider a two-factor model and the same nonlinear transformations to model the spot price. The model captures the mean-reversion, jumps and spikes behavior observed in real market. The model has the advantage over jump-diffusion models that it is Gaussian. Hence we can use the unscented Kalman filter algorithm to estimate the NLMROU model. We calibrate the models to daily EEX market obtaining similar simulated trajectories with the estimated parameters. 6.1 Future work 1. In Chapter 5 we analyzed a NLMROU model of the form = h(t)f(X) where (Xe, L) is given by (5.2)-(5.3). The original model in [3] was justified by considerations of supply and demand curves, and so it might be more realistic to consider a model of the form = f(Xt + h(t)) (6.1) where h(t) is a deterministic seasonal correction. This model would have the merit of generating spikes during periods of high demand, without the need to consider different regimes, as is done in [45, 72]. The main obstacle to parameter estimation for (6.1) is that, since c is unknown, it is no longer possible to estimate h(t) by a least squares method. One possible approach is an iterative method. One would first apply the UKF to the uncorrected series (Si, F(t, T)), to obtain an initial estimate &i. One then applies least squares to the series f’(St), to obtain 102 Chapter 6. Conclusions an estimate h1. Given h1, the spot component of the measurement equation is v(’(t) = log (f(Xt + which can be calculated as in Section 5.2. Using this, and a similar expression for the futures component, one can then apply the UKF to obtain a second round of parameter estimates, and in particular an improved estimate &2 for c. Iterating, one would hope to obtain estimates for o and h(t). This procedure seems feasible to implement, but its convergence prop erties and stability are at this point not clear. 2. Improved parameter estimation using the particle filter. In this thesis we applied the Bayesian approach where an augmented state variable which includes the parameters is processed by the particle filter for the NLMROU model. We adopted the on-line estimation of parameters and state developed by Liu West (2001). The main feature of this approach is that the variance matrix A shrinks step by step and it finally converges towards 0. Hence the parameters could converges towards a wrong value because the A converges towards 0 before reaching the true parameter value. This problem cannot be avoided when we do not have prior knowledge about the parameters. That is why the method is very sensitive to the initial values of the added noise in the parameters. Along the same lines, another method that can we use is the ‘practical filter’ proposed by Polson et al. (2008) [46, 67]. Their approach is based on approximating the target posterior by a mixture of fixed-lag smoothing distributions. According to the authors, unlike the parti cle filter approaches, it provides independent samples from the target distribution, does not suffer from particle degeneracies, and handles outliers and high dimensional problems well. Meyer-Brandis T. & Tankov P. (2008) [62] comment: “Sequential filtering makes less sense when the complete series is avail able for estimation”. 103 Chapter 6. Conclusions Although computationally more intensive, the difficulties encountered in implementing the particle filter suggest this may be the correct ap proach. An approach used by Olsson et al. (2008) [65] and Wills et al. (2008) [81) is an off-line method performing maximum likelihood esti mation via the EM algorithm. An essential component in the E-step is to approximate the ‘smoothing distribution’, that is {po(utlvi:n);t = 1,.. . ,n}. In the general nonlinear and non-Gaussian case various schemes have been proposed. The fixed-lag approximation is the sim pies approach and it was first proposed in [54]. In [65] they apply the particle filter technique to smooth additive functionals based on a fixed-lag smoother. The method exploits the forgetting properties on the conditional hidden chain and is not affected by the degeneracy of the particle trajectories. Approaches of this kind are well worth investigating for jump models of the type considered in Chapter 4. 104 Bibliography [1] L. Anderson and D. Davison. A hybrid system-econometric model for electricity spot prices: Considering spike sensitivity to forced outage distributions. IEEE Transactions on Power Systems, 23(3): 927—937, 2008. [2] M.S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp. A tutorial on particle filters for on-line non-linear/non-Gaussian Bayesian tracking. lEE Trans. Signal Processing, 50(2):174—188, 2002. [3] M. Barlow. A diffusion model for electricity prices. Mathematical Fi nance, 12:287—298, 2002. [4] M. Barlow, M. Gusev, and M. Lai. Calibration of multifactor models in electricity markets. International Journal of Theoretical and Applied Finance, 7(2):101—120, 2004. [5] D.R. Beaglehole and M.S. Tenney. General solutions of some interest rate contingent claim pricing equations. Journal of Fixed Income, 1:69— 83, 1991. [6] F.E. Benth, J.S. Benth, and S. Koekebakker. Stochastic Modelling of Electricity and Related Markets. World Scientific, 2008. [7] F.E. Benth, J. Kallsen, and T. Meyer-Brandis. A non-Gaussian Ornstein-Uhlenbeck process for electricity spot price modeling and derivatives pricing. Applied Mathematical Finance, 14(2): 153—169, 2007. [8] F.E. Benth and S. Koekebakker. Stochastic modeling of financial elec tricity contracts. Energy Economics, 30(3):1116—1157, 2008. [9] T. Bjork. Arbitrage Theory in Continuous Time, 2nd Ed. Oxford Uni versity Press, 2004. 105 Bibliography [10] S. Borovkova and H. Geman. Analysis and modelling of electricity fu tures prices. Studies in Nonlinear Dynamics éi Econometrics, 10(3):Ar- tide 6, 2006. [11] M.J. Brennan and E. Schwartz. Evaluating natural resource investments. Journal of Business, 1985. [12] M. Burger, B. Klar, A. Mueller, and G. Schindimayr. A spot market model for pricing derivatives in electricity markets. Journal of Quanti tative Finance, 4:109—122, 2004. [13] 0. Cappe, E. Moulines, and T. Ryden. Inference in Hidden Markov Models. Springer, 2006. [14] R. Carmona and M. Ludkovski. Spot convenience yield models for the energy markets. AMS Mathematics of Fianance, 351:65—80, 2004. [15] C. Carraro. Square root Kalman algorithms in econometrics. Computer Sciences in Economics and Management, 1:41—51, 1988. [16] A. Cartea and G. Figueroa. Pricing in electricity markets: a mean reverting jump diffusion model with seasonality. Applied Mathematical Fiance, 12(4):313—335, 2005. [17] A. Cartea and P. Villaplana. Spot price modeling and the valua tion of electricity forward contracts. Journal of Banking and Finance, 32(12):2502—2519, 2008. [18] T. Chen, J. Morris, and E. Martin. Particle filters for state and pa rameter estimation in batch process. Journal of Process Control, pages 665—673, 2005. [19] L. Clewlow and C. Strickland. Valuing energy options in a one factor model fitted to forward prices. Working Paper, University of Sydney. [20] L. Clewlow and C. Strickland. Energy Derivatives: Pricing and Risk Management. Lacima Group, 2000. [21] M. Culot, V. Goffin, S. Lawford, S. de Menten, and Y. Smeers. An affine jump diffusion model for electricity, working paper, 2006. 106 Bibliography [22] D. Davison, L. Anderson, B. Marcus, and K. Anderson. Development of a hybrid model for electrical power spot prices. IEEE Transactions on Power Systems, 17(2):257—264, 2002. [23] J.F.G. De Freitas, M. Niranjan, and A.H. Gee. Dynamic learning with the EM algorithm for neural networks. Journal of VLSI Signal Process ing Systems, 26:119 — 131, 2000. [24] S. Deng. Pricing electricity derivatives under alternative stochastic spot price models. 33rd Hawaii International Conference on System Sciences, 4:4025—4034, 2000. [25] J.E. Jr. Dennis and R.B. Schnabel. Numerical Methods for Uncon strained Optimization and Nonlinear Equations. SIAM, 1996. [26] A. Doucet, N. de Freitas, and N. Gordon. Sequential Monte Carlo Meth ods in Practise. Springer-Verlag, 2001. [27] D. Duffie, J. Pan, and K. Singleton. Transform analysis and asset pricing for affine jump-diffusions. Econometrica, 6:1343—1376, 2000. [28] J. Durbin and S.J. Koopman. Time Series Analysis by State Space Methods. Oxford University Press, 2001. [29] C. Erlwein, F. Benth, and R.S. Mamon. Hmm filtering and parameter estimation of electricity spot price model. Revised for Energy Economics, 2009. [30] A. Eydeland and K. Wolyniec. Energy and Power Risk Management: New Developments in Modeling, Pricing, and Hedging. John Wiley & Sons, 2003. [31] H. Geman. Commodities and Commodity Derivatives. John Wiley Sons, 2005. [32] H. Geman and A. Roncoroni. Understanding the fine structure of elec tricity prices. Journal of Business, 79(2):1225—1261, 2006. [331 H. Geman and 0. Vasicek. Forwards and futures on non storable com modities. RISK, August:93—97, 2001. 107 Bibliography [34] Z. Ghahramani and S. Roweis. Learning nonlinear dynamical systems using an EM algorithm. Advances in Neural Information Processing Systems, 11:599—605, 1999. [35] R. Gibson and E. Schwartz. Stochastic convenience yield and the pricing of oil contingent claims. Journal of Finance, 3:959—976, 1990. [36] B. Gopaluni, T.B. Schon, and A. Wills. Particle filter approach to non linear system identification under missing observations with a real appli cation. Proceedings of IFAC Symposium on System Identfication, 2009. [37] N.J. Gordon, D.J. Salmond, and A.F.M. Smith. Novel approach to non-linear/non-Gaussian Bayesian state estimation. lEE Proceeding -F Radars Sonar Navig., 140(2):107—113, 1993. [38] M . Grewal and A.P. Andrews. Kalman filtering: Theory and practice using MATLAB. Wiley, 2001. [39] J.M. Griffin and S.L. Puller. Electricity deregulation: choices and chal lenges. Series in the Economics of Public Policy, 2005. [40] J. Hamilton. Time Series Analisys. Princeton University Press, 1994. [41] S. Haykin. Kalman Filtering and Neural Networks. John Wiley & Sons, 2001. [42] D. Heath, R. Jarrow, and A. Morton. Bond pricing and the term struc ture of interest rates: A new methodology for contingent claim valuation. Econometrica, 60:77—105, 1992. [43] T. Higuchi. Self-organizing time series model. Sequential Monte Carlo Methods in Practise (Doucet, de Freitas and Gordon, Eds), pages 429— 444, 2001. [44] 5. Hikspoors and S. Jaimungal. Energy spot price models and spread op tions pricing. International Journal of Theoretical and Applied Finance, 10(7):1111—1135, 2007. [45] R. Huisman and R. Mahieu. Regime jumps in electricity prices. Energy Economics, 25:425434, 2003. 108 Bibliography [46] M.S. Johannes, N.G. Poison, and J.R. Stroud. Optimal filtering of jump diffusions: Extracting iatent states from asset prices. Review of Finan cial Studies, 22(7): 2559—2599, 2009. [47] S.J. Julier and J.K. Uhlmann. Unscented filtering and nonlinear esti mation. Proceedings of the IEEE, 92(3):401—421, 2004. [48] R.E. Kalman. A new approach to linear filtering and prediction prob lems. Transactions of the ASME - Journal of Basic Engineering, D 82:35—45, 1960. [49] T. Kanamura and K. Ohashi. A structural model for electricity prices with spikes : Measurement of spike risk and optimal policies for hy dropower plant operation. Energy economics, 29(5):1010—1032, 2007. [50] N. Kantas, A. Doucet, S.S. Singh, and J.M. Maciejowski. Overview of Sequential Monte Carlo methods for parameter estimation on general state space models. 15th IFAC Symposium on System Identification (SYSID), In Press, 2009. [51] B.P. Kellerhals. Financial Pricing Models in Continuous Time and Kalman Filtering. Springer, 2001. [52] V.A. Kholodnyi. A non-markovian process for power prices with spikes and valuation of european contingent claims of power. Preprint, 2001. [53] G. Kitagawa. Self-organizing state space model. Journal of the American Statistical Association, 93(443): 1203—12 15, 1998. [54] G. Kitagawa and S. Sato. Monte carlo smoothing and self-organizing state-space models. Sequential Monte Carlo Methods in Practise (Doucet, de Freitas and Gordon, Eds), 2001. [55] 5. Koekebakker and F. Olimar. Forward curve dynamics in the Nordic electricity market. Managerial Finance, 31(6):73—94, 2005. [56] A. Leon and A. Rubia. Testing for weekly seasonal unit roots in the Spanish power pool. Modelling prices in competitive electricity markets, D. W. Bunn Ed, pages 177—189, 2004. 109 Bibliography [57] J. Liu and M. West. Combined parameter and state estimation in simulation-based filtering. Sequential Monte Carlo Methods in Practise (Doucet, de Freitas and Gordon, Eds), pages 197—223, 2001. [58] J.S. Liu and R. Chen. Sequential Monte Carlo methods for dynamic systems. Journal of the American Statistical Association, 93:1032—1044, 1998. [59] J. Lucia and E. Schwartz. Electricity prices and power derivatives - evi dence from the Nordic power exchange. Review of Derivatives Research, 5:5—50, 2002. [60] M.R. Lyle and R.J. Elliott. A simple hybrid model for power derivatives. Energy Economics, 31, 2009. [61] M. Manoliu and S. Tompaidis. Energy futures prices: term structure models with Kalman filter estimation. Apply Mathematical Fiance, 9:21— 43, 2002. [62] T. Meyer-Brandis and P. Tankov. Multi-factor jump-diffusion models of electricity prices. International Journal of Theoretical and Applied Finance, 11(5):503—528, 2008. [63] A. Misiorek, S. Trueck, and R. Weron. Point and interval forecasting of spot electricity prices: Linear vs. non-linear time series models. Studies in Nonlinear Dynamics and Econometrics, 10, 2006. [64] N.K. Nomikos and 0. Soldatos. Using affine jump diffusion models for modelling and pricing electricity derivatives. Applied Mathematical Finance, 15(1):4171, 2008. [65] J. Olsson, 0. Cappé, R. Douc, and E. Moulines. Sequential Monte Carlo smoothing with application to parameter estimation in nonlinear state space models. Bernoulli, 14, 2008. [66] D. Pilipovic. Energy Risk: Valuing and Managing Energy Derivatives. McGraw-Hill, 1998. [67] N.G. Poison, J. R. Stroud, and P. Muller. Practical filtering with se quential parameter learning. Journal of the Royal Statistical Society, Series B, 70(2):413—428, 2008. 110 Bibliography [68] D. Raggi. Adaptive MCMC methods for inference on affine stochastic volatility models with jumps. The Econometric Journal, 8(2):235—250, 2005. [69] B. Ristic, S. Aralampalam, and N. Gordon. Beyond the Kalman filter. Artech House, 2004. [70] C.P. Robert and G. Casella. Monte Carlo Statistical Methods, 2nd Ed. Springer-Verlag, 2004. [71] 5. Roweis and Z. Ghahramani. A unifying review of linear Gaussian models. Neural Computation, 11:305—345, 1999. [72] G. Schindimayr. A regime-switching model for electricity spot prices. Working Paper, 2005. [73] E.S. Schwartz. The stochastic behavior of commodity prices: impli cations for valuation and hedging. Journal of Finance, 52(3):923—973, 1997. [74] R.H. Shumway and D.S. Stoffer. Time Series Analysis and Its Applica tions: With R Examples, 2nd Ed. Springer, 2006. [75] D. Simon. Optimal State Estimation. John Wiley Sons, 2006. [76] A.F.M. Smith and A.E. Gelfand. Bayesian statistics without tears: A sampling-resampling perspective. The American Statistician, 46(2):84— 88, 1992. [77] H. Tanizaki. Nonlinear Filters: Estimation and Applications. Springer- Verlag, 1996. [78] J. Teichmann. A note on nonaffine solution of term structure equations with applications to power exchanges. Mathematical Finance, 15(1):191— 201, 2005. [79] A. van der Merwe, R. Doucet and N. de Freitas. The unscented particle filter. Technical Report CUED, 380, 2000. [80] P. Villaplana. Pricing power derivatives: a two-factor jump-diffusion approach. Business Economics Series (Working Paper), 2003. 111 Bibliography [81] A. Wills, T. Schon, and B. Ninness. Parameter estimation for discrete- time nonlinear systems using EM. Proceedings of the 17th IFAC World Congress, Seoul, Korea, 2008. [82] L. Xiong. Stochastic models for electricity prices in Alberta. University of Calgary, M.Sc. thesis, 2004. 112 Appendix A An Affine Jump-Diffusion process (AJD) is a jump-diffusion process for which the drifts and covariances and jump intensities are linear in the state vec tor U. Duffie et al. (2000) [27] show that AJD processes are analytically tractable in general. Let U be a strong Markov process with realizations {U, 0 t < oo} in some state space D C R’2, which solves the following stochastic differential equation U U0 + f (U5,s)ds + f u(U, s)dW + Z. (A-i) The jump behavior of U is governed by m types of jump processes. Each jump type Z is a pure jump process with a stochastic arrival intensity t) for some : (D, t) —* R and jump amplitude distribution v on R”, where v only depends on time t. The functions : (D, t) F—* and u (D, t) —* are assumed to be Lipschitz continuous in order to guarantee that (A-i) has a unique solution. The process W.. is a standard Brownian motion in IRA. The process U defined by (A-i), is said to be an affine jump-diffusion process if (U,t) = K0(t) +K1(t)U, u(U, t)u’(U, t) = Ho(t) + H(t)Uk, = where for each 0 t < no, K0(t) E RTh, K1(t) R”, H0(t) RT1 and is symmetric, H1(t) e Also, for k = 1,... ,n, H(t), defined to be 113 Appendix A the matrix obtained by fixing the third index of H1 (t) to be k is in R and is symmetric. Finally l(t) R Notice that given an initial condition U0, the tuple (K0,K1,H0,H1, l) can be used to determine a transform ‘J : x [0, cc) x [0, cc) x D C of UT conditional on Ut, 0 <t <T, defined by ‘I(u,t,T,U) := E[exp{u . U}U] (A-2) where E denotes the expectation under the distribution of UT determined by (K0,K1,H0,H1, 1w). If we suppose (K0,K1,H0,H1, 1) is well-behaved at (u, T), then the transform ‘V of Ut, 0 < t <T, defined by (A-2) exists and is given by: P(u, t, T, U) = exp{M(u, t, T) + N(u, t, T) U}. (A-3) Here M(.) and N(.) satisfy the following complex-valued Riccati equations: aM(u,t,T) = -A(N(u,t,T),t), M(u,T,T) = 0, aN(u,t,T) = -B(N(u,t,T),t), N(u,T,T) = u, where, for any c e C’, A(c, t) = Ko(t) . c + c’Ho(t)c + lj(c) — 1), B(c,t) = K1(t)’c+ c’H1(t)c. Here (c) is the “jump transform” for the i-th jump. It is given by = f exp{c. U}dv(U) (A-4) whenever the integral is well defined. 114 Appendix A We define the extended transform 1 : R x C x [0, oc) x [0, oc) x D H—* C of UT conditional on U, by u, t, T, U) E[(v UT) exp{u. UT}IUt]. (A-5) Given sufficient regularity the “extended transform” P can be computed by differentiation of the transform \J1 Hence 1(v, u, t, T, U) = ‘I’(u, t, T, U){C(t) + D(t) U} (A-6) where ‘I’ is given by (A-3), and C(.) and D(.) satisfy 0 = —Ko(t)’D — N’Ho(t)D — 1o(t)Vo(N)D, (A-7) = —K1(t)’D — N’H1(t)D, (A-8) with the boundary conditions C(T)=0, D(T)=v. (A-9) Here V(c) is the gradient of (c) with respect to c C. 115 Appendix B In this appendix we solve the ODEs (3.13)-(3.15) which arise in the calcula tion of futures prices in the AJD model. We begin with the ODEs arising from the MROU model. Recall that (B-i) = —xNi + ALN2, (B-2) = — (uN + uN) —N12pUXUL. (B-3) We begin with (B-i). We have N1 = N1(t, T) and Ni(ui,T,T)=ui. (B-4) If we fix T, so regard N1,N2 as function of t only, then (B-1)-(B-3) are ODEs and (B-4) has solution N1 = etA. Since A = uie_T, therefore Ni(ui,t,T) = uie>(t_T). (B-5) We now treat (B-3): = + N2, N2(u,T,T) = u2. (B-6) i i6 Appendix B Substituting for N1, and multiplying by a factor t we obtain dN2 — /J)LN (B-7) So (N2) dN2 dt [L-—/LALN, N2 PLN2, dp. _ — —ALdE, I1 ij, e. Substituting for in equation (B-?), we obtain d(etN2) _Axuie(t_T)e_t (B-8) dt — _Axuietx_e_T (B-9) so, AXU1 etx_e_T + C, (B-b)eLN2 — — AL e(t_T) + Cet. (B-il)N2 — — AL Now, solving for C An1 ‘u2=N(,T,T)rr_ +CeT (B-12) — AL 117 Appendix B therefore AXU1 e_T. (B-13)C = U2CL + — Substituting C in equation (B-li) gives N2(u,t, — eAx(t_T) + et(u2_T + Xl e_T) (B-14)T)—— Axu1 _____ AX—AL AX—AL Axui Xl et_T). (B-15)e(t_T) +u2e)(t_T) + A — AL= AX — AL Finally we solve (B-3). = —ALN2— + uN) —N12pJXJL (B-16) with M((ui,u2),T,T) 0. Let c —Au1/( x — AL). Replacing the solution for N1 and N2 given by equations (B-5) and (B-15) in (B-16) we have dM = ALL(e + u2e (t—T) — e(t_T)) — [uue2(t_T)dt 2 +(aet_T) + u2e (t—T) — ce)t_T))2] = _ALLe (t—T) — ALLU2e (t_T) + ALLe(t_T) _°e2_T) — U e2A_T) — _Je2t_T) — u u2etT) +uu2e(t_T) +L )(t—T) — puxuLcu1e (t—T) +puxuLulu2e +L)(tT) + puxuLu1e +L)(t-T) 118 Appendix B Integrating both sides gives _____ x(t—T) — LLU2 t_T) + e(t_T)e le2At_T) — _____ — 2A(t—T) — 4\x 4’\x e2(t_T) — U + U2e2t_T) — 4AL u2a__e(t_T) — PxL1 (x+AL)(t-T) + PUXULU1 Xx+L)(t-T) + c.e e Since M((ui,u2),t,T) 0, M = —1)— ALLU2( ____ __ AL(t-T) — 1) + (et_T) — 1) UxU1(e2A(t_T) — 1) — UL (e2At_T) — 1) — ULU2(eL(t_T) — 1) 4’\x 4AX 4’L — (e (t—T) — 1) — (e(t_T) — 1 + (e2At_T) — 1) 4)’L 2-AL + (e(t_T) — 1) — puxuLaUl (e(t_T) — 1) )X+AL PUxULQU1U2(e(Ax+L)ct_T) — 1) PUxULclu1 X + L + X + L (et_T) — 1). Now, taking (ui,u2)= (1,0) and defining m = —)‘x/’x — )L) we obtain 119 Appendix B M(t,T) = mi(e (t_T)_,) +m2(et_T)_1) +m3(e2(t_T)_1) +m4(et_T) — 1) +m5(e2)t_T) — 1), (B-17) where 2 2 m2u m2 Lrn m4 = m crL + mpcrx L m5 = — _______ AX+AL AX+AL 4AL ALLm = — (-k- + m2 mpuxaj\m1=— m3Ax 4Ax 4Ax + 2Ax ) We now turn to the case of the jump model. Here N1 and N2 are as before, but M has two more components due to the jumps. We can write M(t) = M(t) + M1(t), where M(t) is given by (B-17) and M1(t) satisfies: dM1 N1 Mi(u,T,T) = 0. (B-18) dt — N1’ Integrating both sides we have — f— ,u_uie(t_T)dt (B-19) I_ ‘= dt. (B-20)%eAxT_AXt — 1 Ui Applying the formula I dx 1 J = —[cx — ln(a + be)], (B-21)a+be ac with a = 1, b = (?7ueT)/Ul and c = —A we get 120 Appendix B 1 M1 = — in (—i + ((ueT)/ui)e_t)] + C. . (B-22) Since Mi(u,T,T) = 0, C = T + in(—1 + (B-23) U1 Therefore M1 = —(t—T)—-----1n(—1+ x \ nieAx(t_T)) + Ax \ ni) = —(t — T) — —-- (in (_uie(t_T) +N I—u1 ______________ — in AX uieAx(t_T) ) = in 1 (u — uie)(t_T) — u—U1 )• 121
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Filtering and parameter estimation for electricity...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Filtering and parameter estimation for electricity markets Molina-Escobar, Alberto 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Filtering and parameter estimation for electricity markets |
Creator |
Molina-Escobar, Alberto |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | The growing complexity of energy markets requires the introduction of in creasingly sophisticated tools for the analysis of market structures and for the modeling of the dynamics of spot market and forward prices. In order for market participants to use these markets in an efficient way, it is important to employ good mathematical models of these markets. This has proved to be particularly difficult for electricity, where markets are complex, and ex hibit a number of unique features, mainly due to the problems involved in storing electricity. In this thesis we propose three models for electricity prices. All are multifactor models, that is, as well as an observable spot price they assume the existence of an unobservable long term mean’ process. The introduction of such additional processes helps to explain the relation between spot and futures prices. In the first part of the thesis we introduce a two factor Gaus sian model for prices. Using the Kalman filter, and based on both spot and forward prices, we successfully estimate parameters for simulated data. We then estimate parameters for the German EEX market, and compare our fitted model with the observed prices. We find that this model does capture some features of the EEX market, but it fails to exhibit the price spikes which are a prominent feature of true spot prices. We therefore introduce a second model, which includes jumps. The inclusion of jumps has the potential to give a better explanation of the behavior of electricity prices, but it creates difficulties in the estimation of parameters. This is because as the model noise is non-Gaussian the Kalman filter cannot be applied satisfactorily. We implement the particle filter adopting the Liu & West approach for the jump model. This method allows us to identify the hidden process in the model, and to estimate a small number of parameters. The third model is a new model for electricity prices based on the inverse Box-Cox transformation. This model is non-linear with Gaussian noise, and can generate price spikes using fewer parameters than a multi-factor jump-diffusion model. In this context, we successfully applied the Unscented Kalman filter to estimate the parameters. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-03-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 3.0 Unported |
DOI | 10.14288/1.0069305 |
URI | http://hdl.handle.net/2429/21736 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_spring_molina-escobar_alberto.PDF [ 1.98MB ]
- Metadata
- JSON: 24-1.0069305.json
- JSON-LD: 24-1.0069305-ld.json
- RDF/XML (Pretty): 24-1.0069305-rdf.xml
- RDF/JSON: 24-1.0069305-rdf.json
- Turtle: 24-1.0069305-turtle.txt
- N-Triples: 24-1.0069305-rdf-ntriples.txt
- Original Record: 24-1.0069305-source.json
- Full Text
- 24-1.0069305-fulltext.txt
- Citation
- 24-1.0069305.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0069305/manifest