Stability analysis and backward bifurcation on an SEIQR epidemic model with nonlinear innate immunity

1.   Introduction

A weather derivative takes its value from an underlying measure of weather, such as temperature, precipitation, humidity, frost, stream flow or wind speed over a particular period of time, and permits the financial risk associated with climatic conditions to be managed. The first weather derivative transaction appeared over-the-counter (OTC) in US in 1997 and weather derivatives were officially introduced in the major financial markets like Chicago Mercantile Exchange (CME) in 1999. Over the recent years, weather derivatives market has grown rapidly and expanded to Europe and Asia.

Weather derivatives based on temperature are by far the most common type of weather derivatives, accounting for approximately 85% of all weather derivative transactions. Such contracts are typically written on the temperature indices including cumulative average temperature (CAT), heating degree days (HDD) and cooling degree days (CDD). CAT is defined as the sum of the daily average temperature over the period of the contract. HDD or CDD is a measure of the amount by which the daily average temperature during the period deviates from a threshold temperature (typically 18$^{\circ}$C or 65$^{\circ}$F). HDD and CDD are "one-side" metrics and they are critical for creating the risk management tools that utility, agriculture, construction and other firms can use to hedge their activities affected by outdoor climate. For example, a utility company might buy CDD futures if it expects higher temperatures and sell the futures if it expects significantly lower temperatures.

In this study we focus on pricing temperature-based derivatives whose underlying temperature process is modeled as a stochastic process. In the literature, the daily average temperature is often modeled as

where $\Lambda$ is a deterministic function that models the seasonal trend of the temperature process and $Y$, also known as deseasonalized temperature, is a stochastic process to model the random fluctuations from the trend. So far, the mean-reverting OU process is the most popular process used to model $Y$. Queruel (2000) develop a mean-reverting OU model for $Y$, where the diffusion term is a Brownian motion. Benth and Šaltytė-Benth (2005) propose a Lévy-driven OU model with seasonal mean and volatility, where the diffusion term is a Lévy process rather than a Brownian motion. Benth et al. (2007) generalize the OU model to the continuous-time autoregressive (CAR) models with seasonal variance. Swishchuk and Cui (2013) extend the Lévy-driven OU process to the Lévy-driven CAR process for $Y$. Other OU-based models for $Y$ include e.g. fractional OU models of Brody et al. (2002) and Benth (2003), regime switching OU models of Elias et al. (2014) and Evarest et al. (2018).

Temperature-based derivatives pricing has been investigated e.g. in Brody et al. (2002), Benth (2003), Benth and Šaltytė-Benth (2005), Benth et al. (2007), Swishchuk and Cui (2013), Elias et al. (2014) and Evarest et al. (2018). Due to the complex nature of the temperature indices, numerical approximation or Monte Carlo simulations are typically used to obtain the temperature derivative prices. However, in some special cases, the analytical pricing formulas have been derived for some specific temperature derivative products. Benth et al. (2007) obtain the closed-form solutions to CAT future and future option prices when the temperature is modeled by the CAR process driven by Brownian motion. Applying the Fourier technique based on the characteristic function of the temperature process, Swishchuk and Cui (2013) derive the analytical formulas for CAT futures, CAT future options and HDD/CDD futures prices when the underlying temperature is modeled by the general Lévy-driven CAR process. However, for HDD/CDD future options, the analytical formulas for the prices are not available even for OU model driven by Brownian motion; see e.g. Benth et al. (2007) for further discussions.

In this paper, we propose a time-changed OU (TC-OU) model for the deseasonalized temperature process $Y$. We define the time change process $T$ as $T(t) = L(A(t))$, where $L$ is a Lévy subordinator and $A(t) = \int_0^tb(u)du$ with $b(u)$ being a deterministic function of time. Then $Y$ is modeled by the subordinate OU process $Y(t) = X(T(t))$, where $X$ is an OU process driven by a Brownian motion. Our model exhibits several distinct features. First, time changing a continuous time Markov process such as OU process can lead to a much wider class of models than the traditional jump-diffusion models. Depending on whether the Lévy subordinator $L$ has drift or not, the new process will be a jump-diffusion or a pure jump process. Second, for the subordinate process $Y$, its jump measure is state-dependent and mean-reverting, making it ideal for modeling temperature. Third, through the deterministic function $b(u)$, in our model the drift, volatility and jump measure of the process $Y$ are all time dependent and seasonal. Thus, Our model has the potential to capture the salient characteristics of temperature time series. We need to mention that in the time-changed option pricing literature, Carr and Wu (2004) provide option pricing under time-changed Lévy processes, and Cui et al. (2019b) provide a framework for pricing options under general Markov processes, including the subordinate OU process. The subordinate OU process has been successfully applied for callable and putable bonds in Lim et al. (2012), commodities in Li and Linetsky (2014), electricity derivatives in Tong (2017), variance swaps in Tong and Liu (2017), Asian options in Tong and Liu (2018), dual-expiry exotics in Tong et al. (2019) and autocallables in Tong (2019). We also refer to Linetsky and Mitchell (2008) for a brief introduction of subordinate Markov processes and Mendoza-Arriaga et al. (2010), Arriaga and Linetsky (2013), Li et al. (2016b) and Li et al. (2017) for their further applications.

To price temperature derivatives, we assume that under the pricing measure, our model remains in the same form as in the physical measure, but with a different parameter set to account for risk premia. We are able to derive the analytical pricing formulas for CAT, HDD and CDD future and future option prices using eigenfunction expansion method. We need to emphasize that eigenfunction expansion method is particularly suitable for pricing contingent claims written on the subordinate processes. The subordinate process is as analytically tractable as the original process without time change. The subordinate process shares the same eigenfunctions with the original one and the only modification is the replacement of eigenvalues of the original process with the Laplace transform of the time change process. We refer to Linetsky and Mitchell (2008) for the surveys on the eigenfunction expansion method and Linetsky (2004), Mendoza-Arriaga et al. (2010), Lim et al. (2012), Mendoza-Arriaga and Linetsky (2013), Li and Linetsky (2014), Li et al. (2017), Tong (2017), Tong and Liu (2017), Tong and Liu (2018), Tong et al. (2019) and Tong (2019) for its various applications.

In this paper, we apply the eigenfunction expansion method for the pricing of temperature derivatives under the time-changed OU process. There are alternative analytical pricing methodologies that are applicable to the general case when the underlying stochastic process is a time-changed Markov process. A promising method is to approximate the underlying process by a continuous time Markov chain (CTMC). We refer to Cui et al. (2019a) for a review of this recent and growing literature and Li and Zhang (2018) for the relationship between CTMC approximation and the eigenfunction expansion approach applied in this work.

The rest of the paper is organized as follows. In Section 2, we introduce the general framework for modeling temperature as a time-changed OU process, where the time change process is modeled by a Lévy subordinator time changed by a deterministic clock with seasonable activity rate. In Section 3, we introduce the eigenfunction expansion method for our new model and also discuss how to calculate some important integrals that are essential for the determination of the eigenfunction expansion coefficients. In Section 4, we apply the eigenfunction expansion method to obtain the analytical pricing formulas for CAT, CDD and HDD futures and future options. In Section 5, we estimate our new model and compare its performance with the competing models using the temperature data. We also numerically study the impact of market price of risk on the temperature derivatives.

2.   The model set-up

Let $(\Omega, \mathcal{F}, P)$ be a probability space with an information filtration $(\mathcal{F}_t)$. Suppose under the physical measure $P$, the temperature process $S$ is modeled as

where $\Lambda(t)$ is the seasonal trend term that is modeled by a combination of linear and cosine function:

Similar forms can be found in Benth and Šaltytė-Benth (2005) and Swishchuk and Cui (2013).

The deseasonalized temperature process $Y$ is specified as a time-changed OU process, that is

where $T$ is a time change process and $X$ is an OU process

where $\kappa > 0$ and $B(t)$ is a standard Brownian motion.

We choose the time change process $T$ to be of the form $T(t) = L(A(t))$, where $A(t) = \int_0^tb(u)du$ and the activity rate $b(u)$ is modeled as a positive periodical function

where $b_0\ge|b_1|$. Then $A(t)$ can be solved as

We remark that other specifications for $b(u)$ are possible, as long as $A(t)$ is easy to compute.

To introduce jumps into the temperature dynamics, we model the process $L$ as a Lévy subordinator independent of $B$ in (5). The Lévy subordinator $L$ is a nondecreasing process with positive jumps and non-negative drift with the Laplace transform:

where $\phi$ is the Lévy exponent and given by the Lévy-Khintchine formula (see e.g. Sato, 1999)

where $\gamma\ge 0$ and the Lévy measure $\nu$ must satisfy

An important sub-class of Lévy subordinators are the tempered stable subordinators. For such subordinators, the Lévy measure $\nu(s)$ is given by

where $C > 0$ and $\eta > 0$. Important special cases are the Gamma subordinator with $p = 0$, the IG subordinator with $p = \frac{1}{2}$ and the compound Poisson subordinator with $p = -1$. For such subordinators, the Lévy exponent is given by

We can also reparameterize the exponent by setting

where $\vartheta = {\rm{E}}(L(1))-\gamma$ and $\omega = {\rm{Var}}(L(1))$.

Then, the reparameterized exponent takes the following form:

From $T(t) = L(A(t))$, we know $T$ is an additive subordinator, i.e. a non-negative and non-decreasing additive process (see e.g. Li et al., 2016b) and $Y$ is an additive subordinate OU process (ASub-OU), whose drift, diffusion volatility and jumps are all time-dependent and seasonal. For more detailed discussions on general additive subordinators, we refer to Li et al. (2016a) and Li et al. (2017). For our model, we can calculate the Laplace transform of $T$ as

where $\psi(\lambda, u)$ is called the density of the Laplace exponent and is given by

We note that $T(t) = L(t)$ when $b_0 = 1$ and $b_1 = 0$. In this case, $Y$ is a Lévy subordinate OU process (LSub-OU) and it will be a jump-diffusion process with mean-reverting diffusion drift and mean-reverting jumps if $\gamma > 0$ or a pure jump process with mean-reverting jumps if $\gamma = 0$. For this specification, the drift, diffusion volatility and jumps intensity are all time-independent.

3.   Eigenfunction expansion method for the TC-OU model

To price temperature derivatives, we need to specify a risk-neutral probability measure $Q$. However, it is not possible to store or trade temperature and therefore, the temperature market is an incomplete market. Any probability measure $Q$ equivalent to the physical measure $P$ is a risk-neutral measure. In this paper, we will follow Li et al. (2016a) and Li et al. (2016b) and assume that the process for $Y$ under measure $Q$ remains the same as under $P$, but with possibly different parameters to account for risk premium. The generating tuples ($\kappa$, $\theta$, $\sigma$, $b(\cdot)$, $\gamma$, $\nu(\cdot)$) for $Y$ under the measure $P$ then becomes (${\kappa}$, $\bar{\theta}$, ${\sigma}$, ${b}(\cdot)$, ${\gamma}$, ${\nu}(\cdot)$) under the measure $Q$. Therefore, we assume after the change of measure, $\theta$ becomes $\bar{\theta}$ and other parameters remain the same.

Under the risk-neutral measure $Q$, the OU process $X$ defined in (5) becomes

where $B^Q$ is a standard Brownian motion, independent of the time change process $T$.

For the OU process $X(t)$ defined in (17), its infinitesimal generator $\mathcal{L}$ is given by

where $f$ is a transformation function. $f'$ and $f''$ are the first- and the second-order derivatives of $f$, respectively.

The infinitesimal generator $\mathcal{L}$ with domain $dom(\mathcal{L})$ is always self-adjoint on the Hilbert space $L^2((-\infty, \infty), m)$ of functions on $(-\infty, \infty)$ square integrable with the speed measure $m(dx) = m(x)dx$ and endowed with the inner product $(f, g)$, ¹ where

¹An operator $\mathcal{L}$ with domain $dom(\mathcal{L})$ is said to be self-adjoint if $(\mathcal{L}f, g) = (f, \mathcal{L}g)$ $\forall f, g\in dom(\mathcal{L})$.

and

We are interested in the transition semigroup $(\mathcal{P}_t)_{t\ge 0}$ of $X$ defined by

We can use the eigenfunction expansion method for self-adjoint operators in Hilbert spaces to write down the spectral decomposition of $\mathcal{P}_tf(x)$ (see e.g. Linetsky and Mitchell, 2008).

Proposition 1. For any $f\in L^2((-\infty, \infty), m)$, the transition semigroup of the OU process $X$ has an eigenfunction expansion of the form

where $f_n = (f, \varphi_n)$, $\{\lambda_n\}$ are the eigenvalues of $-\mathcal{L}$ and $\{\varphi_n\}$ are the corresponding eigenfunctions satisfying the following Sturm-Liouville equation

The eigenvalues and eigenfunctions of OU process can be summarized in the following result (see e.g. Linetsky and Mitchell, 2008):

Proposition 2. For the OU process $X$ defined in (17), the eigenvalues $\lambda_n$ and the eigenfunction $\varphi_n$, $n = 0, 1, \ldots$, are

and

where $\xi = \frac{\sqrt{\kappa}}{\sigma}(x-\bar{\theta})$, $N_n = \sqrt{\frac{\sigma\sqrt{\kappa}}{2^{n+1}n!\sqrt{\pi}}}$ and $H_n(x)$ is the Hermite polynomial defined as

For the subordinate OU process $Y$ defined in (4), we can also compute the transition semigroup $(\mathcal{P}_t^{\psi})_{t\ge 0}$ of $Y$, where $\psi$ is the density of Laplace exponent for the additive subordinator $T$:

We can employ the eigenfunction expansion method again to compute the semigroup $\mathcal{P}_t^{\psi}f(x)$ using the following result (see e.g. Li and Linetsky, 2014):

Proposition 3. For any $f\in L^2((-\infty, \infty), m)$, the transition semigroup of the subordinate OU process $Y$ has an eigenfunction expansion of the form

where $\{\lambda_n\}$ and $\{\varphi_n\}$ are the eigenvalues and the eigenfunctions of OU process and can be obtained from (24) and (25), respectively. The exponential function $\exp\left(-\int_0^t\psi(\lambda_n, u)du\right)$ is the Laplace transform of the time change process $T$ and can be calculated using (15).

From the above results, it is clear that a key feature of the eigenfunction expansion method is that the temporal and spatial variables are separated. The time variable $t$ enters the expansion only through the exponential function $\exp(-\lambda_n t)$. The eigenfunction expansion of time-changed process $Y$ has the same form as the original process $X$, but with $\exp(-\lambda_n t)$ replaced by $\exp\left(-\int_0^t\psi(\lambda_n, u)du\right)$. As long as the Laplace transform of the time change process $T$ can be calculated in closed form, the time-changed model will be as tractable as the original model.

For our new pricing model, the eigenvalues and the eigenfunctions of the OU process can be calculated easily from Proposition 2. To employ the eigenfunction expansion method to calculate the temperature derivative prices, we still need to obtain the eigenfunction expansion coefficient $f_n$. Here, we provide the formulas for several integrals that will be useful for calculating $f_n$. The following integrals can be computed using the results of Prudnikov et al. (1986) and Lim et al. (2012):

Lemma 1.

(i)

(ii)

where $\Phi(x)$ is the CDF of standard normal distribution.

(iii)

(iv) Define

Then, $\gamma_{m, n}(u)$ can be calculated recursively as follows:

and for $m\neq n$,

We can also prove the following results:

Corollary 1.

(i)

(ii)

(iii)

(iv)

(v)

(vi)

Proof. (ⅰ)-(ⅵ) follow directly from Lemma 1.

4.   Valuation of temperature derivatives

The payoff of a temperature derivative is typically linked to an underlying temperature index such as CAT, HDD or CDD, which can be defined as

and

where $c$ is typically 18$^{\circ}$C or 65$^{\circ}$F.

Assume the time $t$ value of temperature $S(t) = x$. Denote $Ind\in\{CAT, CDD, HDD\}$ and let $F_{Ind}(t, t_{N_1}, t_{N_2}, x)$ with $t < t_{N_1} < t_{N_2}$ be the time $t$ price of temperature futures written on the corresponding index over period $[t_{N_1}, t_{N_2}]$. The future price can be calculated as:

We have the following parity for the temperature futures.

Lemma 2. The prices of CAT, HDD and CDD futures have the following relationship

Proof. The result follows from $\max(x-c, 0)-\max(c-x, 0) = x-c$.

Assume the time $t$ value of temperature $S(t) = x$. Let $P_{Ind}(t, K, t_{N_0}, t_{N_1}, t_{N_2}, x)$ be the time $t$ price of a put option with strike price $K$ and maturity $t_{N_0}$, written on a future contract with payoff $F_{Ind}(t_{N_0}, t_{N_1}, t_{N_2}, S(t_{N_0}))$. Then the put price can be calculated from:

where $r$ is the risk-free interest rate.

Let $C_{Ind}(t, K, t_{N_0}, t_{N_1}, t_{N_2}, x)$ be the corresponding temperature call option. We have the following put-call parity for the temperature future options.

Lemma 3.

Proof. The result follows from $\max(x-K, 0)-\max(K-x, 0) = x-K$.

We can derive the analytical solutions for CAT, HDD and CDD future prices from the following results.

Proposition 4. Assume that the stochastic processes for the asset price $S(t)$ is specified in (2), (4), (17) and the time change process $T(t)$ is an additive subordinator with the density of Laplace transform $\psi$, then the time $t$ prices of CAT, HDD and CDD futures can be calculated as

(i)

where $\lambda_n$ and $\varphi_n$ are in (24) and (25), respectively. Furthermore, $f_n^{CAT} = \bar{b}_n$.

(ii)

where $f_{n, k}^{HDD} = (c-\Lambda(t_k))a_n(c-\Lambda(t_k))-b_n(c-\Lambda(t_k))$.

(iii)

where $f_{n, k}^{CDD} = (\Lambda(t_k)-c)(\bar{a}_n-a_n(c-\Lambda(t_k)))+\bar{b}_n-b_n(c-\Lambda(t_k))$.

Proof.

We will focus on (ⅱ). (ⅰ) and (ⅲ) can be proved in a similar way.

Using eigenfunction expansion and Corollary 1, we have

where

We can also derive the analytical solutions for CAT, HDD and CDD future option prices from the following results.

Proposition 5. Assume that the stochastic processes for the asset price $S(t)$ is specified in (2), (4), (17) and the time change process $T(t)$ is an additive subordinator with the density of Laplace transform $\psi$, then the time $t$ prices of CAT, HDD and CDD future options can be calculated as

(i)

where

where $x^*$ = $x^{**}-\Lambda(t_{N_0})$ and $x^{**}$ is the solution to the equation $F_{CAT}(t_{N_0}, t_{N_1}, t_{N_2}, x) = K$.

(ii)

where

where $x^*$ = $x^{**}-\Lambda(t_{N_0})$ and $x^{**}$ is the solution to the equation $F_{HDD}(t_{N_0}, t_{N_1}, t_{N_2}, x) = K$.

(iii)

where

where $x^*$ = $x^{**}-\Lambda(t_{N_0})$ and $x^{**}$ is the solution to the equation $F_{CDD}(t_{N_0}, t_{N_1}, t_{N_2}, x) = K$.

Proof. To prove (ⅰ), let $x^{**}$ be the solution to the equation $F_{CAT}(t_{N_0}, t_{N_1}, t_{N_2}, x) = K$ and $x^*$ = $x^{**}-\Lambda(t_{N_0})$, then we have

Using eigenfunction expansion, Corollary 1 and Proposition 4.1, we get

where

Similarly, we can prove (ⅱ) and (ⅲ).

5.   Data analysis and numerical study

We carry out an empirical study using the average daily temperature of Toronto, Canada, as measured at Pearson International Airport from January 1, 2003 to December 31, 2012. The data was collected from the website of National Climate Data and Information Archive of Canada. In Figure 1, we plot the average temperature for the whole sample. It is clear that the temperature has a seasonal pattern. As specified in Section 2, we model the temperature $S(t)$ as the sum of the seasonal trend term $\Lambda(t)$ and the deseasonalized stochastic component $Y(t)$. The estimation procedure is performed in two steps.

In the first step, using linear regression we estimate $\Lambda$ according to (3). The estimation results are in Table 1. All the parameters in (3) are statistically significant at the $1\%$ level. In particular, the temperature has a weak but significantly positive linear trend. In Figure 1, we also plot the seasonal component $\Lambda$, which indicates that it captures the seasonal pattern in the temperature very well.

In the second step, we estimate three models for the deseasonalized temperature $Y$. The first model is the OU model based on the process in (5) and $Y(t) = X(t)$. We impose the restriction $\theta = 0$ for identification purpose. The other two models are the LSub-OU and the ASub-OU models. For both models, the Lévy subordinator is taken to be the Gamma subordinator with Lévy exponent given in (14). To uniquely identify the parameters of the models, we set $\vartheta = 1$. For the ASub-OU model, the activity rate function $b$ is specified in (6), with the restrictions $b_0 = 1$ and $b_0\ge |b_1|$. Note if $b_1 = 0$, the ASub-OU model reduces to the LSub-OU model. We also set the drift term of the Gamma subordinator $\gamma = 0$, so $Y$ will be a pure jump process with mean reverting jumps.

We choose maximum likelihood estimation (MLE) method for all the models. For a sample of size $n$ of $\{S(t_i)\}_{i = 1}^n$, the log conditional likelihood function is given by

where $\Theta$ is the set of parameters to be estimated. $p(t_i, t_{i+1}, S(t_i), S(t_{i+1}); \Theta)$ is the transition density function.

Define

and

Then, the transition density function for the OU model is

Let $g_t(du)$ be the transition probability distribution of a Gamma subordinator with zero drift, mean $\vartheta$ and variance $\omega$, then $g_t(du)$ is given by

Then, the transition density function for the ASub-OU model is given by

This integral can be efficiently computed by the Gauss-Laguerre quadrature. Note if we impose the restriction $b_1 = 0$, we can also obtain the transition density function for the LSub-OU model.

The MLE results for the deseasonalized temperature process $Y$ are provided in Table 2. We find that all of the parameters of each model are estimated with high significance. Once the time change is introduced, the volatility parameter of the OU process drops significantly. This is not surprising since in the models with time change, the volatility is not only achieved through the diffusion term of diffusion process, but also can be realized through the contribution of the jumps. In Table 2, we also report the log-likelihood function values for three models. Since OU is nested in LSub-OU and LSub-OU is nested in ASub-OU, we can perform the likelihood ratio test to compare the goodness-of-fit. The likelihood ratio test statistics are 119.8 and 117.44 for OU versus LSub-OU and LSub-OU versus ASub-OU, respectively. Using the critical value of 3.84 for $95\%$ quantile of the Chi-square distribution with one degree of freedom, we can clearly see that the two time-changed models are superior to the OU model and the ASub-OU model outperforms the LSub-OU model.

In Figure 2, we plot the deseasonalized temperature $Y$ together with the activity rate function $b$. It is clear that the seasonality still exists in the process $Y$. The diffusive volatility and jump intensity display seasonality, which are well captured by the function $b$ in the ASub-OU model.

We also apply the eigenfunction expansion method to compute the temperature future and future option prices using the results in Section 4. In order to apply the eigenfunction expansion method to compute the derivative prices, we need to truncate eigenfunctions expansions after a finite number of terms. We price CAT and CDD futures and put options on the futures using the analytical formulas derived in Section 4. We take the model parameters from those reported in Tables 1 and 2. Other parameters are assumed to be the same as those in Figures 3 and 4. In Table 3, we demonstrate that only a small number of eigenfunctions are needed for the prices to converge.

To take the risk premium into consideration, we need to specify the values of $\bar{\theta}$ for OU process $X$ under the risk-neutral Measure $Q$. In Figures 3 and 4, we compare the futures and future options for different values of $\bar{\theta}$ under the OU, LSub-OU and ASub-OU models. As expected, both CAT and CDD future prices are a linear function of $\bar{\theta}$, whereas CAT and CDD future option prices are a decreasing function of $\bar{\theta}$. It is also clear that the market price of risk has a significant role for temperature derivatives. We also find different models can generate different values for temperature derivative prices. For CAT futures and future options, three models produce similar prices. However, the model prices from ASub-OU model can deviate from the other models significantly for CDD futures and future options.

6.   Conclusion

This paper proposes a new model for temperature based on the time-changed OU process, where the time change process is a Lévy subordinator time changed by a deterministic clock with seasonal activity rate. Compared to the previous models, our model allows for state-dependent jumps. Furthermore, under our model, the drift, diffusion volatility and jumps are all seasonal which are empirically supported in the temperature data. Another advantage of our model is that we are able to employ the eigenfunction expansion technique to derive the analytical pricing formulas for the temperature futures and future options. From an empirical study using Toronto temperature data, we show that clearly our model with additive subordination provides a better fit than the competing models. We also demonstrate that the market price of risk has significant impact on the temperature derivative prices.

There are several future works that can be considered. We can introduce a stochastic volatility component into our model by a further time change based on the absolutely continuous time change processes. The combination of a subordinator and an integral of some positive functions of a Markov process with analytically tractable Laplace transform has been studied in Mendoza-Arriaga et al. (2010) and Li and Linetsky (2014). The model based on the combined time change will still be tractable since it still possesses the closed-form Laplace transform for the time change, which enables us to continue to apply the eigenfunction expansion method. In addition, we would like to calibrate our model using the real market data for the temperature derivatives and compare its performance with the existing popular models.

Acknowledgments

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflict of interest

All authors declare no conflicts of interest in this paper.