Pricing Temperature Derivatives under a Time-Changed Levy Model
PPRICING TEMPERATURE DERIVATIVES UNDER ATIME-CHANGED LEVY MODEL
PABLO OLIVARES, RYERSON UNIVERSITY
Abstract.
The objective of the paper is to price weather contractsusing temperature as the underlying process when the later follows amean-reverting dynamics driven by a time-changed Brownian motioncoupled to a Gamma Levy subordinator and time-dependent determin-istic volatility. This type of model captures the complexity of the tem-perature dynamic providing a more accurate valuation of their associateweather contracts. An approximated price is obtained by a Fourier ex-pansion of its characteristic function combined with a selection of theequivalent martingale measure following the Esscher transform proposedin Gerber and Shiu (1994). Introduction:
The objective of the paper is to price weather contracts using temperatureas the underlying process when the later follows a mean-reverting dynamicsdriven by a time-changed Brownian motion coupled to a Gamma Levy sub-ordinator and a time-dependent volatility function. The process reverts to aseasonal periodic deterministic process, while the volatility is considered alsoa periodic function of time, see Dacunha-Castelle, Hoang and Parey (2015)for the later. Temperature models driven by Levy noises and stochasticvolatility have been originally considered in Benth and Benth-S(2009).This type of model captures the complexity of the temperature dynamicproviding a more accurate valuation of their associate weather contracts.On the other hand, the availability of an explicit analytical expression of thecharacteristic function of the process allows for its Fourier expansion withrespect of its characteristic function, which in turn leads to compute the ap-proximated price under an equivalent martingale measure (EMM) obtainedfrom the Esscher transform, see Gerber and Shiu (1994).The combination of these three elements, namely the model, the pricingmethod and the choice of the EMM in the context of weather derivativesoffers a novel methodology for pricing such contracts.Methods based on Fourier expansions of the characteristic function in oneand two dimensions are implement in Fang and Oosterlee (2008)to Europeancontracts and further extended to other derivatives by the same authors, seeFang and Oosterlee (2014).
Key words and phrases.
Temperatures, weather contracts, Fourier expansions, Time-changed Levy subordinators . a r X i v : . [ q -f i n . P R ] M a y Finally, we fit the model to a series of daily average temperatures at Pearsonairport, Ontario, Canada during the period 2014-2019.The organization of the paper is the following:In section 2 we describe the main model for the temperature process andobtain the characteristic function associated with it. In section 3 we discussthe implementation of the Fourier expansion techniques, while in section 4we show the numerical results in the fitting of the model, pricing results andtheir sensitivities to key parameters.2.
Modeling temperature
Let (Ω , A , ( F t ) t ≥ , P ) be a filtered probability space verifying the usualconditions. For a stochastic process ( X t ) t ≥ defined on the space filteredspace above the functions ϕ V t and l X ( u ) = t log ϕ V t ( − iu ) defines its char-acteristic function and the cumulat generating function respectively. Whenthe process has stationary and independent increments the later does not de-pend on t . The σ -algebra F Y t = σ ( Y u , ≤ u ≤ t ) is the σ -algebra generatedby the random variables Y u , ≤ u ≤ t . The changes of the temperature overan interval [ t, t + h ) are denoted ∆ T t = T t + h − T t . For a process ( X t ) t ≥ ,the discounted process ( ˜ X t ) t ≥ is defined as ˜ X t = e − rt X t , where r is thecontstant interest rate.Let ( T t ) t ≥ be the daily average temperature process defined on the filteredspace above. The average temperature is taken as the arithmetic mean be-tween the maximum and the minimum temperature during a given day.We assume the temperature process ( T t ) t ≥ verifies the stochastic differentialequation:(1) dT t = α ( s t − T t ) dt + σ t dV t where ( s t ) t ≥ is a deterministic seasonal process such that: s t = β + β t + β sin (cid:18) π t (cid:19) + β cos (cid:18) π t (cid:19) (2)The parameter α is the mean-reversion rate to the seasonal component. Thebackground noise ( V t ) t ≥ will be specified later on.The solution of equation (1) is given in the following lemma. Lemma 1.
The solution of equation (1) is: T t = e − αt T + αK ( t, α ) + W t (3) with W t = (cid:82) t σ u e − α ( t − u ) dV u and K ( t, α ) = (cid:90) t s u e − α ( t − u ) du = 1 α (1 − e − αt ) β + 1 α (1 − α )(1 − e − αt ) β + 1 α [cos( π t ) − e − αt − α π sin( π t )]1 − α ( π ) β + 1 α [sin( π t ) + α π ( e − αt − cos( π t )]1 + α ( π ) β Proof.
We apply Ito formula to the function f ( x, y ) = xe αy and the process( T t , t ).Hence: T t e αt = T + (cid:90) t e αu dT u − + α (cid:90) t e αu T u − du + (cid:88) u ≤ t [ T u e αu − T u − e αu − ∆ T u − e αu ]= T + α (cid:90) t ( s u − T u − ) e αu du + (cid:90) t σ u e αu dV u + α (cid:90) t T u − e αu du = T + αe αt K ( t, α ) + (cid:90) t σ u e αu dV u Multiplying by e − αt on both sides leads to equation (3). (cid:3) We assume the volatility also follows a deterministic seasonal componentprocess:(4) σ t = c + c t + c sin (cid:18) π t (cid:19) + c cos (cid:18) π t (cid:19) where c j ≥ , j = 0 , , , . We will need to compute the characteristic function of some integrals of thebackground noise process. To this end we will make use of a well-knownresult about functional of a Levy process ( ξ t ) t ≥ and a measurable function f :(5) E ( exp ( i (cid:90) t f ( s ) dξ s )) = exp ( (cid:90) t l ξ ( − if ( s )) ds )In order to select the EMM for pricing purposes we take an Esscher transformof the historic measure P . See Gerber and Shiu(1994) for a rationale in termsof a utility-maximization criteria.For a stochastic process ( X t ) t ≥ we consider its Esscher transform:(6) d Q θt dP t = exp( θX t − tl 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 and l θX ( u ) respectively the characteristic function and moment generating function of a process ( X t ) t ≥ under the probability Q θ obtained by an Esscher transformation as given in equation (6).For consistency we denote ϕ X t := ϕ X t and l X = l V .By analogy with the case of financial underlying assets the risk marketpremium measure Q θ making the discounted temperatures process ( ˜ T t ) t ≥ amartingale for r > Q θ is denoted E θ .We set a subordinator process ( R t ) t ≥ and the time-changed process ( V t ) t ≥ verifying: V t = B R t + µ R t (7)Here µ ∈ R is a parameters in the model and ( B t ) t ≥ is a standard Brownianmotion.The following result describes the characteristic function of the temperatureprocess under the historic measure P . Proposition 2.
Under the model described by equations (1), (2) and (7)the characteristic function of T t under the probability P is: ϕ T t ( u ) = C ( t, α ) exp( (cid:90) t l R ( − iuµ σ s e − α ( t − s ) − u σ s e − α ( t − s ) ) ds (8) where: C ( t, α ) = exp ( iue − αtT + αK ( t, α )) Proof.
By conditioning: ϕ V t ( u ) = E [ E [exp( i ( uV t ) /R t )]] = E [exp( iuµ R t ) E [exp( iuB R t /R t )]]= E [exp( iuµ R t ) exp( − R t u )] = E [exp( i ( uµ + 12 iu ) R t )]= ϕ R t ( uµ + 12 iu )Hence:(9) l V ( u ) = l R ( uµ + 12 u )By lemma 1 and formula (5): ϕ T t ( u ) = E [ e iuT t ] = C ( t, α ) E [ exp ( iu (cid:90) t σ s e − α ( t − s ) dV s )]= C ( t, α ) exp ( (cid:90) t l V ( − iuσ s e − α ( t − s ) ) ds )(10)Combined with equation (9), equation (8) immediately follows. (cid:3) The results below provides the characteristic function of the temperatureprocess under the EMM defined via an Esscher transform.
Proposition 3.
Let ( T t ) t ≥ be the temperature process defined by equations(1)-(7). Then, the characteristic function under the Esscher EMM Q θ is: ϕ θT t ( u ) = C ( t, α ) C ( t, θ ) I t ( u, θ )(11) where C ( t, α ) is defined as in the previous proposition and: C ( t, θ ) = exp( − tl R ( θµ + 12 θ )) I t ( u, θ ) = exp( (cid:90) t l R ( − iuµ σ s e − α ( t − s ) + µ θ + 12 ( − iuσ s e α ( t − s ) + θ ) ) ds ) and for any T > the parameter θ verifies: (12) l (cid:48) V ( θ ) = − e ( α + r ) T (1 − ˜ C ( T, α )) K − ( α, T ) where: K ( α, T ) = (cid:90) T σ u e αu du = c α ( e αT −
1) + c Tα e αT − c α ( e αT − − π c (cid:18) cos( 2 π T ) − (cid:19) + 3652 π c (cid:18) sin( 2 π T ) − (cid:19) Proof.
Notice that: ϕ θV t ( u ) = E θ ( e iuV t e θV t − tl V ( θ ) ) = ϕ V t ( u − iθ ) ϕ V t ( − iθ )and l θV ( u ) = l V ( u + θ ) − l V ( θ ).Then, similarly to proposition 2: ϕ θT t ( u ) = C ( t, α ) exp( (cid:90) t l θV ( − iuσ s e − α ( t − s ) ) ds )= C ( t, α ) exp( − tl V ( θ )) exp( (cid:90) t l V ( − iuσ s e − α ( t − s ) + θ ) ds )from which equation (11) follows.By equation (6) the discounted temperature process ( ˜ T t ) t ≥ verifies:˜ T t = ˜ C ( t, α ) + ˜ W t It is a Q θ -martingale if and only if for any 0 ≤ s < t : E θ ( ˜ T t / F s ) = ˜ T s ⇔ E θ ( ˜ W t − ˜ W s / F s ) = ˜ C ( s, α ) − ˜ C ( t, α ) But: E θ ( ˜ W t − ˜ W s / F s ) = E θ ( e − ( α + r ) t (cid:90) t σ u e αu dV u − e − ( α + r ) s (cid:90) s σ u e αu dV u / F s )= E θ ( e − ( α + r ) t (cid:90) ts σ u e αu dV u + ( e − ( α + r ) t − e − ( α + r ) s ) (cid:90) s σ u e αu dV u / F s )= E θ ( e − ( α + r ) t (cid:90) ts σ u e αu dV u )+ ( e − ( α + r ) t − e − ( α + r ) s ) (cid:90) s σ u e αu dV u On the other hand, from equation (5): ϕ θV t ( x ) = exp ( (cid:90) t l θV ( ixσ u e αu ) du ) = exp ( (cid:90) t ( l V ( ixσ u e αu + θ ) − l V ( θ )) du )(13)Hence: E θ ( e − ( α + r ) t (cid:90) ts σ u e αu dV u ) = e − ( α + r ) t i ( ϕ θV t ) (cid:48) ( x ) | x =0 = − e − ( α + r ) t i ( i (cid:90) ts σ u e αu l (cid:48) V ( − ixσ u e αu + θ ) du | x =0 exp( (cid:90) ts ( l V ( − ixσ u e αu + θ ) − l V ( θ )) du ) | x =0 = − e − ( α + r ) t l (cid:48) V ( θ ) (cid:90) ts σ u e αu du In particular for t = T and u = 0 we have the result in equation (12), thatfollows from elementary calculation. (cid:3) Remark 4.
Notice that the characteristic function under the probability P is obtained from equation (11) taking θ = 0 . Hence we write I t ( u ) = I t ( u, , ϕ Y t = ϕ Y t and Q = P . Example 5.
Gamma subordinatorConsider the subordinator ( R t ) t ≥ is a Gamma process with parameters a > , b > , see Carr and Madan (1999), with respective characteristic functionand Laplace exponent: ϕ R t ( u ) = (cid:18) − iub (cid:19) − at , a > , b > l R ( u ) = − a log (cid:16) − ub (cid:17) , u < b Therefore: ϕ V t ( u ) = ϕ R t ( µ u + 12 iu ) = (cid:32) − i ( µ u + iu ) b (cid:33) − at = (cid:18) − iµ ub + 12 b u (cid:19) − at l V ( u ) = − a log A ( u ) where: A ( u ) = 1 − µ ub − b u Moreover: l θV ( u ) = l V ( u + θ ) − l V ( θ )= − a [log A ( u + θ ) − log A ( θ )]= − a log (cid:18) A ( u + θ ) A ( θ ) (cid:19) = − a log (cid:32) − µ ub − b ( u + 2 θu )1 − µ θb − b θ (cid:33) To compute the characteristic function of the temperature T t under the EMMEsscher transformation given by equation (11) we have: C ( t, α ) = exp ( iue − αtT + αK ( t, α )) C ( t, θ ) = exp( − tl R ( θµ + 12 θ )) = A at ( θ ) I t ( u, θ ) = exp( (cid:90) t l θV ( − iuσ s e − α ( t − s ) ) ds )= exp (cid:32) − a (cid:90) t log (cid:32) A ( − iuσ s e − α ( t − s ) + θ ) A ( θ ) (cid:33) ds (cid:33) To compute the Gerber-Shiu parameter, from the martingale condition givenby equation (12): l (cid:48) V ( θ ) = a ( µ + θ ) bA ( θ ) = − e ( α + r ) T (1 − ˜ C ( T, α )) K ( α, T ) − )= − e ( α + r ) T (1 − e − rT A aT ( θ ) K ( α, T ) − )= − e ( α + r ) T + e αT A aT ( θ ) K ( α, T ) − Therefore, the value θ ∗ that solves: (14) µ + 12 θ + ba e ( α + r ) T A ( θ ) − ba e αT A aT +11 ( θ ) K ( α, T ) − = 0 makes the discounted prices martingales under the Esscher transformation. Pricing weather options
Weather contracts are based on cumulate temperatures (CAT), heating-degrees-days (HDD) or cooling-degrees-days (CDD) over certain period [0 , T ].Futures and option contracts are offered in Chicago Mercantile Exchange.They are respectively defined as: ξ T = CAT = T (cid:88) k =1 T k ξ ,T = HDD = T (cid:88) k =1 ( c − T k ) + ξ ,T = CDD = T (cid:88) k =1 ( T k − c ) + The typical case is c = 18 o Celsius.For concreteness we focus on a CAT index. To this end for convenience werewrite the CAT index as: ξ T = T (cid:88) k =1 T k = T (cid:88) t =1 ( T + t (cid:88) j =1 ∆ T j )= T T + T (cid:88) j =1 γ j ∆ T j (15)where the changes in temperature ∆ T j = T j +1 − T j are independent randomvariables and γ j = T − j + 1.A general payoff of the temperature weather derivative, consisting in a com-bination of a European long put and a long call with different strikes, knownas strangle , is given by:(16) h ( ξ T ) = d ( ξ T − K ) + + d ( K − ξ T ) + , d j > , K > K > , j = 1 , d and d are the costs per unit of temperature below (resp. above)the threshold K (resp. K ) known as tick sizes .The price of a temperature contract over the period [0 , T ] is : p W = d e − r T E Q ( ξ T − K ) + + d e − r T E Q ( K − ξ T ) + = d e − r T (cid:90) R ( x − K ) + f ξ T ( x, θ ) dx + d e − r T (cid:90) R ( K − x ) + f ξ T ( x, θ ) dx (17)where r is the interest rate and f ξ T ( x, θ ) is the p.d.f. of the cumulatedtemperature under the EMM measure. A Fourier expansion of the p.d.f. f ξ T ( x, θ ) on an interval [ b , b ] is given by: f ξ ( x, θ ) = + ∞ (cid:88) k =0 A k ( θ ) cos (cid:18) kπ x − b b − b (cid:19) (18)where the coefficients in the expansion, the first of them divided by two, are: A k ( θ ) = 2 b − b (cid:90) b b f ξ T ( y, θ ) cos (cid:16) kπ y − b b − b (cid:17) dy (cid:39) b − b (cid:90) b b f ξ T ( y, θ ) Re (cid:16) e ikπ y − b b − b (cid:17) dy = 2 b − b Re (cid:16) (cid:90) b b f ξ T ( y, θ ) e ikπ y − b b − b dy (cid:17) = 2 b − b exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) (19)Replacing (19) into (18), then (18) in (17) we have: (cid:90) R ( x − K ) + f ξ T ( x, θ ) dx (cid:39) + ∞ (cid:88) k =0 A k ( θ ) (cid:90) b b ( x − K ) + cos (cid:18) kπ x − b b − b (cid:19) dx (cid:39) N (cid:88) k =0 A k ( θ ) (cid:90) b b ( x − K ) cos (cid:18) kπ x − b b − b (cid:19) dx = N (cid:88) k =0 A k ( θ ) (cid:90) b b xcos (cid:18) kπ x − b b − b (cid:19) dx − K (cid:90) b b cos (cid:18) kπ x − b b − b (cid:19) dx = 2 b − b N (cid:88) k =0 exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) (cid:90) b b xcos (cid:18) kπ x − b b − b (cid:19) dx − K b − b N (cid:88) k =0 exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) (cid:90) b b cos (cid:18) kπ x − b b − b (cid:19) dx = 2 b − b N (cid:88) k =0 exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) (cid:90) b b xcos (cid:18) kπ x − b b − b (cid:19) dx − K b − b N (cid:88) k =0 exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) (cid:90) b b cos (cid:18) kπ x − b b − b (cid:19) dx where b = max ( b , K ) < b and from equation (15): ϕ θξ T ( u ) = e iT T T (cid:89) j =1 ϕ θ ∆ T j ( γ j u ) Moreover, for k > (cid:90) b b xcos (cid:18) kπ x − b b − b (cid:19) dx = ( b − b ) b kπ sin (cid:18) kπ b − b b − b (cid:19) + (cid:18) ( b − b ) kπ (cid:19) (cid:18) ( − k − cos (cid:18) kπ b − b b − b (cid:19)(cid:19)(cid:90) b b cos (cid:18) kπ x − b b − b (cid:19) dx = b − b kπ sin (cid:18) kπ b − b b − b (cid:19) Then, separating the first term in the summation: (cid:90) R ( x − K ) + f ξ T ( x, θ ) dx (cid:39) ( b − b ) b − b )+ 2 b − b N (cid:88) k =1 exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17)(cid:32) ( b − b ) b kπ sin (cid:18) kπ b − b b − b (cid:19) + (cid:18) ( b − b ) kπ (cid:19) (cid:18) ( − k − cos (cid:18) kπ b − b b − b (cid:19)(cid:19)(cid:33) + 2 K N (cid:88) k =1 kπ exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) sin (cid:18) kπ b − b b − b (cid:19) In a similar analysis: (cid:90) R ( K − x ) + f ξ T ( x, θ ) dx (cid:39) K b − b − ( b − b ) − K N (cid:88) k =1 kπ exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17) sin (cid:18) kπ b − b b − b (cid:19) − b N (cid:88) k =1 kπ exp (cid:18) − i kπb b − b (cid:19) ϕ θξ T (cid:16) kπb − b (cid:17)(cid:32) ( b − b ) b kπ sin (cid:18) kπ b − b b − b (cid:19) + (cid:18) ( b − b ) kπ (cid:19) (cid:18) − cos (cid:18) kπ b − b b − b (cid:19)(cid:19)(cid:33) where b = min ( b , K ). The delicate choice of the truncation values b and b as well as the number of terms in the truncated expansion depends on themodel considered, it is discussed in Fang and Oosterlee (2008). For detailederror analysis of the truncation and numerical errors present in the FourierCosine method we refer the reader to the work of Fang and Oosterlee (2008).We address this issue in the next section related to numerical aspects of themethod. Numerical results
We divide the section into three parts. In the first one we do a descriptivestatistical analysis and fit the seasonal component. In the second we discussthe parameter estimation, while in the final part we implement the pricingmethod outlined above and analyze its sensitivities with respect to modeland contract parameters. Partial results in subsections 4.1 and ?? have beenpreviously considered in Porthiyas (2019).4.1. Statistical analysis and parameter estimation.
Daily tempera-ture data (in degree Celsius) at Toronto from January 1st, 2013 to Novem-ber 15th, 2018 have been collected from Environment and Climate Change,Canada. The data is gathered from the Pearson International Airportweather station and yield 2145 data points. Observations consist of anaverage between the daily maximum and minimum temperatures. Missingobservations are replaced by a seven-day moving average around the missingpoint.
Figure 1.
Historic daily average temperature of Torontofrom 1/1/2013 to 15/11/2018A preliminary statistical analysis of the temperature data shows the de-scriptive statistics as in Table 1. As can be seen, the skewness of the datais negative indicating a longer tail to the left. The kurtosis is less than 3indicating more frequent but modest movements of temperature than wouldbe expected under assumptions of normal distribution.It can be observed from both the histogram and the kernel density esti-mate in Figure 2 that the temperature data is bimodal. The left peak iscentered around the mean temperature in winter and right peak is centered Table 1.
Summary of the seriesMean Minimum Maximum Std Dev Skewness Kurtosis9.0483 -22.30 30.45 11.0593 -0.3021 2.1481
Table 2.
Kolmogorov-Smirnov test resultsp-value KSSTAT Critical value0 0.6889 0.0292around the mean temperature in summer.
Figure 2.
Histogram and Kernel density estimate of dailyaverage temperatures in Toronto from 1/1/2013 to15/11/2018Table 2 shows the results of a KolmogorovSmirnov test. This is a goodness-of-fit test to verify whether the data is from a normal distribution. It canbe concluded from the p-value of zero and a KSSTAT value significantlygreater than the critical value, that the temperature data do not seem tofollow a normal distribution.The seasonal component as described in equation (2) is adjusted via aregression model. The results are shown in table 3.As it can be seen from Table 3, the slope term b in the regression fitis small but significantly different for zero, which indicates the existence ofa linear trend in temperature rising, consistent with other climatic studiessignaling the past decade as the warmest one since temperature is recorded. Estimate SE t-Stat Conf. int. pValue b b b -5.8796 0.14176 -41.476 (-6.143, -5.590) 3.103e-276 b -12.866 0.14287 -90.052 (-13.13, -12.57) 0 Table 3.Figure 3.
Seasonal trend for Toronto daily mean temperatureIt must be noted in those cases, a larger set of temperature data for 40 yearsor more was used.4.2.
Parameter estimation.
We base our analysis on the log-return seriesgiven by:(20) X j ∆ = log (cid:18) T ( j +1)∆ T j ∆ (cid:19) = Y ( j +1)∆ − Y j ∆ , j = 1 , , . . . , n where ∆ > Figure 4.
Simulated temperature graph for 2018 is shown,compared with the actual observations
Figure 5.
Simulated trajectories for different values of themean-reverting level5.
Acknowledgments
The author would like to thank the Natural Sciences and EngineeringResearch Council of Canada for its support. Conclusions
A mean-reverting time-changed Levy process with periodic mean-revertinglevel and volatility offers a fair model for temperatures at Pearson Interna-tional Airport temperatures.On the other hand, pricing methods based on Fourier expansions provide analternative algorithm under the models and the underlying series considered.Weather temperature prices are efficiently computed on a PC in reasonabletime.
References [1] Benth, F. and Sgarra, C. (2009)The Risk Premium and the Esscher Transform inPower Markets. Stochastic Analysis and Applications, vol.30(1).[2] F. E. Benth and J. Saltyte-Benth. Stochastic modelling of temperature variationswith a view towards weather derivatives. Applied Mathematical Finance, 12(1):53-85, 2005.[3] Carr P. and Madan D.(1999)Option valuation using the fast Fourier transform. Jour-nal of Computational Finance, vol.2, no.4, pg.61-73[4] Fang Fang and Cornelis W. Oosterlee (2008). A novel pricing method for Europeanoptions based on Fourier-cosine series expansion.
SIAM Journal of Scientific Com-puting .[5] Gerber, H. U. and Shiu, E. S. W. (1994) Option pricing by Esscher-transforms. Trans-actions of the Society of Actuaries 46, 99191.[6] Kallsen Jan and Muhle-Karbe Johannes(2011). Method of moment estimation intime-changed Lvy models, Statistics and Risk Modeling, De Gruyter, vol. 28(2),pages 169-194, May.[7] Porthiyas, J. (2019) Pricing Weather Contracts Under Mean-reverting Levy Pro-cesses. M Sc. Thesis. Department of Mathematics. Ryerson University.[8] Chunfa Wang (2017). Pricing European Options by Stable Fourier-Cosine Series Ex-pansions https://arxiv.org/abs/1701.00886[9] Samuel Asante Gyamerah, Philip Ngare, and Dennis Ikpe (2018)Regime-Switching Temperature Dynamics Model for Weather Derivatives. InternationalJournal of Stochastic Analysis. Volume 2018, Article ID 8534131, 15 pages.https://doi.org/10.1155/2018/8534131[10] A. Swishchuk and K. Cui (2013) Weather derivatives with applications to canadiandata. Journal of Mathematical Finance, 3(1):81-95, 2013.[11] Bowen Zhang and Cornelis W. Oosterlee. Efficient pricing of european-style asianoptoins under exponential