A Stochastic Harmonic Oscillator Temperature Model for the Valuation of Weather Derivatives

: Stochastic processes are employed in this paper to capture the evolution of daily mean temperatures, with the goal of pricing temperature-based weather options. A stochastic harmonic oscillator model is proposed for the temperature dynamics and results of numerical simulations and parameter estimation are presented. The temperature model is used to price a one-month call option and a sensitivity analysis is undertaken to examine how call option prices are affected when the model parameters are varied.


Introduction
Weather influences our daily lives and affects the volumes of sales and everyday trading in a wide range of businesses. For instance, the powerful typhoons that hit Tokyo and nearby regions in the fall of 2019 caused millions of dollars in damage to infrastructures and extensive damage to agriculture, forestry and fisheries industries [1,2]. Meanwhile, on the other side of the Pacific, the historic cold wave and winter storm of 2021 in Northwest, Central and Eastern America caused an estimated loss of $20.4 billion and represents the most costly U.S. winter storm event on record surpassing (nearly doubling the inflationadjusted cost) the 93 Superstorm [3]. These two examples are just amongst the many other illustrations of how adverse weather could continuously affect the economy and routinely pose significant financial risks to 70% of all economic sectors worldwide [4].
In the United States, studies show a 1.2% decline in the annual gross domestic product (GDP) for every 1 • C increase in temperature by the end of this century [5]. While at a global scale, according to the International Monetary Fund, failure to meet a temperature rise to 1.5 • C may curb the world's real GDP per capita by roughly 7% by 2100 [6]. Both demonstrate the grave repercussions of Mother Nature that could potentially affect already stressed economies.
To hedge weather-related financial risks, financial instruments called 'weather derivatives' were introduced in the financial markets in 1997 in the form of futures, options on futures and swaps [7,8]. The value of such a derivative is based on some underlying weather-dependent variable such as wind speed, amount of snowfall or rainfall and temperature levels. Financial products have been developed for industries facing business uncertainty as a result of weather fluctuations, allowing companies to trade and hedge those risks caused by weather volatility.
In these exotic weather derivative contracts, the payoffs depend on the chosen weather index, reflecting the weather conditions against which protection is needed. The largest portion of traded weather derivatives are written on temperature-based indices. Clearly, the temperature does not have a price in itself and therefore is not a tradeable asset. As a consequence, following the standard derivative pricing theoretical framework, the valuation of these derivatives relies heavily on the appropriate modelling of the index. This context motivates the purpose of this work, that is, to propose an easy-to-calibrate stochastic temperature model that does not require specific meteorological background and inputs. Hence, the model can be straightforwardly implemented by practitioners in the financial industry.
We follow a stochastic modelling approach in this paper. Previous results in this direction are due to Benth and Saltyte-Benth [14], Dornier and Querel [15], Brody et al. [16] and Alaton et al. [17]. All of these papers put forward a mean reverting Ornstein-Uhlenbeck process to describe the temperature fluctuations, some of which are driven by fractional Brownian motion or generalised Lévy processes instead of standard Brownian motion.
It can be argued whether in the aforementioned results it is sufficient to describe the temperature dynamics by a single stochastic process. Since observed temperature data indicate a seasonality trend, we will deviate from the previous approaches and consider as our starting point the classical harmonic oscillator model [18], which is described by a second-order ordinary differential equation (ODE). We will convert this model to a system of two coupled stochastic differential equations (SDEs) with two driving Wiener processes. The system will turn out to be still analytically tractable and we fit our stochastic harmonic oscillator model to actual temperature data. Using a risk-neutral framework, we will price weather derivatives, more specifically, a call option on an underlying variable which is a temperature-based index.
The outline of this paper is as follows: in Section 2, we derive a stochastic temperature model based on the classical harmonic oscillator and compare it to a popular stochastic temperature model. Section 3 presents numerical simulations of the proposed model, as well as results of parameter estimation using theoretical and real temperature data. In Section 4, we use the temperature model to price a one-month call option and perform a sensitivity analysis to examine how call option prices are affected when the model parameters are varied. We provide brief concluding remarks in Section 5.

Derivation of a Stochastic Harmonic Oscillator Temperature Model
Let T(t) ∈ R denote the daily average temperature at time t ≥ 0. We begin by considering a general (deterministic) sinusoidal function of time, namely where A, B, C, D ∈ R and A, B = 0. It can be shown that T(t) satisfies the ODE where k = B and a = B 2 D. Without loss of generality, we may take k > 0 and a ∈ R. Let x 1 (t) = T(t) and x 2 (t) = T (t). Then, (1) can be written as a system of two first-order ODEs which in differential form is equivalent to We introduce randomness in (2) as follows. Consider the following system for the stochastic processes X 1 = {X 1 (t) : t ≥ 0} and X 2 = {X 2 (t) : t ≥ 0}: where σ 1 , σ 2 ≥ 0 are parameters that control the size of the fluctuations and W 1 = {W 1 (t) : t ≥ 0} and W 2 = {W 2 (t) : t ≥ 0} are (possibly correlated) standard Wiener processes. Equation (3) can be expressed in matrix form as where From Evans [19], we can write the solution of (3) as where e tA is the matrix exponential of A and X 0 = X(0). The eigenvalues of A are easily seen to be ±ik. Recalling the explicit formula for the matrix exponential [20], we obtain e tA = cos(kt) 1 k sin(kt) −k sin(kt) cos(kt) and therefore e tA B = a k sin(kt) a cos(kt) , e tA C = σ 1 cos(kt) σ 2 k sin(kt) −σ 1 k sin(kt) σ 2 cos(kt) .
Since in the deterministic model (2) we have x 1 (t) = T(t), for the stochastic model (3), we associate X 1 (t) = T(t). Therefore, taking the first component in (4) and setting θ 1 = X 1 (0) and θ 2 = X 2 (0), we get Thus, our stochastic harmonic oscillator temperature model is By renaming the parameters in (5), we can write We wish to estimate the parameters k and c j for j = 0, 1, 2, 3, 4 by calibrating (6) to temperature data.
Remark 1. Dornier and Querel [15] and Alaton et al. [17] assumed that the daily average temperature is modelled by The term S(t) is the deterministic seasonal component of the form where a, b, c j and d j for j = 1, 2, 3 are constants to be determined. The stochastic deseasonalised component is described by the stochastic process and W = {W(t) : t ≥ 0} is a standard Wiener process. The parameters λ, µ and σ denote the mean reversion rate, long term average and volatility, respectively, of the process X. It is well known that X is an Ornstein-Uhlenbeck process and Thus, the expected value of T(t) in (7) is Note that this is not periodic in t. Comparing (7) with (5), we see that in (5) there are two driving Wiener processes W 1 and W 2 . Moreover, the expected value of T(t) in (5) is which is a periodic function of t and shows that the deterministic component of the temperature would also satisfy a harmonic oscillator ODE. The model we propose here is different from other two-factor models in the literature as (5) is obtained from the system (3), which in turn is a stochastic analogue of the harmonic oscillator model (1). In contrast, Hell et al. [21] established a general theoretical framework for the pricing of temperature derivatives by consistent factor models. Groll et al. [22] provided an empirical back-up for the theoretical framework introduced by Hell et al. [21] and put their pricing approach into practice.

Numerical Simulations of a Stochastic Harmonic Oscillator Temperature Model
In this section, we show the results of various numerical simulations of the models (5) and (3). All simulations are done in Scilab (see https://www.scilab.org, accessed on 12 October 2021), a free and open-source, cross-platform programming language.

Simulation of a Stochastic Integral as a Stochastic Process
If f is a deterministic function of t and u, where t, u ∈ [0, t max ], and W = {W(t) : 0 ≤ t ≤ t max } is a Wiener process, then we use the Riemann-Stieltjes approximation As For example, suppose that we take n = 400 and t max = 4. Then, sample paths of and are shown in Figures 1 and 2, respectively. We observe that the sample paths of Y 1 are more 'jagged' than those of Y 2 .

Simulation of a Stochastic Harmonic Oscillator Temperature Model
To simulate the stochastic process (5), take n = 1000 and t max = 4. Then, choose We generate two correlated Wiener processes and ρ ∈ [−1, 1] is the correlation coefficient. Choose ρ = 0.5 for definiteness. A sample path for the process T is represented by the red graph in Figure 3. Next, we solve the system (3) using the Euler-Maruyama algorithm. We remark that the Milstein and Euler-Maruyama schemes are equivalent for this model. Define the functions Setting t j = (j − 1)∆t for j = 1, 2, . . . , n, we see that X 1 (t 1 ) = X 1 (0) = θ 1 and X 2 (t 1 ) = X 2 (0) = θ 2 . For j = 1, 2, . . . , n − 1, we compute recursively Plotting the first component of the approximate solution yields the blue graph in Figure 3. As expected, the plots are relatively close to each other since we identify X 1 with T. The L 2 -norm of the difference between the sample paths of X 1 and T is 16.77. Choosing a larger n may decrease the error at the cost of a slower computational speed.

Parameter Estimation for a Stochastic Harmonic Oscillator Temperature Model Using Theoretical Data
The next step is to fit (5) using theoretical data. As before, let n = 1000 and t max = 4, while assuming the parameters in (11). We use (5) to generate a set of theoretical data {T j : j = 1, 2, . . . , 1000}, where T j = T(t j ) = T((j − 1)∆t) for j = 1, 2, . . . , 1000. The profile of the theoretical temperature curve is given by the red curve in Figure 4, which is actually hidden behind the blue curve. Then, with the exception of k = 6 in (11), we 'pretend' that all of the parameters a, θ 1 , θ 2 , σ 1 and σ 2 are unknown and we attempt to fit (5) to the theoretically generated data.  Since the data set {(t j , T j ) : j = 1, 2, . . . , n} is known and the Wiener processes W 1 and W 2 can be generated, and we assume for now that k is given, then the stochastic integrals t j 0 cos(k(t j − u)) dW 1 (u), t j 0 sin(k(t j − u)) dW 2 (u) appearing in (6) are also known for each j = 1, 2, . . . , n because they can be computed as in Section 3.1. It is therefore convenient to use (6) since this equation is linear in the parameters c 0 , c 1 , c 2 , c 3 and c 4 if we further replace T(t) by T j , cos(kt) by cos(kt j ) and sin(kt) by sin(kt j ). Then, a standard multi-dimensional linear regression can be implemented to estimate these parameters. That is, if Y = T 1 · · · T n * (where (·) * denotes the matrix transpose) and 1 cos(kt 1 ) sin(kt 1 ) t 1 0 cos(k(t 1 − u)) dW 1 (u) t 1 0 sin(k(t 1 − u)) dW 2 (u) 1 cos(kt 2 ) sin(kt 2 ) t 2 0 cos(k(t 2 − u)) dW 1 (u) The result isĉ 0 = 0.06,ĉ 1 = −0.06,ĉ 2 = −0.17,ĉ 3 = 2.50,ĉ 4 = 0.17.
Using these estimated parameters in (5) generates a set {T j : j = 1, 2, . . . , 1000} and is represented by the blue curve in Figure 4. As the original parameters in (11) have been recovered, this explains why the blue curve lies on top of the red curve. The L 2 -norm of the difference between {T j : j = 1, 2, . . . , 1000} and {T j : j = 1, 2, . . . , 1000} is 0.00. For comparison, the green curve in Figure 4 is obtained from (5) but without the stochastic integral terms, i.e., (8).

Remark 2.
Earlier, we assumed that k was known. Now suppose that k is unknown. Assuming for the moment that the stochastic integral terms in (5) are negligible, the 'period' P can be estimated through P = 2π k or k = 2π P .
In Figure 3, we see that there is a 'crest' that starts at around t = 2.60 and ends at around t = 3.60. Taking P 3.60 − 2.60 = 1.00, then k 6.28. This could be taken as the k fixed in the above linear regression. However, a more reasonable approach is to assume that k lies within some range, say 5 ≤ k ≤ 7, and test different values of k ranging from k = 5 to k = 7. For each such fixed k, we apply linear regression as above and compute the L 2 -error between {T j : j = 1, 2, . . . , 1000} and {T j : j = 1, 2, . . . , 1000}. When the L 2 -error is less than some predetermined tolerance, the corresponding k is chosen and the rest of the parameters a, θ 1 , θ 2 , σ 1 and σ 2 are estimated.

Parameter Estimation for a Stochastic Harmonic Oscillator Temperature Model Using Real Data
In this section, we use a data set of daily mean temperatures observed at the Toronto Pearson International Airport located in Mississauga, Ontario, Canada and compiled by the National Climatic Data Center (NCDC) (see http://www.ncdc.noaa.gov/cdo-web/ datasets, accessed on 12 October2021) of the United States Department of Commerce. The period considered is from January 2005 to December 2008, for a total of 1461 data points. Each daily mean temperature is obtained by taking the average of the lowest and highest temperature readings for that day and converting from • F to • C.
Note that t j will be expressed in years and t n = (n − 1)∆t = n 365 . In this case, n = 1461 and t max = t n = t 1461 = 4. Fix k = 6.3 (an explanation for this choice will be provided in Remark 3 below). Following the procedure described in Section 3.3, the estimated parameters arê a = 358.55,θ 1 = −4.14,θ 2 = −33.69,σ 1 = 0.45,σ 2 = 1.64. (12) In Figure 5, the red graph is the observed temperature data {T j : j = 1, 2, . . . , 1461} while the blue curve (behind the green curve) is the fitted temperature data {T j : j = 1, 2, . . . , 1461}, i.e., (5) with the parameters replaced by the estimated parameter values in (12). The green curve is (5) but without the stochastic integral terms, which is equivalent to {E(T j ) : j = 1, 2, . . . , 1461} from (8).  Unlike in Figure 4, the fitted temperature data in Figure 5 is not too 'jagged' but is actually close to E(T(t)). This can be seen from the fact that 0 <σ 1 <σ 2 . Hence, in (5), the contribution from the stochastic integral with the cos integrand is not as pronounced as that from the stochastic integral with the sin integrand. As pointed out previously for Figures 1 and 2, the stochastic integral with the sin integrand tends to produce less 'jagged' sample paths. Remark 3. The L 2 -norm of the difference between {T j : j = 1, 2, . . . , 1461} and {T j : j = 1, 2, . . . , 1461} is 167.20. Following the explanation given in Remark 2, the value k = 6.3 was chosen because it produced the smallest L 2 -error amongst all the values of k taken from k = 5 to k = 7 in increments of ∆k = 0.1.

Basics of a Weather Derivative Contract
A standard weather derivative contract typically has the following elements [8]: • contract period (traded contracts are seasonal and monthly); • weather-based variable (contracts are based on temperature, snowfall, rainfall, etc.); • measurement station (e.g., La Guardia Airport, New York City); • index (e.g., HDD or CDD); • payoff function (makes it possible to convert the index to cash flow); • tick (multiplier for each weather index point).
The underlying indices of these contracts, from which these financial instruments derive their value, represent the biggest difference with respect to other financial derivatives. There are two major types of indices used for temperature. The first is the cumulative average temperature (CAT), which is the sum of the daily average temperature values throughout the contract period. The other major type is an index based on degree days. Heating degree days (HDD) is a measurement of the amount the daily temperature falls below a set baseline B (typically, B = 18 • C = 65 • F), accumulated over the contract period. Cooling degree days (CDD) is similar to HDD but measures the amount the daily temperature exceeds B.
Suppose that n is the number of days over which the HDD/CDD index is considered. Let T j denote the average of the maximum and minimum temperatures observed on day j, where j = 1, 2, . . . , n. Then, These indices were originally created by the energy industry as a proxy for the domestic demand for heating and cooling. The motivation for the above formulas is that, when the temperature is above (respectively, below) the baseline, people would consume energy on cooling (respectively, heating) appliances but not on heating (respectively, cooling) appliances.
At the Chicago Mercantile Exchange (CME), HDD contracts are traded from January to April and from October to December, while CDD contracts are traded from April to October. Months in which both indices are traded are known as shoulder months.

Remark 4.
As HDD/CDD indices depend on the number of days over which they are considered, and max(x, 0) ≥ 0 for any x ∈ R, we deduce that Similarly, CDD(n + 1) ≥ CDD(n). Hence, HDD/CDD indices rise with increasing n.
For example, let us return to the results in Section 3.4 which are summarised in Figure 5. HDD months for Toronto at the CME are from January to April and from October to December.  Figure 6. The data set {T j : 1 ≤ j ≤ 1461} was used to compute the observed HDD (red graph), while the data set {T j : 1 ≤ j ≤ 1461} was used to compute the expected HDD (blue graph). As anticipated, HDD(n + 1) ≥ HDD(n).

Oct 05
Apr 06 Oct 06 Apr 07 Oct 07 Apr 08 2000 4000 6000 8000 10000 12000 14000 No. of days HDD Observed HDD Expected HDD Figure 6. Plots of the observed HDD and expected HDD using (13) and the data from Section 3.4, respectively.

Call Option Pricing: An Example
In this section, we illustrate how to use our stochastic harmonic oscillator temperature model to price a call option, where the underlying variable is the HDD index. Other types of options or trading strategies may also be considered.
Suppose that the contract period is from 1-31 December (n = 31), the strike HDD is K = 550, the tick size is s = $20 and the risk-free rate is r = 0.0206. The latter was determined from the 10-year Canadian government bond yields compiled by Bloomberg (see http://www.bloomberg.com/quote/GCAN10YR:IND, accessed on 12 October 2021).
The stochastic harmonic oscillator temperature model we use iŝ T j =â 6.3 2 + θ 1 −â 6.3 2 cos(6.3t j ) +θ 2 6.3 sin(6.3t j ) where the parametersâ,θ 1 ,θ 2 ,σ 1 andσ 2 are as in (12) Using a risk-neutral framework, the price C of a call option is therefore To evaluate the expectation in (16), we take the average of 5000 sample paths of the process (14) and substitute it into (15). The resulting call price is C = $3196.84. Figure 7 shows the variation of the call price over a range of strike HDDs. The figure is generated by varying K in (16) from 500 to 800.

Sensitivity Analysis
Following Siu et al. [23], we perform a sensitivity analysis for the call option price when the parameters of the model vary. The aim is to identify the key parameters that have significant impact on the call price.
To explain the methodology, let us suppose that we wish to investigate the sensitivity of C with respect to a by looking at how C behaves with respect to K as a is varied. More precisely, we assume thatθ 1 ,θ 2 ,σ 1 andσ 2 are fixed as in (12) and letâ vary through five different values on an interval of unit radius centred at the value 358.88 ofâ in (12). For eachâ in this interval, we perform the calculations similar to those for Figure 7. Figure 8 shows the impact of varying a, Figure 9 of varying θ 1 , Figure 10 of varying θ 2 , Figure 11 of varying σ 1 and Figure 12 of varying θ 2 . We conclude that the only parameter that has a significant effect on the call option price is θ 1 = T(0). This result reveals that, under this stochastic harmonic oscillator temperature model, we must monitor the actual state of the temperature dynamics because of its considerable effect on the value of the weather-based option.

Conclusions
In this article, we proposed a new stochastic model for the temperature dynamics based on the classical harmonic oscillator (strictly speaking, a driven harmonic oscillator because of the term a in (1)). The temperature model is derived from a system of linear SDEs with two driving Wiener processes. We benchmarked our model by first using theoretical temperature data and then using actual temperature data observed at the Toronto Pearson International Airport. Substituting our temperature model in a temperature-based index such as HDD, we priced a call option under a risk-neutral framework and performed a sensitivity analysis. Numerical simulations suggest that the parameter related to the temperature at time zero was the only significant one that affected the call price as the strike HDD was varied. Of course, other temperature-based indices such as CDD or CAT can also be used, as well as other trading strategies besides a call option can also be considered following a similar procedure.
An extension of this model that is being currently investigated is a stochastic harmonic oscillator model with damping. The motivation behind this is the consideration of a temperature model that incorporates the gradual increase in temperature (due to the exponential factor introduced by damping) brought about by global warming, for example. Another possible extension is the modification of the stochastic terms in (3), e.g., from σ j dW j (t) to σ j X j (t) dW j (t) for j = 1, 2. The resulting system remains linear and therefore a similar analysis can be performed.