A stochastic volatility model for the valuation of temperature derivatives

Abstract

Accepted by: Konstantinos Nikolopoulos

This paper develops a new stochastic volatility model for the average daily temperature. It is a natural extension of a Gaussian model in which the temperature returns to a seasonal trend with a deterministic time-dependent volatility. The new model allows to be more conservative regarding extreme events while keeping tractability. We give a method based on conditional least squares to estimate the parameters on daily data and estimate our model on eight major European cities. We then show how to calculate efficiently the average payoff of weather derivatives both by Monte-Carlo and Fourier transform techniques. This new model allows to better assess the risk related to temperature volatility.

1. Introduction

With the increased awareness on climate risk and its tremendous consequences on all economic sectors, the demand for financial tools that enable to hedge this weather-related perils has significantly increased in the last decades. The development of weather derivatives coincides with this trend. First launched in 1997 with over the counter (OTC) contracts, the Chicago Mercantile Exchange introduced standardized contracts for American cities in 1999. Until date, the open market remains relatively small and therefore lacks of liquidity. However, OTC market has developed between industrialists and large insurance and finance companies leading to the necessity to develop satisfactory risk valuation methodologies for those derivatives. These derivatives usually bring on a Heating Degree Day (HDD) index or a Cooling Degree Day (CDD) index or a Cumulative Average Temperature (CAT) index. Since most of the deals are OTC, the market is incomplete and thus no arbitrage free pricing. Weather derivative sellers are more interested in understanding the distribution of the payoff, and how this distribution may change under stressed conditions. These information allow them to determine their pricing depending on their risk appetite as well as to assess their maximum losses.

The literature on the valuation of weather derivatives mostly relies on the development of temperature models. On the one hand, there are continuous time models. Brody et al. (2002) suggest to use an Ornstein–Uhlenbeck process driven by a fractional Brownian motion with periodic time-dependent parameters. Benth & Šaltytė Benth (2007) consider the same dynamics with a classical Brownian motion with a further expansion in periodic functions and including a trend for the temperature. Šaltytė Benth & Benth (2012) propose a continuous time version of autoregressive models. More recently, Groll et al. (2016) have developed a continuous-time model with factors for the temperature forecasting curve. All these papers use their model to calculate the average payoff of standard weather derivatives on HDD, CDD or CAT. On the other hand, Tol (1996) and Franses et al. (2001) have proposed discrete-time GARCH models for the temperature. Cao & Wei (2004) and Campbell & Diebold (2005) use this approach in view of pricing derivatives. Šaltytė Benth & Benth (2012) complete this latter study and show that a low order of autoregression is enough to fit well temperature data. Recently, Meng & Taylor (2022) have proposed an extension of this family of models that gives a joint modeling of the daily minimum and maximum temperature.

In this paper, we focus on the temperature-related derivatives which have been, by far, the most studied tools in the literature. Following the studies of Benth et al. (2007) and many others on temperature dynamics, we develop a time inhomogeneous affine stochastic volatility model inspired from the celebrated Heston (1993) model for equity. This model integrates additional flexibility to the temperature dynamics enabling a better time-continuous modelling of the temperature while keeping tractability. In particular, it is more conservative regarding extreme events than the corresponding Gaussian model. Following the recent work of Bolyog & Pap (2019), we develop a conditional least squares estimation method for its parameters, which is easy to implement. Besides, we propose two different pricing algorithms. The first one based on simulation enables to sample the payoff distribution and then to compute empirically quantities such as the average payoff or quantiles of this distribution. The second one takes advantage of the affine structure and is an adaptation of the Fast Fourier Transform method introduced by Carr & Madan (1999) to calculate the average payoff of some weather derivatives. Besides, this second method, combined with control variates, enables to reduce considerably the computational time (up to |$10^{5}$|⁠). We can then easily calculate the sensitivity of the average payoff with respect to the different parameters, including those on the temperature volatility. Thus, an important contribution of our model is to better assess the risk behind the volatility of the temperature. Finally, we compare the results of this approach to price weather derivatives with common business practices.

This paper can also be leveraged by practitioners. First, our model contributes to a deeper understanding of daily temperature models. This topic is particularly sensitive in the context of global warming and increased extreme events. This paper reviews different approaches to trend modeling and volatility capture. Second, the study describes the step-by-step pricing of temperature derivatives from calibration to computational implementation. This is particularly useful for reproductive purposes and has already been implemented, with some adjustments, in the industrial framework. Finally, the last section provides a comparison with standard business pricing practices to bridge the gap between research and application.

The paper is organized as follows. The first section presents our model for the average daily temperature. In particular, we explain the rationale behind this model and how it goes beyond the dynamics studied in the literature up to date. The second section focuses on the estimation of the different parameters of our model. We develop a conditional least squares approach and check the robustness of the fitting on simulated data. The third section concentrates on the pricing of weather derivatives. It develops both Monte Carlo and Fast Fourier Transform (FFT) algorithms and studies the sensitivity to model parameters of the average payoff. It also compares these pricing methods with current business approaches.

2. Temperature models

2.1 A stochastic volatility model for temperature dynamics

Temperature dynamics have largely been analyzed in the statistics literature and they are of particular interest in the field of weather derivatives. These analysis apply to average daily temperature models that are defined as the average of maximum and minimum daily temperature, i.e. |$T = \frac{T_{max} + T_{min}}{2}$|⁠. This choice comes from the insurance contracts that typically use this average for the daily temperature.

Different models have been suggested for the associated |$(T_{t})_{t \geq 0}$| process. In the present paper, we investigate the interest of applying a stochastic volatility model for the temperature process |$(T_{t})_{t \geq 0}$|⁠. Our model extends the existing models proposed in the literature, while it enables to capture important fluctuations that were not illustrated by former models. Namely, we introduce Model (M)

where |$(W_{t})_{t \geq 0}$| and |$(Z_{t})_{t \geq 0}$| are independent Brownian motions, |$\kappa , \eta ,K>0$|⁠, |$\rho \in [-1,1]$|⁠, |$\sigma ^{2}$| is a non-negative function and the functions |$s$| and |$\sigma ^{2}$| will be taken as in (2) and (5). We will denote |$(\mathcal{F}_{t})_{t\ge 0}$| the filtration generated by |$(W,Z)$|⁠, so that the processes |$T$| and |$\zeta $| are adapted to it. This model integrates the following three components.

• First, the function |$s$| represents the trend and seasonality of the average temperature |$(T_{t})_{t \geq 0}$|⁠. This function is deterministic, bounded and continuously differentiable. In this paper, we will consider the following parametric form:

  $$ \begin{align*} &s(t)=\alpha_0 +\beta_0 t+\alpha_1 \sin\left( \frac{2\pi}{365}t \right)+\beta_1 \cos\left( \frac{2\pi}{365}t \right),\ t\ge 0, \end{align*} $$

  see Subsection 2.2.1 for a discussion on this.

• Second, the detrended and deseasonalized temperature process |$(\tilde{T}_{t})_{t \geq 0}$| follows a mean-reverting process which enables to include memory effects into the model. The parameter |$\kappa $| tunes the mean-reversion speed.

• Third, the volatility of the |$(\tilde{T}_{t})_{t \geq 0}$| process denoted |$(\zeta _{t})_{t \geq 0}$| follows a Cox–Ingersoll–Ross (CIR) process (Cox et al., 1985) with time-dependent parameters (Hull & White, 2015). The process |$(\zeta _{t})_{t \geq 0}$| includes a seasonal deterministic component |$\sigma ^{2}(\cdot )$| which is supposed deterministic, bounded, continuously differentiable and non-negative, so that the process |$\zeta $| is well defined and remains non-negative if |$\zeta _{0}\ge 0$|⁠. In this paper, we will the following parametric form for |$\sigma ^{2}$|⁠:

  $$ \begin{align*} & \sigma^2(t)=\gamma_0+ \sum_{k=1}^2 \gamma_k \sin\left( k\frac{2\pi}{365}t \right)+\delta_k \cos\left(k \frac{2\pi}{365}t \right).\end{align*} $$

  The parameter |$K$| tunes the mean-reversion of the volatility process and |$\eta $| corresponds to the volatility of the volatility.

Here and through the paper, we note |$f(t)$| (time |$t\ge 0$| in parenthesis) a deterministic function and |$F_{t}$| (time in index) a stochastic process.

This model is somehow an adaptation of the celebrated Heston (1993) model. It enables to go beyond Ornstein–Uhlenbeck models (Franses et al., 2001, Taylor & Buizza, 2004, Benth & Šaltytė Benth, 2007, Benth et al., 2007) and is an alternative to GARCH volatility models (Franses et al., 2001, Taylor & Buizza, 2004, Campbell & Diebold, 2005, Šaltytė Benth & Benth, 2012) while keeping some flexibility as we will show in the next sections.

In the following sections, we study daily average temperature data series for 8 major European cities: Stockholm, Paris, Amsterdam, Berlin, Brussels, London, Rome and Madrid. The data series present daily data from 1st January 1980 to 31st December 2020. After removing 29 February of leap years, this gives a time series of 14 965 observations coming from weather stations. Table A1 in  Appendix A summarizes the characteristics of the weather stations. These time series were extracted from Speedwell, the main historical and settlement data provider for weather derivatives. The company performs quality checks including physical consistency comparison, statistical consistency tests and comparison with neighboring sources and corrects the data if necessary. For instance, gap filling is mentioned for Rome Ciampino. The positive aspect of using these data is that it can be considered as cleaned and reliable and has been used to determine weather derivatives payoff in the past. The downside is that it is private with restricted access.

2.2 Some background on temperature models

2.2.1 Ornstein–Uhlenbeck models

Interest on temperature models for weather temperature derivative pricing have arisen in the last years. Even though there exists some alternative modeling (Schiller et al., 2012, Elias et al., 2014, Gülpinar & Çanakoglu, 2017), there seems to be an overall agreement on the capacity of Ornstein–Uhlenbeck processes to model daily average temperatures. These models have largely been studied by Benth and Benth (Benth & Šaltytė Benth, 2007, Benth et al., 2007) and are written as follows:

where |$(W_{t})_{t \geq 0}$| is an independent Brownian motion, |$s$| a seasonality deterministic function, |$\kappa $| a nonnegative mean-reverting parameter and |$\sigma ^{2}$| is a deterministic non-negative function.

There is no clear agreement on the form of |$s$|⁠. Different functions have been suggested from step functions (Taylor & Buizza, 2004) to polynomial functions Franses et al. (2001). However, there is a clear preference for Fourier decomposition (Benth et al., 2007, Šaltytė Benth & Benth, 2012, Hølleland & Karlsen, 2020) which gives the following expression:

where |$\xi = \frac{2 \pi }{365}$|⁠, |$\xi _{k} = k \xi $| and |$K_{s} \in{{\mathbb N}}^{*}$|⁠. Here and through the paper, the time unit is the day.

In this model, |$s$| integrates a trend component which enables to introduce climate change phenomenon into our model and a periodic component through trigonometric functions. There is a large discussion on the convenient |$K_{s} \in{{\mathbb N}}^{\star }$| to choose as well as on the pertinence to introduce interaction terms (Meng, 2018). For simplicity purposes, we take |$K_{s}=1$| which is the most commonly choice—see, e.g. (Benth et al., 2007, Šaltytė Benth & Benth, 2012)—and shows significance at 5% confidence level on our experiments.

The different parameters of the Ornstein–Uhlenbeck process can be determined by conditional least squares estimation (CLSE) developed by Klimko & Nelson (1978), which boils down to minimize

Since we have daily data, we mostly use |$\varDelta =1$| through the paper, unless specified. The quantity

is called the residual of the regression at time |$i\varDelta $|⁠. Note that the same formula for the residual holds true for Model (M). We can equally note that the autoregressive factor |$e^{-\kappa \varDelta }$| derives from the integration of the Ornstein–Uhlenbeck dynamics which is explicit in Equation (C.5) in Appendix  C. The CLS estimators have been studied by Overbeck & Ryden (1997) for the CIR process, Li & Ma (2015) for the stable CIR process and Bolyog & Pap (2019) for Heston-like models. This approach has been used to estimate the parameters of Model (M), See section 3.1.

Similarly, |$\sigma $| represents the deterministic volatility function. For flexibility and periodicity reasons, it can also be modeled thanks to a Fourier decomposition

The coefficients |$\gamma $|’s and |$\delta $|’s are assumed to be such that |$\sigma ^{2}$| is indeed a non-negative function, and we also exclude the trivial case |$\sigma ^{2}(t)=0$| for all |$t\ge 0$|⁠. For |$K_{\sigma ^{2}}=1$|⁠, a straightforward necessary and sufficient condition for having |$\sigma ^{2}$| non-negative is |$\gamma _{0}\ge \sqrt{\gamma _{1}^{2}+\delta _{1}^{2}}$|⁠. For |$K_{\sigma ^{2}}\ge 2$|⁠, a sufficient condition is to assume |$\gamma _{0}\ge \sum _{k=1}^{K_{\sigma ^{2}}} \sqrt{\gamma _{k}^{2}+\delta _{k}^{2}}$|⁠. This decomposition is used for example by Benth & Šaltytė Benth (2007) or Meng & Taylor (2022) in a more evolved model. We also consider it in Model (M) as the function to which the stochastic volatility mean reverts. The annual periodicity is a quite natural feature. Besides, this choice gives a bounded continuous function whose non-negativity is easy to check, which is required for the definition of Model (M). In this paper, |$K_{\sigma ^{2}}$| will be taken equal to 2, but |$K_{\sigma ^{2}}=1$| will also be considered.

2.2.2 Limits of Ornstein–Uhlenbeck models

Ornstein–Uhlenbeck processes, and their discrete form corresponding to Autoregressive Models, present some important limitations.

First, they remove any potential long memory effects as today’s temperature will only depend on the previous day’s. This weakness has been challenged through literature. In particular, Brody et al. (2002) suggest the introduction of Fractional Brownian motions. They suppose an Ornstein–Uhlenbeck model driven by a Fractional Brownian motion of Hurst parameter |$H$|⁠. Such models generate temperature paths with |$(H- \varepsilon )$|-Hölder regularity, |$\varepsilon>0$| being an arbitrary small real number. Following the work of Gatheral et al. (2018), we perform an analysis on our data to check whether we could observe such regularity and consider the below metric

For a |$H$|-Hölder dynamics, |$m(q,\varDelta )$| should behave as |$\varDelta ^{qH}$|⁠. Figure 1 shows that |$\log (m(q,\varDelta ))$| is well approximated by an affine function of |$\log (\varDelta )$| for different values of |$q$|⁠. Therefore, we can compute the coefficient of the regression of |$\log (m(q,\varDelta ))$| by |$\log (\varDelta )$| and estimate |$H$|⁠, see Table 1. We observe that this parameter is very close to |$0.5$|⁠, which corresponds to the regularity of a diffusion driven by a standard Brownian Motion. This justifies why we still consider in Model (M) a diffusion model.

Second, we also contemplated autoregressive models with higher than one autoregressive order. Figure 2 shows the partial autocorrelation plots of |$(\tilde{T}_{t})_{t \geq 0}$| and the residuals |$(Res_{t})_{t \geq 0}$| (see Eq. (4)) for the city of Paris in the first |$1,000$| observed days. While we can consider a second significant correlation term, its coefficient should be pretty small. Higher than two autoregressive terms are clearly non-significant. Therefore, we studied the possibility of having two time delay Ornstein–Uhlenbeck dynamics. However, the second-order autoregressive coefficients were small and unstable. We hence decided to keep a model with a unique autoregressive terms which is coherent with Franses et al. (2001), Campbell & Diebold (2005) and Taylor & Buizza (2004) findings.

Third, as can be observed in Fig. 3, daily temperature present erratic noises with possible volatility autoregression or clustering.

To analyze this, we focus on the residuals (4) and compare the historical ones to those obtained on 9 simulations of the same 40 years period. More precisely, for each simulation, we have plotted the curve |$(Res^{sim}_{(i)},Res^{obs}_{(i)})_{1\le i\le 40\times 365}$|⁠, where |$Res^{sim}_{(i)}$| (resp. |$Res^{obs}_{(i)}$|⁠) is the ordered statistic of the residuals obtained with the simulated (resp. observed) temperature. In the center of the distribution, the points are on the line |$y=x$| in red, which indicates a very good fit by the model. Instead, we see that for small (resp. large) values, the curves are slightly below (resp. above) the red line which indicates that the extreme events produced by the Ornstein–Uhlenbeck model are smaller than the ones observed. The residuals given by Model (1) present noticeable deviations from the historical ones as can be seen in the qqplots: Fig. 4 shows slight skewness but the main observation remains that the observed tails are heavier than the simulated ones, especially for the left tail. Different authors emphasize skewness deviation of residuals as well as volatility clustering (Franses et al., 2001, Meng, 2018). This skewness and heavy tail issues are coped with different methods either thanks to the fitting of a skew-t distribution on residuals (Erhardt et al., 2015) or through GARCH models that enable to capture dependencies on volatility and can lead to skewed residual qqplots. Here, anticipating on our estimation results, we have plotted in Fig. 5 the qqplot on the residuals between those observed and 9 simulations of Model (M) with estimated parameters. We observe again a very good fit of the center distribution, but the curves are slightly above (resp. below) the red line |$y=x$| for small (resp. large) values, which means that Model (M) produces larger extreme events than the ones observed. Thus, compared with the Ornstein–Uhlenbeck model, we note that Model (M) produces heavier tails that are slightly heavier to those observed on our 40 years data set. Thus, Model (M) is more conservative on extreme events, which is an interesting feature when dealing about risk quantification.

2.2.3 GARCH models

The Generalized Autoregressive Conditional Heteroskedasticity (GARCH) models (Bollerslev, 1986) represent a first response to these limits. They enable to integrate stylized features of temperature times series such as skewness, tail heaviness and volatility clustering. Let us define the following temperature dynamics for |$t\in{{\mathbb N}}$|⁠:

with |$\lambda _{k}\ge 0$| for |$1\le k\le p$|⁠, |$\theta _{l}\ge 0$| for |$1\le l\le q$| and |$s$| and |$\sigma ^{2}$| are defined, respectively, by (2) and (5) with parameters |$\gamma $|’s and |$\delta $|’s such that |$\sigma ^{2}(t)\ge 0$| for all |$t\ge 0$|⁠. Thus, |$(\zeta _{t})$| is indeed a non-negative process.

Different authors have worked on the application of these GARCH models particularly on GARCH(1,1) to temperature process fitting (Franses et al., 2001, Taylor & Buizza, 2004, Campbell & Diebold, 2005, Šaltytė Benth & Benth, 2012, Meng, 2018). While these models enable to integrate the complexity of the series, we will see that they also present considerable limits.

First, the analysis of autocorrelation and partial autocorrelation plots on squared corrected residuals shows evidence of the necessity to integrate first-order autocorrelation on the volatility model. This is coherent with the previously cited papers. However, it should be noted that the significance of the first-order autocorrelation term remains small even for subsamples. We can therefore expect that it has a limited impact into the quality of multi-day forecasts and pricing.

Second, the GARCH volatility model and its continuous equivalents integrate an important hypotheses: they suppose that temperature and volatility are driven by the same noise (with a time shift). In Model (M), this would formally correspond to have |$\rho \in \{-1,1\}$|⁠. Given the complexity of meteorological dynamics, there is no reason explaining why temperature and volatility on temperature should be driven by the same noises. Model (M) represents a generalization of average temperature dynamics that enables to integrate more flexibility on the temperature’s volatility.

Finally, Model (M) presents a clear advantage as it provides closed form formula for the pricing of temperature indices deriving from average temperature models. While these closed formulas can only be used for some derivatives depending on the payoff structure, they can be combined with Monte Carlo approaches to speed up derivative’s pricing. In particular, this document will develop the use of control variates methods to reduce the variance of Monte Carlo approaches, see Subsection 4.3.4.

To sum up, Model (M) is aligned with the average temperature models that have been considered in the literature. First, from an Ornstein–Uhlenbeck process, it integrates a seasonal component corresponding to the natural climatology and climate change trend. The Ornstein–Uhlenbeck model also enables to include an autoregressive component upon which agree both literature and statistical tools. However, the residuals of this first model show deviation from normal hypotheses with skewness, tail heaviness and volatility clustering patterns. GARCH models partly answer to these limits and integrate an autoregressive component on the volatility process which is statistically observed. Our model is another natural extension to the Ornstein–Uhlenbeck dynamics that gives a larger flexibility on the volatility process. Model (M) presents two additional advantages. First, unlike ARMA and GARCH models, it is a time-continuous model. While this might not be necessary when we only have daily data, it can be a significant advantage if this model is coupled with a model for energy or commodities that are traded continuously. We could foresee, for example, combining this model with other commodity models to identify hedging opportunities or price hybrid derivatives. Second, Model (M) is an affine model for which efficient pricing methods based on Fourier techniques can be implemented. In particular, we show in Section 4 how to use the Fast Fourier Transform for pricing, which translates into a significant competitive advantage.

3. Fitting Model (M) to historical data

The previous section has motivated the interest of the stochastic volatility model (M) for the temperature dynamics. The pertinence of this new model is however related to our capacity to well estimate its parameters. This section focuses on this challenge as well as on the robustness of the estimation. Contrary to financial derivatives, we do not consider the calibration of our model to market prices. First, as pointed by Weagley (2019), the market volume of weather derivatives is quite low and the transactions are mostly Over The Counter (OTC). Market prices are thus arguable. Besides, the underlying of these contracts is the real temperature, which is not a traded asset. There is therefore no justification of a risk-neutral pricing of weather derivatives. For these reasons, we prefer to estimate our model to historical data and then use it to determine the distributions of weather derivatives payoffs under the real probability.

3.1 Parameter estimation

Such as Overbeck & Ryden (1997) and Bolyog & Pap (2019), we implement the CLSE which consists in minimizing the observed value against the predicted conditional expectation. Under general assumptions, Klimko & Nelson (1978) have shown that CLS estimators are strongly consistent, with a speed of convergence close to |$O(N^{-1/2})$|⁠, where |$N$| is the number of observations.

It can be noted that Bolyog and Pap suggest in Bolyog & Pap (2019), on a time-homogeneous model that is similar to (M), to simultaneously minimize the conditional temperature and volatility expectations, i.e. to minimize

with respect to |$\kappa $|⁠, |$K$| and the (parameterised) functions |$s(\cdot )$| and |$\sigma ^{2}(\cdot )$|⁠. As already remarked in Bolyog & Pap (2019), the estimators of the parameters |$\kappa ,K, s(\cdot ), \sigma ^{2}(\cdot )$| do not involve the values of |$\eta $| and |$\rho $|⁠. We can thus do the estimation without knowing these values and then estimate separately the parameters |$\eta $| and |$\rho $|⁠. Besides, Model (M) corresponds to a special case of the dynamics considered by Bolyog and Pap where the volatility is absent from the mean-reverting term of the temperature process (this corresponds to |$\beta =0$| in (Bolyog & Pap, 2019)). In this case, |$\mathbb{E} [T_{(i +1) \varDelta } | T_{i\varDelta },\zeta _{i\varDelta } ] =\mathbb{E} [T_{(i +1) \varDelta } | T_{i\varDelta }]$|⁠, and the minimization problem is equivalent to minimize |$\sum _{i=0}^{N-1} \left ( T_{(i+1)\varDelta } - \mathbb{E} [T_{(i +1) \varDelta } | T_{i\varDelta } ] \right )^{2} $| and |$\sum _{i=0}^{N-1} \left ( \zeta _{(i+1)\varDelta } - \mathbb{E} [\zeta _{(i +1) \varDelta } | \zeta _{i\varDelta } ] \right )^{2}$| separately, which we do here.

Thus, Paragraph 3.1.1 deals with the estimation of |$\kappa $| and |$s(\cdot )$|⁠, Paragraph 3.1.2 brings on the estimation of |$K$| and |$\sigma ^{2}(\cdot )$|⁠, while Paragraph 3.1.3 focuses on the estimation of |$\eta $| and |$\rho $|⁠.

3.1.1 Parameter estimation for κ and s(·)

To estimate |$\kappa $| and |$s(\cdot )$|⁠, we are thus interested by the following minimization problem:

that brings on temperature process |$(T_{t})_{t \geq 0}$|⁠. Proposition C.1 solves this problem and gives explicit formulas for |$\hat{\kappa }$|⁠, |$\hat{\alpha }_{0}$|⁠, |$\hat{\beta }_{0}$|⁠, |$\hat{\alpha }_{1}$| and |$\hat{\beta }_{1}$|⁠. Table 2 shows the estimated parameters for the 8 European cities. All the parameters are significant at 5% confidence level (their 95% confidence interval does not cross zero). We can also add that there is a certain coherence between the different European cities for the mean-reverting parameter. In addition, given the range of values close to |$0.25,$| we can also derive that the temperature memory effect lasts for around 4 days. These results are aligned with the literature, see, e.g. Benth et al. (2007). We also notice that the estimated values of |$\beta _{0}$| are around |$0.00013$|⁠, which corresponds to a warming of about |$0.5^{\circ }$|C every 10 years.

3.1.2 Parameter estimation for K and σ (·)

An important challenge we face when estimating the parameters of the process |$(\zeta _{t})_{t \geq 0}$| is that the instantaneous volatility process is, per se, unobservable. There is a large literature on this issue, we mention here Aït-Sahalia & Kimmel (2007) and Azencott et al. (2020) that deal with the Heston model. Following Azencott et al. (2020) approach, we approximate the unobservable volatility |$\zeta $| by the series of realized volatilities |$\hat{\zeta }$|⁠. These realized volatilities |$\hat{\zeta }$| correspond to the observed volatility on a time window of |$Q$|-days such that we get

Here, |$\hat{\zeta }_{iQ \varDelta }$| corresponds to the realized volatility on |$[iQ\varDelta , (i+1)Q\varDelta ]$| and we have thus |$I = \lfloor N/Q\rfloor $| different values. The autoregressive factor comes from the integration of the temperature dynamics of Model (M) which is explicit in Equation (C.5) in Appendix  C. The correction factor |$\frac{2 \hat{\kappa }}{1-e^{-2 \hat{\kappa } \varDelta }}$|⁠, which is not present in Azencott et al. (2020), is related to the mean-reverting behaviour of the temperature and is justified by Remark 3.1. Since |$\frac{2 \hat{\kappa }}{1-e^{-2 \hat{\kappa } \varDelta }}\approx _{\varDelta \to 0} \frac 1 \varDelta $|⁠, |$\hat{\zeta }_{iQ \varDelta }$| is close to the usual quadratic variation, but the difference is not negligible as |$\varDelta $| cannot be smaller than one day on observed data.

 

Estimation of |$\sigma ^{2}$| and |$K$|⁠. Once again, we use the CLSE to simultaneously compute the parameters of the deterministic component of the volatility and the mean-reversion coefficient. For this, we would like to minimize

and use Proposition D.1. The convergence of such kind of estimators for a time inhomogeneous CIR is given by Theorem  F.1. However, since the volatility is not directly observed, we minimize the difference of the realized volatility and its conditional expectation given by the previously realised volatilities. Namely, we apply Proposition D.1 to minimize

replacing the volatility |$\zeta _{iQ \varDelta }$| by the realized volatility |$\hat{\zeta }_{iQ \varDelta }$|⁠. This leads to the following estimators of the volatility dynamics of Model (M):

where |$k \in \{1,2\}$| and

and

Table 3 summarizes the numerical implementation of the parameter estimation for our eight cities. We can again observe a coherence between the different cities. It can also be noted that |$\hat{\gamma }_{0}$| has more importance in the |$\sigma ^{2}$| than the trigonometric components. Finally, the mean reverting parameter |$K$| is more unstable than |$\kappa $| through the cities. Its influence can differ from 1 to 7 days depending on the city.

We also tested the choice of setting |$K_{\sigma ^{2}}=1$|⁠. Table 4 summarizes the estimation of the different parameters for |$K_{\sigma ^{2}}=1$|⁠. We can see that the impact on the estimations of |$\gamma _{0}$|⁠, |$\gamma _{1}$|⁠, |$\delta _{1}$| and |$K$| is small or nill. Hence, while all the parameters |$\gamma _{2}$| and |$\delta _{2}$| are significant, we find that setting |$K_{\sigma ^{2}}=1$| or |$2$| has a rather small impact in the reasoning that follows.

In addition, Fig. 6 shows the plots corresponding to the observed volatility and estimated seasonality on volatility. We can see that while the seasonality does not seem negligible |$\sigma ^{2}$| is far from completely explaining the observed volatility. Indeed, we observe important fluctuations around |$\sigma ^{2}(t)$| on the dynamics of |$\zeta $|⁠.

3.1.3 Estimation of the parameters η² and ρ

Estimation of the volatility of volatility |$\eta ^{2}$| While in the previous section we use the conditional expectation to estimate |$\sigma ^{2}$| and |$K$|⁠, here the idea is to implement a similar approach based on conditional variance. Namely, we want to solve

Note that Li & Ma (2015) or Bolyog & Pap (2019) do not study the properties of such conditional second moment estimators. Here, we show that in Theorem  F.2, the convergence of such estimator for a time-dependent CIR process, when the values |$\zeta _{i\varDelta }$| are directly observed.

Again, we have to work with the estimated volatility |$\hat{\zeta }_{iQ\varDelta }$| since we do not directly observe the process |$\zeta $|⁠. We then apply Proposition E.1 with the previously estimated parameters |$\hat{K},\hat{\gamma },\hat{\delta }$|⁠. This leads to the following estimator:

where |$\hat{\vartheta }^{T}$| and |$\hat{X}^{\prime}_{iQ\varDelta }$| are defined by (12) and

with

and

and |$(\hat{\gamma }_{0}, \hat{K}, \hat{\gamma }_{k}, \hat{\delta }_{k} )$| defined by Equation (11).

The numerical application of the above formula applied to our dataset gives the values collected in Table 5. First, comparing Tables 3 and 5 enables to observe the two components of the volatility dynamics, the mean reversion component and the pure noise. They both have comparable magnitude and therefore both contribute significantly to the volatility dynamics. The mean reverting coefficients are rather small, which is consistent with the fluctuations observed in Fig. 6. Second, Table 5 shows a certain coherence in terms of magnitude for all the European cities.

Estimation of the correlation |$\rho $| The last step consists in estimating the correlation |$\rho $|⁠. The idea this time is to use conditional covariance to estimate this parameter, and Proposition E.2 gives the minimizer of the following problem:

Again, since we do not observe the volatility, we use Proposition E.2 with the estimated volatility |$\hat{\zeta }_{iQ\varDelta }$| and the previously estimated parameters |${\hat{\kappa }},\hat{\alpha },\hat{\beta },{\hat{K}},\hat{\gamma },\hat{\delta },\widehat{\eta ^{2}}$| given by Proposition C.1, (11) and (13). This leads to

where |$\hat{\lambda }$| and |$X_{i\varDelta }$| are defined by (C.2) with |$N:= I-1$|⁠, |$\hat{\vartheta }$| and |$\hat{X}^{\prime}_{iQ\varDelta }$| are defined by (12), and |$\hat{Y}^{\prime}_{i Q\varDelta } = \theta _{0}^{\prime\prime} + \phi _{0}^{\prime\prime} \hat{\zeta }_{iQ\varDelta } + \sum _{k} \theta _{k}^{\prime\prime} \sin (\xi _{k} i Q\varDelta ) + \sum _{k} \phi _{k}^{\prime\prime} \cos (\xi _{k} i Q \varDelta )$|⁠, with

and

The numerical application of the above formula applied to our data sets gives the values collected in Table 6. We observe that the correlation is close to zero for all the cities. Therefore, for simplification purposes, we will consider on the following of this document that |$\rho $| equals zero. This finding also questions the pertinence of GARCH model (6) that corresponds to |$\rho \in \{-1,1\}$| since the temperature and its volatility are driven by the same noise.

3.2 Robustness of estimators

In the previous subsection, we have obtained conditional least squares estimators for the different parameters of Model (M). All the above expressions have been computed by discretizing the processes |$(\tilde{T}_{t})_{t \geq 0}$| and |$(\zeta _{t})_{t \geq 0}$|⁠. Theoretically speaking, Overbeck & Ryden (1997) and Bolyog & Pap (2019) have proven the convergence of the CLS estimators for the CIR and a generalized Heston model. Their proof is mainly based in ergodic arguments. In Appendix  F, we prove the consistency of CLSE for a time inhomogeneous CIR process.

However, to estimate the parameters of the volatility dynamics, we have approximated the unobservable volatility |$\zeta $| by the realized volatility |$\hat{\zeta }$|⁠. Azencott et al. (2020) have deeply studied the convergence modes of the volatility process |$\hat{\zeta }$| to the instantaneous process |$\zeta $| as well as the estimated realized estimators |$(K, \sigma , \eta )$| under a classic Heston framework. Under boundary and continuity hypothesis on |$T$| and |$\zeta $|⁠, uniform convergence of |$\hat{\zeta }$| to |$\zeta $| over |$[0,T]$| in |$L^{2}$| is proven. Probability convergence of estimators is also proven for moments-based estimators. The extension of these convergence properties to our particular model is left as a further work. In this section, we test numerically the robustness of our estimators and check their accuracy on simulated data.

3.2.1 Methodology

Robustness of the estimators is checked through simulated data. Essentially, we simulate data series with the model that have the same length as our data set (40 years) and we check that we find back the parameters by using the CLS estimators. The detailed methodology is presented in Appendix  B.

As a first qualitative check, we have represented in Fig. 7 an example of simulated temperature for Paris. The simulated paths look similar to the observed ones. In the following paragraphs, we will test our capacity to estimate the values of the parameters and discuss the choice of averaging-time windows |$Q$|⁠.

3.2.2 Estimation of parameters κ, s(·) related to the simulated temperature T

The estimation of |$\kappa $| and |$s(\cdot )$| is a priori easy since it relies on the temperature that is directly observable. Table 8 summarizes the estimators of the parameters related to the temperature. Figure 8 corresponds to the related temperature plot for Stockholm and Paris. We can see that all the estimated parameters remain very close to the original values. We can thus conclude that the estimation of the parameters of the temperature is robust enough to be reliable.

3.2.3 Estimation of K, σ²(·) and η² on the simulated realized volatility |$\hat{\zeta }$|

Let us recall that instantaneous volatility |$\zeta $| is not observable and we estimate it by |$\hat{\zeta }$| defined in Equation (8). In this paragraph, we apply the same process as in Section 3.1 to simulated data.

We first focus on the effect of |$Q$|⁠, i.e. the size of the averaging window (8), on the estimation of parameters. For simplicity purposes, we start by setting the trigonometric coefficients of |$\sigma ^{2}$| to |$0$| and study the impact on |$K$| and |$\eta ^{2}$|⁠. For each time window |$Q$| and city, we perform |$50\,000$| simulations of daily volatility |$\zeta $|⁠, by using (B.1). We average this series trough the time window |$Q$| and then estimate the corresponding |$(K, \eta ^{2})$|⁠. This exercise enables to analyze the influence of the time window |$Q$| on the capacity to well estimate |$(K, \eta ^{2})$|⁠.

Tables 9 and 10 represent different estimates of |$(K, \eta ^{2})$| depending on window width |$Q$|⁠. First, we can observe that for small |$Q$|⁠, |$(K, \eta ^{2})$| is overestimated due to the preponderance of the noise related to volatility estimation. On the contrary, for large values of |$Q$|⁠, |$(K, \eta ^{2})$| is underestimated due to the averaging effect. Therefore, there is a need to find a trade-off between the effect of the noise and the effect of the window averaging. Depending on the city the most efficient |$Q$| can range from 5 to 10.

We now focus on the impact of |$Q$| on the estimation of the function |$\sigma ^{2}$|⁠. We have plotted in Fig. 9 the estimated trend with the original one for |$Q=5$| and |$Q=12$|⁠. We see that |$\sigma ^{2}$| is correctly estimated, independently of |$Q$|⁠. From this numerical study, the choice of |$Q=10$| appears to be quite reliable. Of course, as one may expect, the parameters are not as well estimated as for |$\kappa $| and |$s(\cdot )$|⁠. They still however give the correct magnitude of the parameters, which is acceptable for risk management. Since we are then interested in evaluating derivative products, we will then analyze the sensitivity to these parameters, see Section 4.4, which can be done efficiently in Model (M).

4. Application to pricing weather derivatives

The previous sections concentrate on modeling the daily average temperature and on the estimation of the parameters of the model. However, the final objective of our model is to better assess the risk related to weather temperature derivatives. This section will focus on how we evaluate the average payoff of these derivatives and how Model (M) improves our capacity to understand their risk.

4.1 Temperature derivatives

4.1.1 Average temperature indices

Temperature derivatives are financial products used to hedge weather risk. The covers are often based on an index corresponding to a proxy of the buyer’s financial risk. This index corresponds to an aggregate of a more granular meteorological parameter which in this document is the average daily temperature. There exist different possible indices, the main ones being HDD, CDD and CAT

where |$T_{t}$| corresponds to the average daily temperature on day |$t$|⁠, |$T_{b}$| to a base temperature, |$t_{1}$| to the inception date and |$t_{2}$| the exit date of the contract. We will call risk period the time period between |$t_{1}$| and |$t_{2}$|⁠.

Physically speaking, the HDD corresponds to the cumulative degrees needed to heat a given building. Hence the base temperature |$T_{b}$| corresponds to the temperature at which heating is probably switched on and the cumulative HDD measures the demand on energy of the building, the base temperature |$T_{b}$| varies depending on the country. In the EU, |$T_{b}$| is often taken as equal to |$15.5^{\circ }$|C, while in the USA, it corresponds to |$65^{\circ }$|F. Symmetrically, CDD corresponds to the energy demand for air conditioning. Finally, CAT corresponds to cumulative temperature degrees which is related to the energy demand between |$t_{1}$| and |$t_{2}$|⁠.

In the following of this document, we will focus on the HDD index; however, the methodology presented can be applied to all average temperature related indices. Particularly for options on the CAT, the FFT approach would enable to get a very efficient pricing method.

4.1.2 Payoff function

Weather derivatives are used to hedge weather risk. They trigger a payment depending on an aggregate temperature index. The payment is defined given a payoff structure. Standard payoff structures correspond to capped put or call options applied to the aggregate index. For simplification purposes, we will suggest the payoff structure

As we want to price this kind of instruments, our objective is to understand the characteristics of the payment distribution (expectation, VaR and CVaR) under the real-world probability. In particular, we consider the average payoff of the derivative

where |$D(t_{0}, t_{2})$| is a discount rate that will be taken equal to |$1$| in this paper. Note that this is not a fair price: there is no market dealing HDD continuously, and therefore, the classical pricing theory of Black and Scholes does not apply. The calculation of the average payoff (18), as well as other indicators on the distribution of |$\min ((HDD - HDD_{strike})^{+}, L)$| such as the variance and quantiles, is used in practice to propose a price over the counter. Thus, a very accurate evaluation of (18) for a given model is not really at stake: one is more interested in evaluating risk and how the average payoff may change under stressed parameters.

Finally, there is no consensus on how to choose |$HDD_{strike}$|⁠. However, it is a market practice to use quantiles, and particularly historical quantiles of the index, to define this strike. In the present paper, we consider the 90% quantile which is within market practices.

4.2 Monte-Carlo approach

A first pricing approach to identify the distribution of payments is to simulate temperature paths based on the discretization schemes in (B.1) for |$\varDelta =1$|⁠. We proceed as follows.

• Simulate temperature paths starting from the pricing date |$t_{0}$|⁠, the day until which we can observe temperature data, to the expiration date |$t_{2}$|⁠.

• Compute simulated HDD between |$t_{1}$| and |$t_{2}$| for each of the paths and obtain an HDD distribution.

• Either fix an arbitrary |$HDD_{strike}$| or choose a quantile to select the moneyness of the structure.

• Deduce the payment distribution.

• Compute payment distribution characteristics: mean, VaR and CVaR.

Figure 10 shows the results of this method for Paris temperature in 2019. Contracts last 1 month and are computed 30 days in ahead, i.e. |$t_{1}-t_{0}=30$|⁠. We consider the payoff function (17) with |$L=+\infty $| and |$HDD_{strike}$| set to the 90% empirical quantile of the |$HDD$| distribution obtained with Model (M). We perform |$50\,000$| Monte Carlo simulations for the Ornstein–Uhlenbeck model and for Model (M). All these choices are challenged in the following sections.

From Fig. 10, we can see that both Model (M) and the Ornstein–Uhlenbeck model lead to similar expected payoffs for winter months while Model (M) states higher expected payoffs for summer months. Model (M) tends to have heavier left tails which are particularly visible in summer months. This explains slightly higher mean payoffs for Model (M) during these months. Nevertheless, it should be noted that these derivatives are mainly sold for winter months to cover against cold waves. During these months, both models give similar mean payoffs.

In terms of risk metrics, we compute conditional value at risks at 95% for both models. We can see that Model (M) presents again heavier tails. Thus, as already noticed from Figs 4 and 5, Model (M) is more conservative and reduces the problem of underestimating rare events related to the Ornstein–Uhlenbeck Gaussian framework.

4.3 F‌FT approach

This section explores an alternative methodology for HDD pricing thanks to the FFT approach developed by Carr & Madan (1999). This method has been widely used in the literature for pricing, we just mention here the recent work of Benth et al. (2021) for an application close to ours.

4.3.1 The characteristic function

In order to apply FFT pricing, we first calculate the characteristic function of |$(\tilde{T}_{t},\zeta _{t},\int _{0}^{t}\tilde{T}_{s} ds)$|⁠. A semi explicit formula is available because of the affine structure of Model (M).

 

 

 

4.3.2 Approximation of the characteristic function

We now discuss the approximation of the characteristic function (19). To do so, we consider a time step |$\delta>0$|⁠, and we will assume that |$t = t_{k} = k \delta $| and |$t^{\prime} = t_{l} = l \delta $|⁠. Note that the function |$a_{1}$| is fully explicit and does not need to be approximated. We use the trapezoidal rule to integrate the function |$a_{0}$|

The main issue may come from the discretization of the Riccati equation which may lead to instabilities if it is not well handled. Here, we take advantage of the fact that an explicit solution of (20) is known for |$\kappa =0$| and |$u_{3}=0$|⁠, see, e.g. (Alfonsi, 2015, p. 101),

with

Thus, to solve (20), we freeze on each interval |$[t_{k},t_{k+1}]$| the value of the time inhomogenous term to its value at |$t=\frac{t_{k}+t_{k+1}}{2}$| and use the explicit formula. This is the midpoint method that leads formally to a convergence of order |$O(\delta ^{2})$|⁠. This leads to

where

We implement the three functions in (22) which enable us to deduce the characteristic function of |$(u_{1} \tilde{T}(t)+u_{2}\zeta _{t})_{t\geq 0}$| where |$(u_{1},u_{2}) \in \mathbb{C}^{2}$|⁠. Figure 11 shows the characteristic function of |$(T_{t})_{t\geq 0}$| calculated with the approximation (22). It is compared with the Monte-Carlo estimator obtained with simulated path using (B.1) (we have used here independent simulations for each values of |$t^{\prime}$|⁠). We can see that both methods give close results, which validates the relevance of the approximation.

4.3.3 F‌FT for pricing HDD and related options

Once we have the characteristic function of |$(T_{t})_{t\geq 0}$|⁠, we can price HDD by using Fourier inverse transform techniques. Here, we adapt the approach of Carr & Madan (1999) that uses the FFT in order to calculate the cumulative distribution function of |$T_{t}$|⁠. This allows us to calculate then easily the average value different types of bespoke options.

We first focus on the calculation of |${{\mathbb E}}[(T_{b}- T_{t})^{+}]$| (resp. |${{\mathbb E}}[\min ((T_{b}- T_{t})^{+},L)]$|⁠) that can be seen as the average price of a ‘daily HDD’ (resp. capped daily HDD). We will use this naming later on. The characteristic function |$\varPhi (u)={{\mathbb E}}[e^{iu \tilde{T}_{t}}]$| is given by Proposition 4.1. We have

We approximate this cumulative density function by using the Gil-Pelaez inversion formula (Alfonsi, 2015, Theorem 4.2.3 p. 104)

Let |$\delta _{x},\delta _{v}>0$| be such that |$\delta _{x} \delta _{v}=\frac{2\pi }{N}$| (one can take for example |$\delta _{x}=\delta _{v}=\sqrt{\frac{2\pi }{N}}$| or |$\delta _{x}=\frac{L}{N-1}$| for the capped option so that |$x_{0}$| defined below is equal to |$T_{b}-s(t)-L$|⁠). We define

so that |$x_{N-1}=T_{b}-s(t)$|⁠, and use the following approximation:

since |$ v_{j+1/2}x_{k}= 2\pi \frac{jk}{N}+(j+1/2)\delta _{v}x_{0} + \frac 12 \delta _{v} k \delta _{x}$|⁠. This amounts to use the midpoint rule and to truncate the integral at |$M=N\delta _{v}$|⁠. Other choices of quadrature are possible but have to be taken in compliance with the FFT. Using (24), we can obtain |$(\mathbb{P}( \tilde{T} \leq x_{k}), 0\le k\le N-1)$| by applying the FFT to |$\left (\frac{ e^{-i (j+1/2)\delta _{v}x_{0}}\varPhi (v_{j+1/2})}{iv_{j+1/2}},0\le j\le N-1\right )$|⁠: these |$N$| values are obtain with a time complexity of |$O(N\log (N))$| (instead of |$O(N^{2})$| with the naive calculation of the sums). We finally approximate the expectation of the daily HDD by

Figure 12 shows the characteristic function and the expected HDD in Paris during January 2019 comparing Monte Carlo simulation and FFT approach. Both graphs display a clear coherence between Monte Carlo simulations and the FFT approach. In this case, Monte Carlo simulations show precision given that we simulate |$50\,000$| scenarios. However, FFT pricing is more precise, smooth and faster. We perform FFT with |$N=2^{17}$|⁠.

We now focus on the pricing on options on HDD. We observe from the distribution in Fig. 12 that we almost always have |$T_{t}\le T_{b}$| in January, otherwise we would notice a Dirac mass at |$0$|⁠. In fact, with the standard strike |$T_{b}=15.5^{\circ }$|C, we mostly have |$T_{t}\le T_{b}$| during winter, and therefore, |$HDD\approx (t_{2}-t_{1}+1)T_{b}-CAT$|⁠, so that the average value of the option (17) can be approximated by

The problem of computing the right hand side is then similar to the pricing of daily HDD (23): the underlying is now |$CAT$| instead of |$T_{t}$|⁠. We can thus calculate the average payoff, provided that we know the characteristic function of the CAT |$\varPhi (u)=e^{iu \sum _{t=t_{1}}^{t_{2}} s(t)} {{\mathbb E}}[e^{iu \sum _{t=t_{1}}^{t_{2}} \tilde{T}_{t}}]$|⁠. To do so, it is possible to use formula (19) inductively with |$u_{3}=0$| (using Remark 4.1) in order to calculate |${{\mathbb E}}[e^{iu \sum _{t=t_{1}}^{t_{2}} \tilde{T}_{t}}|{{\mathcal F}}_{t_{2}-\ell }]$| for |$\ell =1,\dots ,t_{2}-t_{1}$| and then |$\varPhi $|⁠. This is however cumbersome, and we prefer to make the following approximation:

We apply Proposition 4.1 with |$u_{1}=u_{2}=0$| and |$u_{3}=u$|⁠, and get |${{\mathbb E}}[e^{iu \int _{t_{1}}^{t_{2}+1} \tilde{T}_{t} dt }|{{\mathcal F}}_{t_{1}}]=\exp (a_{0}(t_{1},t_{2}+1)+i u \frac{1-e^{-\kappa (t_{2}+1-t_{1}) }}\kappa \tilde{T}_{t_{1}}+a_{2}(t_{2}+1-t_{1})\zeta _{t_{1}})$|⁠. Hence, for |$t_{0} \leq t_{1} \leq t_{2}$|⁠,

To obtain |$\check{a}_{0}, \check{a}_{1}, \check{a}_{2}$| and the above characteristic function, we apply a second time Proposition 4.1 with |$u_{1}= u \frac{1-e^{-\kappa (t_{2}+1-t_{1}) }}\kappa $|⁠, |$u_{2}=-ia_{2}(t_{2}+1-t_{1})$| and |$u_{3}=0$|⁠. Figure 13 compares CAT distribution obtained with Monte Carlo and FFT inverse methods for the month of January 2019. We can observe a good fit between both methods.

This section has focused on the capacity of FFT method to compute explicit formulas for options on |$T$| and |$CAT$|⁠. This methodology enables to get rid of the computation burden of Monte Carlo simulations. However, not all the indices or payoff functions can be explicited with FFT. In particular, derivatives that integrate double non-linearities, like put or call payoff functions applied to HDD, cannot be explicitly computed with the FFT method. The next section will focus on how to use FFT results to increase the performance of Monte Carlo simulation for such cases.

4.3.4 Control variates method for Monte-Carlo

In practice, the approximation (26) is precise when |$\mathbb{P}(T_{t}>T_{b})$| is close to zero. However, when this probability is small but not negligible, the approximation may not be enough precise. However, we can use the calculation above to run a Monte-Carlo method with the control variable |$\min ((HDD-HDD_{strike})^{+},L) - \lambda \min ( ((t_{2}-t_{1}+1)T_{b} -HDD_{strike}-CAT)^{+},L) $| in order to calculate the average payoff |${{\mathbb E}}[(HDD-HDD_{strike})^{+} ]$|⁠. Namely, we write (we take here |$L=+\infty $| for simpler notation)

and we chose |$\lambda $| that minimizes |$Var\left [ (HDD-HDD_{strike})^{+} - \lambda ((t_{2}-t_{1}+1)T_{b} -HDD_{strike}-CAT)^{+} \right ]$|⁠, i.e.

The first term of the right hand side of (27) is calculated by using the FFT, while the second one is calculated by Monte-Carlo.

Figure 14 and Table 11 show the results of the implementation of the control variates method. First, from Table 11, we can see that the control variates method enables to decrease the variance up to |$ 2.41 \times 10^{5}$| times, leading to a price computation |$10^{5}$| time faster. Second, we can observe that the performance of the method depends on the correlation between |$(\sum _{t=t_{1}}^{t_{2}}(T_{b}-T_{t})^{+} - HDD_{strike})^{+}$| and the control variable: the more correlated they are, the more VR we obtain (and the more approximation (26) is valid). Hence for the winter months, the computation is more effective, which coincides with the months for which such options are sold. In Fig. 14, we see that the confidence interval for winter months is considerably narrower.

To sum up, this section has focused on exploring alternative pricing methodologies. The FFT pricing method enables to bypass Monte Carlo simulations and to get analytical expressions for the expected payoffs of some derivatives. This enables direct expectation computations for some derivatives such as |$CAT$|⁠. For derivatives integrating non-linear indices and payoffs, this method can be combined with the control variates method to decrease the computational cost of the Monte Carlo simulations. In our case, this method enables to considerably decrease the number of required simulations.

4.4 Sensitivity study

This section studies the sensitivity of the pricing to the different parameters that were either imposed or estimated in the previous sections.

Sensitivity to |$\kappa $| We first analyze the sensitivity to the mean reverting parameter of the temperature dynamics. Figure 15 shows that increasing |$\kappa $| has an important effect on the average payoff and the simulated HDD distribution and therefore on the pricing. Indeed when |$\kappa $| increases the volatility loses its importance, the HDD distribution becomes more certain and therefore peaks around an expected average value. Similarly, the quantiles are less spread, and therefore, the strikes based on initially estimated |$\kappa $| are less frequent and mean payoffs shrink.

Sensitivity to |$\eta ^{2}$|  Figure 16 shows different average payoffs and confidence intervals for |$50\,000$| simulations and different values of the volatility of the volatility |$\eta ^{2}$|⁠. We can see that increasing |$\eta ^{2}$| creates more peaked distributions for HDD and heavier tails. This also leads to usually higher prices. Nevertheless, it should be noted that the impact of |$\eta ^{2}$| is marginal in winter months when this product is meant to be sold. In summer months, HDDs only capture extreme temperature left tails and this is when we can see a real impact of the volatility of Model (M). Besides, let recall that the estimation of |$\eta ^{2}$| is sensitive to the choice of |$Q$|⁠. Wrongly estimating this parameter would therefore mainly impact the pricing of derivatives on summer months where the demand of such derivatives is much lower.

Sensitivity to |$K$|  Figure 17 shows different average payoffs and confidence intervals for |$50\,000$| simulations. We can see that increasing |$K$| creates less peaked distributions for HDD and lighter tails. This leads to usually lower mean payoffs when |$K$| increases. Likewise, the impact is marginal in winter months when this product is meant to be sold. This phenomenon is intuitive as we increase the mean reverting term of the volatility, the volatility of the volatility losses weight in the dynamics and the extreme HDDs decrease.

Sensitivity to |$t_{1}-t_{0}$| Now we suppose that we compute the price of the derivative different possible times ahead. Figure 18 shows temperature paths for different observation times |$t_{0}$|⁠, but the same observed temperature and volatility at |$t_{0}$|⁠. The derivative applies between the black vertical lines. We can observe that all the paths end up following the seasonality |$s$|⁠.

From simulated densities in Fig. 18, we can first observe shifts depending on |$t_{1}-t_{0}$|⁠. This significantly impacts the quantiles of these densites and therefore the |$HDD_{strike}$|⁠. Second, we can note that the more ahead we forecast the less information we have. In this case, we can observe that while pricing 20 and 30 days ahead lead to similar average payouts during the risk period, forecasting 5 days ahead significantly impacts the average payoffs and therefore pricing. This element is key to answer the risk of antiselection.

Sensitivity to the moneyness of the product The moneyness of the product has a direct impact on the payoffs distribution as can be observed in Fig. 19. The lowest the |$HDD_{strike}$|⁠, the more HDD we capture in the payoff and the higher the mean payoffs becomes.

To sum up, this section has focused on the sensitivity of the pricing to the different parameters. We particularly show that the parameters related to the temperature dynamics such as |$\kappa $| as well as the moneyness of the payoff function are the ones affecting the most the mean payoffs. Parallely, the distribution of the payoffs and the strike based on quantiles show important sensitivity to the time interval |$t_{1}-t_{0}$|⁠. Finally, the parameters related to the volatily have a relatively lower impact on the payoffs distribution and hence on the pricing.

4.5 Comparison of our pricing methodology with business practices

This section aims to make the bridge between the pricing methodology exposed in this document and the current market practices. In particular, we will compare a pricing based on modeling of the underlying meteorological parameter with a pricing based on historical index modeling.

As described in Schiller et al. (2012) and Jewson et al. (2005), the index-based pricing methodology consists in modeling the independent yearly indices, in this case cumulative HDD.

For computation ease, the index pricing is automatised following the below algorithm.

• (1) Compute historical index, cumulative |$HDD$|⁠, from 1980 to 2018.

• (2) Remove the linear trend in this time series.

• (3) Fit a gamma distribution to these observations through a maximum likelihood method.

• (4) Compute the expected payoff of the fitted index distribution.

While this approach can be simplistic, as we can consider other probability distributions for the index, it syntheses the common market practices.

Figure 20 represents the expected payoffs computed with the two methodologies. |$HDD_{strike}$| corresponds to the 90% simulated quantile with the Monte Carlo method and to the 90% historical quantile for the index modeling method. First, we can observe that there is a coherence between the approaches that give expected payoffs in the same ranges. However, the index modeling approach introduces important instabilities. These instabilities affect both the average payoffs as well as the strike |$HDD_{strike}$|⁠, which is estimated from less than 40 observations.

We also studied the possibility of using the same strikes with both methods. First, using the 90% quantile of |$50\,000$| Monte Carlo simulations leads to slightly more volatile average payoffs for the index model method. However, we feel that it is counter-intuitive to use a yearly index modeling for pricing and a more cumbersome daily index modeling just to get the strikes. Second, we can use historical quantiles for both methodologies as on the left of Fig. 21. In this case, we can see that the winter months seem to be completely overpriced by the business practice, while the summer months are underpriced. On the right of Fig. 21, we compare simulated and historical quantiles. We can see that the more difference we have between these quantiles, the higher the risk of over or underpricing.

To sum up, there exist a clear coherence between the pricing methodology followed in this document and the current market practices as both give prices in the same ranges. However, our approach enables to better quantify the sources of risk, which is crucial in this growing market.

Data availability

The temperature data underlying this article were provided by Speedwell under licence.

References

A. Weather station data description

B. Model (M) simulation and estimator testing

This section explains the algorithm to simulate Model (M). Simulation is first performed to test the robustness of the estimation through the following steps.

• (1) Estimate |$\kappa $| and the seasonality |$s$| from temperature data.

• (2) Estimate the parameters |$K$| and |$\sigma ^{2}$| and then |$\eta ^{2}$|⁠.

• (3) Fix these parameters at the estimated values.

• (4) Generate a simulated instantaneous volatility series |$\zeta $| based on a generalized Ninomiya–Victoir scheme for Cox–Ingersoll–Ross (CIR) processes, and the corresponding temperature series |$T$|⁠.

• (5) Estimate |$(\alpha _{0}, \beta _{0}, \alpha _{1}, \beta _{1}, \kappa )$| on the simulated data and compare with the fixed values.

• (6) Compute realized volatility |$\hat{\zeta }$| for different time lags |$Q$|⁠.

• (7) Estimate |$(\gamma _{0}, \gamma _{1}, \delta _{1}, \gamma _{2}, \delta _{2}, K)$| then |$\widehat{\eta ^{2}}$|⁠, and compare with the fixed values.

The generation of the combined discrete series |$(T_{i\varDelta }, \zeta _{i\varDelta })_{i \in{{\mathbb N}}}$| is performed thanks to the recurrence formula

where |$(Y_{i},Z_{i})_{i\ge 0}$| is an i.i.d. sequence of two independent standard normal variables. Remember that we assume here and in the sequel of the paper that |$\rho =0$|⁠. The first row of (B.1) corresponds to a discretization of the integral following the temperature dynamics in Model (M) for a step |$\varDelta $|⁠. The use of the trapezoidal rule comes from the operator splitting method and allow to get a second order scheme as in (Alfonsi, 2015, Eq. (4.31)). The second row of (B.1) is the Ninomiya–Victoir scheme for CIR processes Ninomiya & Victoir (2008) when freezing the time-dependent coefficients at time |$(i+1/2)\varDelta $|⁠, which preserves the convergence of order 2 of this scheme, see (Alfonsi, 2015, Paragraph 3.3.4). In general (i.e. when |$K \sigma ^{2}((i+1/2)\varDelta ) \geq \frac{\eta ^{2}}{4}$|⁠), |$\phi $| corresponds to

where |$\psi _{K}(t) = \frac{1- e^{-K t}}{K}$|⁠. The case where |$K \sigma ^{2}((i+1/2)\varDelta ) < \frac{\eta ^{2}}{4}$| is handled as in Alfonsi (2010), so that (B.1) is a second-order scheme for the weak error, and thus an accurate approximation of the exact law.

 

C. CLS estimators of the temperature process

Let us consider the following dynamics for the temperature |$(T_{t})_{t \geq 0}$|⁠:

where |$\kappa>0$|⁠, |$s(t) = \alpha _{0} + \beta _{0} t + \alpha _{1} \sin (\xi t) + \beta _{1} \cos (\xi t) $|⁠, |$W$| and |$Z$| are two independent Brownian motions and |$\zeta _{t}$| is a non-negative adapted process such that |${{\mathbb E}}[\int _{0}^{t} \zeta _{s} de]<\infty $| for all |$t>0$|⁠. The goal of this appendix is to compute the conditional least squares estimators of |$(\kappa , \alpha _{0},\beta _{0}, \alpha _{1}, \beta _{1} )$| and to prove the next proposition.

 

 

Let us note here that |$\hat{\lambda }^{T} X_{i\varDelta }$| can then be seen as the estimation of |$\mathbb{E} [T_{(i+1) \varDelta } | T_{i \varDelta } ]$|⁠.

D. CLS estimators of the volatility process

Let consider the volatility of the temperature |$(\zeta _{t})_{t \geq 0}$| follows the below dynamics:

where |$K>0$|⁠, |$\sigma ^{2}$| is a non-negative function with the parametric form given by (5) and |$W$| is a Brownian motion. The goal of this appendix is to compute the conditional least squares estimators of the parameters |$(\gamma _{0}, K,\gamma _{1},\dots , \gamma _{K_{\sigma ^{2}}}, \delta _{1}, \dots ,\delta _{K_{\sigma ^{2}}})$|⁠.

 

For the proof of Proposition D.1, we will need Lemma D.1, whose proof is straightforward.

 

 

Let us note here that |$\hat{\vartheta }^{T} X^{\prime}_{i\varDelta }$| can be seen as the estimation of |${{\mathbb E}}[\zeta _{(i+1)\varDelta }|\zeta _{i\varDelta }]$|⁠.

E. Computation of the CLS estimators of η² and ρ

Let consider the volatility of the temperature |$(\zeta _{t})_{t \geq 0}$| follows the below dynamics for |$t \geq 0$| and |$\zeta _{0} \geq 0$|

as in Model (M) with |$ \sigma ^{2}(t) = \gamma _{0} + \sum _{k=1}^{K_{\sigma ^{2}}} \gamma _{k} \sin (\xi _{k} t ) + \sum _{k=1}^{K_{\sigma ^{2}}} \delta _{k} \cos (\xi _{k} t )$|⁠. We first focus on the conditional least squares estimator of the volatility |$\eta ^{2}$| and assume that the coefficients |$K$| and |$\sigma ^{2}(\cdot )$| are known.

 

For the proof of Proposition E.1, we first state Lemma E.1, which is a straightforward generalization of (D.6).

 

 

We now focus on the conditional least squares estimator of the correlation |$\rho $| for Model (M) and assume that the coefficients |$\kappa $|⁠, |$s(\cdot )$|⁠, |$K$|⁠, |$\sigma ^{2}(\cdot )$| and |$\eta ^{2}$| are known.

 

Let us note that we do not know a priori that |$\hat{\rho } \in [-1,1]$|⁠.

 

F. Strong consistency of CLS estimators for the time-dependent CIR processes

We study in this appendix the strong consistency of CLS estimators of a time-dependent CIR process. This process is implemented in this paper to represent the temperature volatility dynamics. Let us consider the following process:

with |$K,\gamma ,\eta>0$| and |$\theta :{{\mathbb R}}_{+}\to{{\mathbb R}}_{+}$|⁠. We assume that process is observed at discrete times |$(\zeta _{k\varDelta })_{k \in{{\mathbb N}}}$|⁠.

The goal of this appendix is twofold. First, we prove in Theorem  F.1 the consistency of the CLS estimator of |$\gamma $| when other parameters are known and give the rate of convergence. This result complements the one of Overbeck & Ryden (1997) in a time inhomogeneous case. This is a simplification with respect to the estimation of |$K$|⁠, |$\gamma $|’s and |$\delta $|’s in model (M) given by Proposition D.1: we only estimate one drift parameter instead of |$2(K_{\sigma ^{2}}+1)$| drift parameters. This avoids cumbersome calculations, but the same behaviour is expected for the CLS estimators of these |$2(K_{\sigma ^{2}}+1)$| parameters. Second, we prove in Theorem  F.2 the consistency of the CLS estimator of |$\eta ^{2}$| when other parameters are known. This result complements the results of Bolyog & Pap (2019) that only focus on the CLS estimation of the drift part.

By straightforward calculations, we have for |$0\le s \le t$|⁠,

The CLS estimator of |$\gamma $| consists in minimizing |$\sum _{i=0}^{N-1} \left (\zeta _{i\varDelta }- {{\mathbb E}}[\zeta _{(i+1)\varDelta }|\zeta _{i\varDelta }]\right )^{2}$|⁠, i.e.

which leads to

In this appendix, we note |$\hat{\gamma }_{N,\varDelta }$| instead of |$\hat{\gamma }$| to remind the dependence on |$N$| and |$\varDelta $|⁠. This makes clearer the statements of Theorems  F.1 and  F.2 that involve these two quantities.

In the particular case |$\theta \equiv 1$|⁠, we have

The second term is negligible and, following Overbeck & Ryden (1997), we get that the estimator |$\hat{\gamma }_{N,\varDelta }$| is strongly consistent (i.e. |$\hat{\gamma }_{N,\varDelta }\to \gamma $| a.s.) and asymptotically normal (i.e. |$\sqrt{N}(\hat{\gamma }_{N,\varDelta }-\gamma )$| converges in law to a normal random variable) by using the ergodic theorem.

When |$\theta $| is not constant, we can no longer use the ergodic theorem. We will make the proof of consistency under the assumption that |$\theta $| is a bounded function. We lose the asymptotic normality but still have a convergence rate of |$\sqrt{N}$|⁠. We will use the following lemma.

 

 

 

 

We now turn to the Conditional Least Squares estimation of |$\eta ^{2}$| for the process (F.1). We now assume that |$K,\gamma>0$| and |$\theta (\cdot )$| are known. Without loss of generality, we assume that |$\gamma =1$| and consider the minimization problem of

with respect to |$\eta ^{2}$|⁠. By using Equation (F.2) and Fubini theorem, we get

The minimization of

then leads to the following estimator:

with |$\varTheta ^{1}_{i}={K \gamma } \int _{i\varDelta }^{(i+1)\varDelta } \theta (v) e^{-K((i+1)\varDelta -v)}dv$| and |$\varTheta ^{2}_{i}=\int _{i\varDelta }^{(i+1)\varDelta } \theta (v) e^{-K((i+1)\varDelta -v)}(1-e^{-K((i+1)\varDelta -v)}) dv$|⁠.