Pricing Energy Contracts under Regime Switching Time-Changed models
PPRICING ENERGY CONTRACTS UNDER REGIMESWITCHING TIME-CHANGED MODELS
KONRAD GAJEWSKI AND SEBASTIAN FERRANDO AND PABLO OLIVARES
Abstract.
The shortcomings of the popular Black-Scholes-Merton (BSM)model have led to models which could more accurately model the behavior ofthe underlying assets in energy markets, particularly in electricity and futureoil prices. In this paper we consider a class of regime switching time-changedLevy processes, which builds upon the BSM model by incorporating jumpsthrough a random clock, as well as randomly varying parameters accordingto a two-state continuous-time Markov chain. We implement pricing methodsbased on expansions of the characteristic function as in [6]. Finally, we esti-mate the parameters of the model by incorporating historic energy data andoption quotes using a variety of methods. Introduction
A main goal of the paper is to develop a methodology to price European option’scontracts on electricity and future oil prices. The approach is based on Fourierexpansions and implements models that capture specific stylized features of theunderlying assets such as stochastic volatility and random jumps. In particular,we consider a switching time-changed Levy process as an alternative to the BSMmodel and implement a pricing algorithm based on the expansion of its character-istic function as considered in [6].We compare the prices we obtain with those obtained via a computationallycostly, but accurate, Monte Carlo method and study sensitivities with respect torelevant parameters, e.g. maturity, strike price and initial price. Moreover, we con-trast prices under the regime switching model to those given by the Black-Scholesequation and show that the prices agree when the switching model is reduced tothe Black-Scholes model.We rely on the Esscher transformation, see [8], to obtain an equivalent martin-gale measure(EMM) and in order to work on a risk-neutral setting. To calibrateparameters, we use market option prices and minimize the mean squared error. Onthe other hand, and in order to estimate parameters, we implement the methodof moments, minimum distance estimation and maximum likelihood estimationtechniques. For simplicity an Expectation Maximization algorithm (EM) is notconsidered.Although most of these elements have been previously implemented, to the bestof our knowledge, the combination consisting of the selected model class, the pricingmethodology, the risk-neutral framework and the estimation/calibration approachhave not been studied before in energy markets or elsewhere. It is worth noticing a r X i v : . [ q -f i n . P R ] M a y KONRAD GAJEWSKI AND SEBASTIAN FERRANDO AND PABLO OLIVARES that non-switching models with Levy noises have been introduced in [3]. On theother hand, results for switching Levy models and their characteristic functions canbe found in [4].Many financial time series, including commodity futures seem to exhibit dra-matic breaks in their behaviour, for example in the events of political changes orfinancial crises. Different intervals sharing similar characteristic can be groupedtogether under a single regime. Models that can capture such behaviour are regimeswitching Levy. Under such a model, the process switches randomly between dif-ferent Levy processes according to an unobservable Markov chain. The regimeswitching time-changed Levy process is a pure jump process which captures twokey features of the market: the existence of regimes and price jumps.The organization of the paper is the following: Section 2 introduces the regimeswitching time-changed Levy model, we then derive its characteristic function underGamma and Inverse Gaussian subordinators. In section 3 we discuss Monte Carloand Fourier Cosine pricing methods. For Monte Carlo, we develop an algorithmto simulate trajectories of the process, as well as for pricing European call optionsby simulating many regime switching processes simultaneously. In section 4 weuse calibration and various methods to estimate values of model parameters fromoption quotes and historical prices of oil and electricity commodities.2.
Model, contract and characteristic function
Let (Ω , A , ( F t ) t ≥ , P ) be a filtered probability space verifying the usual condi-tions. For a stochastic process ( X t ) t ≥ defined on the space filtered space abovethe functions ϕ X t ( u ) and Ψ X ( u ) = t log ϕ X t ( u ) define its characteristic functionand characteristic exponent respectively. When the process has stationary and in-dependent increments the later does not depend on t . A T is the transpose of matrix A unless it is specified differently.Let the process { S t } t ≥ represent the price of the underlying asset at any time t > X t = log S t represents the logarithm of the prices.We define a continuous-time Markov chain { s t } ≤ t ≤ T driving the changes betweenregimes with the state space E = { , } . The switching times are described by asequence of independent random variables ( τ j ) j ∈ N in a way that:lim t → τ − k s t (cid:54) = s τ k for k = 1 , . The infinitesimal generator matrix of the continuous-time Markov chain is given by(1) Q = (cid:20) − λ λ λ − λ (cid:21) . Hence: P { s t + h = 2 /s t = 1 } = λ h + o ( h ) P { s t + h = 1 /s t = 2 } = λ h + o ( h )Next, consider a collection of independent subordinators L j = { L jt } t ≥ for j ∈ E ,where each subordinator L j is also independent of each process X i , for i, j ∈ E .Each subordinator is characterized by two parameters α j , β j > L j is a pure RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS3 jump process and each process X j has almost surely continuous pathsWe define the collection of time-changed Levy processes Y j = { Y jt } t ≥ where Y jt = µ j L jt + σ j B L jt . with µ j ∈ R and σ j > Z = { Z t } t ≥ as:(2) Z t ≡ Y s t t where Y s t t = µ s t L s t t + σ s t B L stt or, in differential terms: dZ t ≡ dY s t t where dY s t t = µ s t dL s t t + σ s t dB L stt The regime switching time-changed Levy process Z is assumed to be the log-priceprocess of the underlying asset and the stochastic process of the asset price itself { S t } t ≥ is defined as: S t = S exp( Z t ) . For simplicity we assume that the process will always start out at state 1 withprobability 1.Following [4], for the process Z defined in equation (2) along with a Markov chainwith the infinitesimal generator matrix Q defined in equation (1), its characteristicfunction is given by:(3) ϕ Z ( u ) = exp( iuy ) E Q { [1 ,
1] exp( t Φ( u ))[1 , T } where y = log S and Φ( u ) is the matrix:Φ( u ) = (cid:32) − λ + Ψ L (1) t ( µ u + iσ u ) λ λ − λ + Ψ L (2) t ( µ u + iσ u ) (cid:33) . Notice that conditionally on s t = j the characteristic function of Z j , i.e. thecharacteristic function of Z conditionally on s t = j , is: ϕ Z t ( u ) = ϕ L jt ( µ j u + 12 iσ j u ) . To compute the exponential matrix Φ( u ) we use a scaling and squaring algorithm ,see [10], based on the following approximation: e Φ( u ) = ( e − s Φ( u ) ) s ≈ r m (2 − s Φ( u )) s , where r m ( x ) is the [ m/m ] Pade approximant of e x and the nonnegative integers m and s are chosen in such a way as to achieve minimum error at minimal cost. Atable of errors as a function of s and m is given in [1]. The [ k/m ] Pade approximantfor the exponential function is: r km ( x ) = p km ( x ) /q km ( x )where: p km ( x ) = k (cid:88) j =0 ( k + m − j )! k !( k + m )!( k − j )! x j j ! , q km ( x ) = m (cid:88) j =0 ( k + m − j )! m !( k + m )!( m − j )! ( − x ) j j ! . KONRAD GAJEWSKI AND SEBASTIAN FERRANDO AND PABLO OLIVARES
The discounted price process is denoted ˜ S = ( ˜ S t ) t ≥ where ˜ S t := exp ( − rt ) S t . Wehave that under an EMM Q , the discounted price process ˜ S is a martingale under Q if and only if the following equation is satisfied:(4) Ψ Z ( − i) = r. See [14] for details.
Example 2.1.
Case of Inverse Gaussian and Gamma subordinators.Inverse Gamma and Gamma subordinators have been studied in [2] and [12] .When the subordinator L j is an Inverse Gaussian process with shape parameter α j and rate parameter β j , we have: (5) Ψ Z j ( u ) = α j ( (cid:114) µ j u + 12 iσ j u ) + β j − β j ) , α j > , β j > , j = 1 , . In Figure 1 the real and imaginary parts of the characteristic function and thecharacteristic function of Z = { Z t } t ≥ with an Inverse Gaussian subordinator areshown. Parameters are obtained from estimating procedures for future oil prices asexplained in Section 4. A similar result is obtained under a Gamma subordinator.When the subordinator L j is a Gamma process with shape parameter α j and rateparameter β j , we have: (6) Ψ Z j ( u ) = − α j log (cid:16) µ j u + iσ j u ) + β j β j (cid:17) , α j > , β j > , j = 1 , To have discounted prices being a martingale, according to equation (4): ψ jGamma ( − i ) = t − log (cid:34)(cid:32) i µ j u − ( σ j ) u β j (cid:33) − α j t (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u = − i = r leading to: µ j = − β j (exp( − rα j ) −
1) + ( σ j ) When the process Z j is a time-changed process subordinated by an Inverse Gaussianprocess with parameters α j , β j , the characteristic function is given by equation (5)and therefore for each state j ∈ E we solve for µ j : ψ jIG ( − i ) = t − log (cid:2) exp( − α j t (cid:16)(cid:114) − i µ j u + ( σ j ) u β j ) − β j (cid:17) ) (cid:3)(cid:12)(cid:12) u = − i = r It leads to: µ j = 12 [( β j − rα j ) + ( β j ) ] + σ j . Holding all the other parameters constant, the drift verifies equation (4). The valueof µ j is such that the j -th discounted price process, when the subordinator is anInverse Gaussian process, is a martingale. RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS5 (a)(b)
Figure 1.
The real part of the function ϕ Z t ( u ) under an IG sub-ordinator (top) and its imaginary part(bottom). Parameters areobtained from estimating procedures for future oil prices as ex-plained in Section 4.3. Pricing European options under switching Levy models
We turn back to pricing, consider a European call contract with maturity at T and strike price K . Its payoff, written in terms of the log-returns, is: h ( Z T ) = ( S e Z T − K ) + = K ( exp ( x + Z T ) − + where x := log( S /K ).The price of the contract at a time t < T and x = log( S t K ) is denoted as v ( x, t ) andverifies:(7) v ( x, t ) = e − r ∆ t E θ [ v ( y, T )] = e − r ∆ t (cid:90) R v ( y, T ) f Z T ( y | x ) dy, Notice that v ( y, T ) is the payoff at maturity time T and y := log( S T /K ).The value ∆ t = T − t is the time to maturity and r is the risk-neutral interest rate.The function f Z T ( y | x ) is the probability density function (p.d.f.) of Z T given the KONRAD GAJEWSKI AND SEBASTIAN FERRANDO AND PABLO OLIVARES value x = log( S /K ). E θ is the expectation value operator with respect to an EMM Q θ , θ ∈ R whichis determined by an Esscher transform of the historic measure P . See [8] for arationale in terms of a utility-maximization criteria.For a stochastic process ( X t ) t ≥ the latter is defined as:(8) d Q θt dP t = exp( θX t − t l X ( θ )) , ≤ t ≤ T, θ ∈ R where P t and Q θt are the respective restrictions of P and Q θ to the σ -algebra F t . We define by ϕ θX t , Ψ θX t and l θX ( u ) respectively the characteristic function,characteristic exponent and moment generating function of a process ( X t ) t ≥ underthe probability Q θ obtained by an Esscher transformation as given in equation (8).We follow a pricing approach via Fourier- Cosine Series expansion of the p.d.f. f Z T .The method has been proposed in [6].The solution to (7) is obtained by expanding the p.d.f. f Z T ( ./x ) in terms of aFourier basis under Gamma and Inverse Gaussian subordinators introduced in theprevious section .The domain of integration is truncated to a finite interval [ a, b ] for the purposes ofnumerical integration, for a discussion about selecting the truncation interval andits associated error, see [11].The Fourier-cosine expansion of f Z T is given by:(9) f Z T ( y | x ) = ∞ (cid:88) k =0 A k ( x ) cos (cid:16) kπ y − ab − a (cid:17) , where the first term of the summation is weighted by one-half.The coefficients of the Fourier expansion, denoted by A k ( x ) , are approximated by: A k ( x ) = 2 b − a Re (cid:110) ϕ Z T (cid:16) kπb − a /x (cid:17) exp (cid:16) − i kπ ab − a (cid:17)(cid:111) . Hence, substituting (9) into equation (7) we have: v ( x, t ) = 12 ( b − a ) e − r ∆ t ∞ (cid:88) k =0 A k ( x ) V k , where V k is the Fourier coefficients of v ( y, T ) given by: V k := 2 b − a (cid:90) ba v ( y, T ) cos (cid:16) kπ y − ab − a (cid:17) dy. In particular, for a European call option we have: V Callk = 2 Kb − a (cid:90) b ( e y −
1) cos (cid:16) kπ y − ab − a (cid:17) dy. For a European put option denoted by V P utk a similar expression is found. Finally,truncating the infinite series to N terms, we obtain: v ( x, t ) ≈ e − r ∆ t (cid:88) ≤ k 1) 18.9554 9.31 19.0401 0.0234(2 , 1) 19.9612 17.54 20.3456 0.1433(1 , 2) 17.9942 10.01 18.2164 0.1339(2 , 2) 19.0166 11.40 18.5523 0.193 Table 1. Comparison of European Call option Payoffs usingMonte Carlo Simulation and Fourier-Cosine Pricing, as well astheir computational times.In the parametric set considered (the next section discusses how to estimate theparameters), the Fourier-cosine method is on average 100 times faster than a stan-dard Monte Carlo approach, while producing similar level of precision.On the other hand, it is found that the difference between pricing European calloptions using Monte Carlo and Fourier-Cosine pricing remains constant for differ-ent strike prices. The error however grows linearly for increasing maturity times.To compare with the price obtained via Monte Carlo we discuss the simulationprocedure employed.It is well-known that the Markov chain { s t } t ≥ spends independent and exponen-tially distributed random times between regime transitions. The k − th switchingtime is determined by: τ k = k (cid:88) j =1 ∆ τ j KONRAD GAJEWSKI AND SEBASTIAN FERRANDO AND PABLO OLIVARES where ∆ τ j are the times between regime changes. The parameters of the exponen-tial random variables alternate between λ and λ .The differential equation is approximated numerically through finite differences us-ing the Euler-Maruyama Method. Hence, if at time t the process Z is under regime j ∈ E , the increment of the process Z during the time interval [ t, t + ∆ t ) is givenby: ∆ Z jt = µ j ∆ L jt + σ j (cid:113) ∆ L jt N (0 , . Notice that, for small ∆ t , the process remains under regime j with a probabilityclose to one. (a)(b) Figure 2. Trajectories of a switching time-changed model. Pa-rameters: T = 1 , µ = 0 . , µ = − . , σ = 1 , σ = 5 , α = α = 0 . , β = 0 . , β = 10 , λ = 5 , λ = 2 . , λ = λ = 4, λ = , λ = 4Figure 2(a) shows a single realization of a switching time-changed Levy process(top) under an IG subordinator, as well as the underlying Markov chain (bottom). RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS9 Figure 2(b) shows trajectories under a Gamma subordinator.We devise an algorithm which computes the payoff of a European Call option bysimulating m independent realizations of the process Z simultaneously and thenestimating the price according to: v ( x, T ) (cid:39) exp( − rT ) 1 m m (cid:88) j =1 ( S exp( Z j,T ) − K ) + where Z j,T is the j -th simulated log-price at maturity.We can also estimate the price using confidence intervals. The confidence interval is (a)(b) Figure 3. Monte Carlo price and its confidence interval vs thenumber of simulations under Gamma (top) and an IG (bottom)subordinators. Parameters are the same than in previous figureexcept β = 0 . , β = 0 . h ± z . s √ m , where ¯ h is the sample mean of the simulated payoffs, s is the sample standarddeviation, m is the sample size and z . is the 95% -percentile of the normal dis-tribution.The confidence interval decreases as the number of simulations approaches infinity,however in the case of the Inverse Gaussian subordinator, the interval is larger ateach simulation because Inverse Gaussian random variables have a larger variancethan Gamma random variables when β < m = 10 , the confidence interval when the subordinator is Inverse Gaussianis [18 . , . . , . (a)(b) Figure 4. European call payoff as a function of T and K undertwo different subordinators with risk neutral drift. The other pa-rameters are identical for each figure: σ = 0 . , σ = 0 . , α = α = 0 . , β = 1 , β = 1 . , λ = 2 . , λ = 1. For an IG sub-ordinator µ = 0 . , µ = 0 . µ = − . , µ = 0 . RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS11 Figure 4 indicates the behaviour of the price of a European call option for aregime switching time-changed Levy process model under Inverse Gaussian (top)and Gamma subordinators (bottom). Both payoff models are monotone in T and K . Moreover as T increases, the expected payoff increases. For K >> S theprobability that S T ≥ K is very small, therefore the payoff is close to 0 . (a)(b) Figure 5. Payoff as a function of parameters α , α . Payoff as afunction of parameters β , β In each figure parameters between states are held constant, except the onesindicated in the graph: µ = µ = 0 . , α = α = β = β = 0 . , σ = σ =0 . , λ = λ = 0 . , r = 0 . , T = 1 , S = 20 and K = 1 . Setting the parameters λ = λ implies the process spends an equal amount of time in each regime, onaverage. The only exception made is in Figure 5, where λ = 1000 and λ = 1 / β , β , the payoff approaches an asymptote because Gamma andInverse Gaussian random variables both have mean α/β, which approaches infinity (a)(b) Figure 6. Payoff as a function of parameters α , β . Payoff as afunction of parameters λ , λ for β → + . In Figure 6 changes in the price with respect to the intensity parametersand the parameters of the underlynig subordinators are shown.Parameters Reduced Switching Levy Model Black Scholes formula( T, K, r, σ ) ( N = 10 )(1,1,0.04,0.5) 19.0463 19.0392(3,1,0.1,1) 19.2955 19.3139(2,30,0.5,0.001) 8.963608 8.96361 Table 2. European call option payoff comparison between theBlack Scholes formula and Monte Carlo simulation of the reducedswitching Levy process at different parametersNotice that we can reduce the regime switching Levy processes to the Black-Scholes model by defining the subordinator L t = t and setting the parameters RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS13 equal across both regimes. Setting the drift such that the discounted Black Scholesprice process is a martingale µ = r − σ . we find that the price by a Fourier-cosineapproach is consistent with the price given by the Black-Scholes formula. See Table2. (a) (b) Figure 7. Daily price series and log-returns for WTI futures (left)and log-returns of WTI futures in NYMEX. Source: BlooombergTerminal, April 20184. Parameter estimation and numerical pricing We implement two approaches of fitting the parameters of the underlying modelto financial historical data: calibration and estimation depending when option’sprices or the underlying electricity and oil prices are considered. In a calibrationapproach, the parameters are fitting by minimizing the quadratic error between theprices obtained numerically and option quotes. The option quotes are taken fromBloomberg’s data base for a variety of strike prices and times to maturity. The pa-rameters are fitted using European call option quotes of West Texas Intermediate(WTI) crude oil.In a parameter estimation, we implement the following techniques based on histori-cal prices of oil and electricity: method of moments, minimum distance method andmaximum likelihood estimation, combined with empirical estimation of the switch-ing parameters. Specifically, we use daily historical NYMEX WTI crude oil futures(11-16-2012 to 06-05-2018) and IESO Ontario (Canada) Zone 24H electricity av-erage spot prices (06-06-2008 to 06-05-2018). Spikes and stochastic volatility areobserved in the series, Similar phenomena have been reported in other electricitymarkets, see [3, 7]. We assume that there are 250 trading days in a year with eachtrading day corresponding to ∆ t = 1 / 250 of unit time.Figures 7 and 8 plot the historic WTI oil futures and average Ontario electricityprices, as well as their log-returns. Electricity spot prices sometimes move belowzero, implying a surplus of electricity produced during low demand. Because elec-tricity produced by power suppliers must be consumed immediately, the supplier (a)(b) Figure 8. Price series and log-returns for Ontario daily averageelectricity spot price. Source: Blooomberg Terminal, April 2018.pays wholesale customers to buy the surplus energy, see [15]. All negative pricesare arbitrarily modified to CAD $0 . 01 for estimation purposes.We also compare the empirical density function of the log-returns of each commod-ity to the normal distribution under the same mean and variance parameters as thehistorical log-return process. To this end we implement a kernel smooth technique.Hence, the estimated p.d.f. is:ˆ f ( x ; θ ) = 1 nh n (cid:88) j =1 K (cid:16) x − z j h (cid:17) , where the function K is the Gaussian kernel.As the p.d.f. f ( x k ; θ ) of the log-returns is not available in a close form we simu-late the model under a parametric set θ to get the empirical p.d.f. ˆ f ( x k ; θ ) usinga kernel smoothing technique, see for example [9], and continue the optimizationprocedure.The value of h, is chosen to equal Silverman’s quantity h = 1 . σn − / , where σ RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS15 is the standard deviation of the log-return series, see [13]. See Figure 9 for empir-ical p.d.f.’s corresponding to WTI futures (left) and Ontario electricity log-returnprices(right).The kurtosis of WTI future and Ontario electricity log-return series are respec- (a) (b) Figure 9. Empirical density functions vs. normal distributionsoil log-returns.electricity log-returnstively 6.34 and 23.947, much larger than that of the normal distribution, suggestingthe presence of heavier tails and extreme behavior. Significant negative skewnessis also reported on both series, respectively -0.1314 and -0.1377In our model the parameters are described by vectors θ j = ( µ j , σ j , α j , β j ) for j = 1 , , while λ and λ reflect the parameters of the hidden Markov chain.Hence the times the chain remains in regime j are independent and exponentiallydistributed with mean λ j = λ jk , with k = mod (2) + 1. We set Θ j ⊂ R to bethe set of all feasible parameters for θ j . We assume that the two sets of parame-ters belong in different parameter spaces i.e. Θ (cid:54) = Θ . The two parameters of thesubordinator and the diffusion coefficient are required to be positive, therefore weadd the natural constraints σ j > , α j > , β j > . We assume that the duration of the j -th observed historic regime is the most prob-able value i.e. it is equal to the expectation value λ j . For each commodity, the j -thholding-rate parameter is given by: λ j = total number of days in regime j number of occurrences of regime j where it is assumed that there are 250 trading days every year.By simple inspection of the log-return process data of oil futures we set the processto be in regime one between 11-16-2012 and 11-16-2014 as well as between 02-06-2017 and 06-05-2018; otherwise, we assume that the process is in regime two. Inthe case of the log-returns of electricity spot prices; we set the process be in regimetwo whenever the absolute value of the log-returns exceeds 3 and in regime oneotherwise.Table 3 shows the average time and daily standard deviation in the two regimes of the WTI futures and Ontario electricity series. We also include the variance of thelog-returns within each regime; the different orders of magnitude between regimesjustifies the use of a switching model.By having defined the location of the regime changes and therefore estimated theCommodity ˆ λ ˆ λ St. Dev. (regime 1) St. Dev. (regime 2)Oil 0.900 3.80 0.0052 0.0148Electricity 0.2618 0.0081 0.6020 6.1086 Table 3. Holding-rate parameters estimation for each commodityas well as the variance in each regimevalues of λ , λ , the historic log-returns are separated into two sets of data, onecontaining all the data points for each regime.To calibrate the parameters within each regime we minimize the mean square error between the numeric option payoffs and European call option quotes. When theoption is out of the money the option price is obtained using a Monte Carlo approachbecause the Fourier Cosine method exhibits significant error.Thus, the objective function in regime j is J ( θ j ) = (cid:115) n (cid:88) T,K ( V ( T, K ; θ j ) − ˆ V ( T, K ; θ j )) , j ∈ { , } . where, by a convenient change in notation to emphasize the dependence on theparameters we write V ( T, K ; θ j ) and the option quotes ˆ V ( T, K ; θ j ), taken over arange of strike prices K and maturity times T .The optimal parameter is ˆ θ j = arg min θ ∈ Θ j J ( θ j ). It is calculated using a gradientdescent method.The stopping criteria is taken to be step tolerance, taken to be equal to 1e-10. Thek-th step tolerance is a lower bound on the size of the step || θ t − θ t − || . The solverstops if the stopping criteria is reached, or if the maximum number of iterations(fixed to 1000 steps) is exceeded. Different initial starting points are found to givesimilar estimation of the parameters. Table 4 gives the estimated calibration in thecase when the subordinator is a Gamma process or an Inverse Gaussian process. Asexpected, the volatility is higher in regime two in the case of both subordinators.Commodity (subordinator) ˆ µ ˆ σ ˆ α ˆ β Oil log-return Regime 1 (Gamma) -0.03387 0.0030 2.640710 1.007e-8Oil log-return Regime 2 (Gamma) -0.01445 1.116184 2.56567e-5 10.32441Oil log-return 1 (Inverse Gaussian) -0.04976 0.130011 0.24788 92.6926Oil log-return 2 (Inverse Gaussian) -0.04950 0.515891 8.531e-4 8.43091 Table 4. Parameter Calibration using a Mean Square Error criteriaTo choose the initial set of parameters we use the Method of Moments . Theo-retical moments are computed from the characteristic function of the model underboth subordinators considered. Matching both, empirical and theoretical moments RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS17 up to order forth leads to the following non-linear system of equations, in the caseof a model under an Inverse Gaussian subordinator:ˆ µ = αµ ∆ t/β ˆ µ = α ∆ t ( µ + β σ + αβµ ∆ t ) /β ˆ µ = αµ ∆ t (3 µ + 3 β σ + 3 αβµ ∆ t + α β µ ∆ t + 3 αβ σ ∆ t ) /β ˆ µ = ( α ∆ t (15 µ + 3 β σ + 18 β µ σ + 15 αβµ ∆ t +6 α β µ ∆ t + α β µ ∆ t + 3 αβ σ ∆ t + 6 α β µ σ ∆ t + 18 αβ µ σ ∆ t )) /β . where ˆ µ k is the empirical k-th moment.The system of equations is solved separately for each regime using the function fsolve in MATLAB based on the trust region algorithm. The results are summa-rized in Table 5.A similar result is obtained for the model under a Gamma subordinator.The method encounter difficulties to find a global minimum in the case where theCommodity ˆ µ ˆ σ ˆ α ˆ β Oil log-return Regime 1 0.1624 0.7213 0.3238 1.6971Oil log-return Regime 2 -0.0354 1.3402 0.0400 1.9584Electricity log-return 1 0.1111 3.1233 28.4386 3.4862Electricity log-return 2 -3.7405 19.5346 0.0132 0.3539 Table 5. Parameter Estimation using Method of Moments underInverse Gaussian Subordinatorempirical moments were calculated using electricity log-return prices. Changing theinitial starting points resulted in varying results, which indicates the presence oflocal minima. Despite these shortcomings the method of moments can be used asan initial solution for a Minimum Distance approach based on the distance betweenthe theoretical and empirical characteristics function of the log-returns, the laterdefined as: ˆ ϕ Z j ∆ t ( u ) = 1 n n (cid:88) k =1 exp(i uz k )for a sample z , z , . . . , z n of n log-returns of the underlying series. See [17] fordetails.The objective function under regime j is defined by: || ϕ Z j ∆ t ( u ; θ ) − ˆ ϕ Z j ∆ t ( u ) || := (cid:32) (cid:90) ∞−∞ | ϕ Z j ∆ t ( u ; θ ) − ˆ ϕ Z j ∆ t ( u ) | w ( x ) dx (cid:33) / , where w is the weight function w ( x ) = (1 / √ π ) exp( − x / θ is the minimum distance estimate of θ if || ϕ Z j ∆ t ( u ; ˆ θ ) − ˆ ϕ Z j ∆ t ( u ) || = inf θ ∈ Θ {|| ϕ Z j ∆ t ( u ; θ ) − ˆ ϕ Z j ∆ t ( u ) || } Again, we apply the algorithm to each regime separately.The integral is computed numerically using a global adaptive quadrature algorithm,where the interval of integration is subdivided and the integration takes place on each subdivided interval. Intervals are further subdivided if the algorithm deter-mines that the integral is not computed to sufficient accuracy. In table 6 estimatesCommodity ˆ µ ˆ σ ˆ α ˆ β Oil log-return Regime 1 0.01736 0.11675 31.648 8.0554Oil log-return Regime 2 -0.4956 2.0078 2.2260 10.141Electricity log-return 1 0.00813 02.0139 67.456 0.00154Electricity log-return 2 5.7435 4.48714 76.004 6.871e-4 Table 6. Parameter Estimation using Minimum Distance Methodunder Inverse Gaussian subordinatorof the model parameters under both regimes and for the two series of underlyingassets are shown.Given a random sample x = ( x , ..., x n ) of a random variable X with an associateddensity function f ( x ; θ ) of the data x under the real world and unknown parameters θ , maximum likelihood estimation (MLE) is a method used to estimate the vectorvalued parameter θ of the model by maximizing the likelihood function: l ( θ ; z ) = n (cid:88) k =1 log f ( z k ; θ ); θ ∈ Θ , with respect to θ . The value of θ is constrained to Θ ⊂ R , the space of allfeasible values of the parameters. The maximum likelihood function L is primarilya function of the unknown parameters θ . The maximum likelihood estimator isgiven by: ˆ θ = arg max θ ∈ Θ l ( θ ; x )In tables 7 and 8 estimations based on the empirical m.l.e. for both regimes andmodels under Gamma and Inverse Gaussian subordinators are shown. In eachCommodity ˆ µ ˆ σ ˆ α ˆ β Oil log-return Regime 1 0.0023 0.0431 42.928 11.9960Oil log-return Regime 2 -0.372 0.52851 17.3008 88.556Electricity log-return 1 5.844e-3 1.5002 93.271 2.1903Electricity log-return 2 -0.0148 7.543 90.5900 0.01770 Table 7. Parameter Estimation using Maximum LikelihoodMethod under Gamma subordinatorCommodity ˆ µ ˆ σ ˆ α ˆ β Oil log-return Regime 1 -0.4883 0.5058 0.64603 63.709Oil log-return Regime 2 0.1201 2.9707 0.00014 9.993Electricity log-return 1 -0.1781 0.25873 0.91878 20.0860Electricity log-return 2 -0.0191 4.9752 5.512e-5 11.016 Table 8. Parameter Estimation using Maximum LikelihoodMethod under Inverse Gaussian subordinator RICING ENERGY CONTRACTS UNDER REGIME SWITCHING TIME-CHANGED MODELS19 case, the values of the volatility σ are found to be higher in the second regime,hence justifying the use of a regime switching model. In nearly every method,the value of | µ | was found to be very small, which is expected as the long termdeterministic contribution to the process is expected to be near zero.In choosing constraints, we set the lower bound of σ, α, β to be some small number (cid:15) = 10 − . We set the upper bound of σ to 5 as the diffusion is expected to besmaller than 1 and for α, β, we set the upper bound to be 100, as the expectedvalue of both Inverse Gaussian and Gamma random variables depends on the ratio α/β rather than any particular value for α and β . The drift µ is expected to besmall, so in most cases, it was constrained to the set [ − , Conclusions We price European-style options with oil and electricity prices as underlyingassets under a switching Levy time-changed noise. These are realistic models thatallow to incorporate stylized features in the dynamic of the prices. Our findingsshow that under this framework a pricing method based on a Fourier-cosine offers anefficient and accurate result when compared with a standard Monte Carlo approach.In addition, we address the problem of parameter estimation and calibration. Tothis end we successfully tried different methods based on both, historic and risk-neutral measure. 6. Acknowledgments This research is supported by the Natural Sciences and Engineering ResearchCouncil of Canada. References [1] Awad H. Al-Mohy and Nicholas J. Higham. A New Scaling and Squaring Algorithm for theMatrix Exponential. SIAM J. Appl. , 31(3):970–989, 2009.[2] O. E. Barndorff-Nielsen. Normal inverse gaussian distributions and stochastic volatility mod-elling. Scandinavian Journal of Statistics , 24:1–14, 1997.[3] Fred Benth, J. Kallsen, and T. Meyer-Brandis. A non-gaussian ornstein-uhlenbeck processfor electricity spot price modeling and derivatives pricing. Applied Mathematical Finance ,14(2):153–169, 2007.[4] Kyriakos Choudarski. Switching l´evy models in continuous time: Finite distributions andoption pricing. Center for Computational Finance and Economic Agents , 2005.[5] Kyriakos Choudarski. Multinomial method for option pricing under variance gamma. Centerfor Applied Mathematics and Economics - ISEG , 2018.[6] F.Fang and C.W. Oosterlee. A novel pricing method for european options based on fourier-cosine series expansions. Journal on Scientific Computing , 2008.[7] MG Figueroa and A Cartea. Pricing in electricity markets: A mean reverting jump diffusionmodel with seasonality. Applied Mathematical Finance , 12(2):313–335, 2005.[8] Elias S.W. Shiu Hans U. Gerber. Option pricing by esscher transforms. Transactions ofSociety of Actuaries, Vol. 46 , 1994.[9] Nathaniel E. Helwig. Density and Distribution Estimation .[10] Nicholas J. Higham. The scaling and squaring method for the matrix exponential revisited. SIAM J. Appl. , 26(4):11791193, 2005.[11] Boyd J.P. Chebyshev and Fourier spectral methods . Springer Verlag, Berlin, 1989.[12] D.B. Madan, P. Carr, and E.C. Chang. The variance gamma process and option pricing. European Finance Review , 79(2), 1998.[13] A.I. McLeod and B. Quenneville. Mean Likelihood Estimators . Statistics and Computing 11,57-65, 2001. [14] Andrea Pascucci. PDE and Martingale Methods in Option Pricing . Bocconi University Press,2011.[15] Stefan Trck Rafa Weron, Michael Bierbrauer. Modeling electricity prices: Jump diffusion andRegime Switching . Physica A, Hugo Steinhaus Center, Wroclaw University of Technology,2004.[16] Kyle Smith. Estimation methods of reference evapotranspiration at unsampled locations. NewMexico Institute of Mining and Technology .[17] Jun Yu. Empirical characteristic function estimation and its applications.