Valuation and Risk Management of Weather Derivatives:
The Application of CME Rainfall Index Binary Contracts
天氣衍生性商品之評價與風險管理：CME 雨量指數二元 式合約之應用
Shih-Kuei Lin, Department of Money and Banking, National Chengchi University 林士貴 / 國立政治大學金融學系
Ming-Che Chuang, Department of Finance, Feng Chia University 莊明哲 / 逢甲大學財務金融學系
Dong-Jie Fang, Department of Money and Banking, National Chengchi University 方東杰 / 國立政治大學金融學系
Received 2019/8, Final revision received 2020/11 Abstract
In this paper, we explore and analyze the valuation and risk management of rainfall index binary contracts, a type of precipitation derivative issued by the Chicago Mercantile Exchange (CME). We describe the underlying rainfall index with the occurrence model, which is built on a first-order, two-state Markov chain, and with the magnitude model based on mixed exponential distribution. To capture the seasonality characteristics, we describe the parameters of these two models with the truncated Fourier series. Since the weather derivatives market is incomplete due to the essence of its product, we value the rainfall index binary options with the Esscher transform and calibrate the market price of risk (MPR) with real market data. After analyzing the temporal behavior of the MPR, we find that the investors could have more accurate estimations of the rainfall index when approaching the end of the contract period or when entering the accumulation period. We also find that rather than speculators, the market participants are mainly hedgers, which may explain the shrinking of the precipitation derivatives market.
【Keywords】precipitation derivatives, Markov chain, truncated Fourier series, Esscher transform, market price of risk
摘 要
本文探索並分析了芝加哥商品交易所(CME) 發行之雨量衍生性商品—雨量指數二元
式合約的評價與風險管理。作為標的之雨量指數由兩種模型所刻畫：以一階兩狀態馬 可夫鏈建立的降雨發生模型，以及以混合指數分配建立的雨量強度模型。上述兩模型 之參數均以截斷傅立葉級數擬合，以捕捉降雨的季節性特徵。鑒於天氣衍生性商品市
場之不完備性（肇因於商品本質），本文利用Esscher 變換對選擇權進行評價，並以
CME 之真實市場價格校準獲得市場風險價格。在進一步分析市場風險價格的時間行 為後，本文發現雨量指數二元式合約之市場參與者會隨著合約到期時間迫近或進入契
約的累積期間，而對標的指數有著日益準確的估計; 另本研究也發現這個市場的參與
者大多為避險者而非投機者，這或許能解釋雨量衍生性商品市場萎縮的原因。
【關鍵字】雨量衍生性商品、馬可夫鏈、截斷傅立葉級數、Esscher 變換、市場風險價格
領域主編：黃泓人教授
1. Introduction
Beginning from the late twentieth century, countries around the world started to recognize the urgency and impact of climate change. Moving into the first half of twenty- first century, extreme climate change has become the “new normal” of our daily life.
The key findings of the New Climate Economy (NCE) 2018 report shows that “disasters triggered by weather-and climate-related hazards were responsible for thousands of deaths and US$320 billion in losses in 2017.” (NCE 2018). António Guterres, United Nations Secretary-General also mentions in the 2017 United Nations Climate Change Annual Report that “climate change is the defining challenge of our time.” (UNFCCC 2018). This challenge gradually becomes the “new normal” in the world. On September 23, 2019, government heads, business leaders, subnational actors, youth, indigenous people, and other civil society stakeholders all attended the 2019 UN Climate Action Summit, where bold announcements and commitments were made to reduce emissions, strengthen climate resilience, and mobilize political will for the Paris Agreement. Besides, many scholars pay attention to the climate finance recently, such as Hong, Karolyi, and Scheinkman (2020), Barnett, Brock, and Hansen (2020), and Choi, Gao, and Jiang (2020).
How can investors deal with the extreme climate challenge in finance? Are there any financial instruments to hedge the weather risk? The weather derivatives, which were first developed by energy companies.1 In the financial realm, the market responded to the impact of climate change by developing weather derivatives as early as in 1997, when companies from the old economy started to report loss due to extreme weather conditions.
According to Jeucken (2010), and Müller and Grandi (2000), weather derivatives are widely used for hedging purposes by industries whose profits are highly affected by weather conditions, or by hedge funds for speculative purposes. The payoff of a weather derivative depends on the underlying weather indices such as temperature, wind speed, or precipitation (including rainfall and snowfall).
Since 1999, the Chicago Mercantile Exchange (CME) has played an important role in the weather derivatives market, which has issued temperature, frost index, and snowfall derivatives. Moreover, the CME also issued rainfall index derivatives to provide an
1 Three over-the-counter (OTC) transactions between Willis Group Holdings, Koch Industries and Enron Corporation in 1997 were generally viewed as the beginning of the weather derivatives market.
Among these transactions, two were weather futures traded between Koch and Enron, while the other was traded between Koch and PXRE Reinsurance Company with Willis as the intermediary.
instrument to hedge the risk of rain; and the CME Hurricane Index (CHI) derivatives due to the devastating impact of Hurricane Katrina.
Because of the great loss caused by heavy rainfall disasters2 in Asia in the summer of 2018, the research about the precipitation derivatives can help better deal with these situations in the future. Besides, Lin, Chung, and Yeh (2016), who summarize several pricing methods of weather derivatives published in Taiwan, consider it worthwhile to discuss the valuation of weather derivatives. Therefore, we focus on the valuation and risk management of precipitation derivatives in this paper, hoping to provide reference for the future development of weather derivatives in Taiwan. Though the precipitation can be roughly divided into rainfall and snowfall, the modeling of them is similar and we choose rainfall as the underlying index and take CME Rainfall Index Binary Contracts as a valuation example.
Usually, the pricing of weather derivatives includes two steps: first, modeling the behavior of underlying weather index and secondly, pricing the derivatives with a variety of methodologies.
In this paper, we choose the occurrence model to describe whether it rains or not by using a first-order, two-state Markov chain with four transition probabilities p00, p01, p10, and p11, which is a widely-used model in valuation and actuarial science as in Lin, Hsu, and Chen (2012). We also choose the magnitude model to describe the daily amount of precipitation with mixed exponential distribution having α, β1 and β2 parameters.
For all parameters we mentioned above, including four transition probabilities and three parameters of mixed exponential distribution, we use the truncated Fourier series to describe the seasonal characteristics. The coefficients are estimated by maximum likelihood estimation (MLE). Past literature such as Cao, Li, and Wei (2004), Woolhiser and Pegram (1979), Wilks (2011) and Cabrera, Odening, and Ritter (2013) all witness the application of this type of model.
Compared with other models, like the binomial generalized linear models (GLMs) in Shah (2017), the jump Markov process based on the pulse Poisson process in Carmona
2 On August 27, the hourly cumulative precipitation in Setagaya, Tokyo reached 110 mm, causing the subway being forced to stop. The heavy precipitation in Chungcheong-do, South Korea in August was observed at 159 mm, causing many roads and houses to flood. The torrential rainfall disaster on August 23made cities of southern Taiwan drown in water for days. The Kaohsiung Water Resource Bureau reports an estimation of NT$5.56 billion cost for flood control and rainwater drainage planning.
and Diko (2005), the Ornstein-Uhlenbeck process driven by the Hougaard Lévy process in Noven, Veraart, and Gandy (2015), and the pure jump model embedded in an enlarged filtration framework in Hess (2016), the model we use in this paper is straightforward and comprehensive, and also captures the significant seasonal features in precipitation.
Since the weather market is incomplete, we cannot simply derive the pricing formula like that in the Black-Scholes model (1973). Several methodologies have appeared in prior literature. As mentioned in Dorfleitner and Wimmer (2010), Odening, Musshoff, and Xu (2007) and Shah (2017), burn analysis is a simple but rough approach, which directly prices the weather derivatives with the historical average of the underlying indices. Cao et al. (2004) and Shah (2017) simulate rainfall paths by the Monte Carlo simulation and price the rainfall derivatives with the average of simulated payoffs. Brockett, Wang, Yang, and Zou (2006), Xu, Odening, and Musshoff (2008), and Carmona and Diko (2005) all propose the applying of the utility indifference approach, which is widely used in pricing derivatives with non-tradable underlying assets. However, we choose to follow the approach in Cabrera et al. (2013) and Noven et al. (2015), that is, pricing the precipitation derivatives with the Esscher transform, which was first introduced by Esscher (1932) and been widely used as a risk premium principle in actuarial science.
Compared with other approaches, the Esscher transform has great advantages as shown by several past studies. Bühlmann (1980) proves that the Esscher transform could be derived as the Pareto-optimal solution in a market situation where all participants are characterized by an exponential utility function, and all risks are stochastically independent. Kremer (1982) discusses a characterization of the Esscher transform and proves it yields the closest distance and discrimination (usually called the Kullback- Leibler divergence) between distributions Q and P when calculating net premium. Frittelli (2000) shows that the application of the Esscher transform in option pricing minimizes the relative entropy between the physical measure ℙ and an equivalent martingale measure ℚ, and also equivalently maximizes the exponential expected utility.
There are also past studies consider other factors in pricing precipitation derivatives.
Hess (2016) prices precipitation swaps and futures with customized approximation procedures and extends a multi-location model. Härdle and Osipenko (2017) develops a dynamic programming approach in pricing baskets of weather derivatives under default risk on the issuer side in the OTC market. Cramer, Kampouridis, Freitas, and Alexandridis (2017) compares seven machine learning methods in pricing rainfall derivatives.
Building on the basis of past studies, this paper makes several contributions to the
precipitation derivatives literature. First, we extend the methodologies in Cabrera et al.
(2013) to price CME rainfall index binary contracts, which have never been valued before.
Second, we discuss several observations about the temporal behavior of the MPR θ during the trading period. We especially probe into the risk attitude of the market participants by interpreting the MPR with the definition of the Arrow-Pratt measure of Absolute Risk Aversion (ARA) and sensitivity analysis. Accordingly, we find that the investors could have a more accurate estimation of the rainfall index during the approaching of trading day or when entering the accumulation period. Besides, since many papers like Wu and Chung (2010), Chang, Lin, and Yu (2011), Lo, Lee, and Yu (2013), and Lo, Chang, Lee, and Yu (2020) have focused on the catastrophe risk management, we also find that the influence of local climate conditions and extreme catastrophes plays an important role in the precipitation derivatives market. Last but not least, we find that the investors mainly enter the market for hedging the weather risk, rather than speculating, which may explain the shrinking of this particular derivatives market.
In the remaining of this paper, we present the models of precipitation in section 2, and the model of pricing derivatives in section 3. We then discuss the results of empirical analysis in section 4, and draw our conclusions in the last section.
2. The Models of Precipitation
Since models for different types of precipitation are similar, we mainly take the daily rainfall model as an example, which is widely applied in Woolhiser and Pegram (1979), Cao et al. (2004), Leobacher and Ngare (2011), and Cabrera et al. (2013). Based on hydrological study, there are four characteristics of daily rainfall:
(1) The probability of rainfall occurrence obeys a seasonal pattern;
(2) The sequence of rainy and sunny days follows an autoregressive process;
(3) The amount of precipitation on a rainy day varies with the seasons;
(4) The volatility of the rainfall amount also changes seasonally.
Straightforwardly, the first two characteristics describe the occurrence of rainfall, while the other two characteristics describe the magnitude. To model the dynamic process, we build an occurrence model to describe the presence or absence of rainfall, and a magnitude model to feature the amount of rainfall. The daily rainfall index Rt3 at time t is
3 The daily rainfall index Rt is defined in the CME Rainfall Index Binary Contracts in Table A.1, Appendix A.
measured as the product of the occurrence model Xt and the magnitude model Mt.
Rt = Xt × Mt . (2.1)
2.1 Occurrence Model
The occurrence model is defined as a zero-one process, in which “0” means no rain and “1” means rain:
0, day t it does not rain, Xt =⎧
⎨⎩ (2.2) 1, day t it rains.
As shown in Figure 1, Xt is a first-order, two-state Markov chain. The Markovian property implies that the transition probability of whether it rains or not on day t only depends on the rainfall situation on day t-1, as listed from (2.3) to (2.5).
P{Xt | Xt-1, Xt-2, …, X2 } = P{Xt | Xt-1 }, (2.3) p00 = P{Xt = 0 | Xt-1 = 0 }
p01 = P{Xt = 1 | Xt-1 = 0 } = 1 - p00, (2.4) p10 = P{Xt = 0 | Xt-1 = 1 }
p11 = P{Xt = 1 | Xt-1 = 1 } = 1 - p10. (2.5)
Figure 1 First-order, Two-state Markov Chain
State 0 (no rain)
State 1 (rain) p
00p
11p
01p
10For example, p10 is the probability that it does not rain on day t but it rains on day t-1.
It should be noticed that in (2.4) and (2.5) p01 = 1 - p00 and p11 = 1 - p10 since there are only two states. According to Wilks (2011), the precipitation time series usually result in p01
< p00 and p10 < p11, which means a rainy day is more likely to be followed by a rainy day rather than a sunny day. Similarly, a sunny day is more likely to be followed by a sunny day rather than a rainy day. From now on, we use px0 as the abbreviation for p00 and p10, x
= 0 or 1.
According to the characteristic (1), the time-varying transition probabilities px0 will not be constants during a year, so px0 is approximated by a truncated Fourier series like (2.6), which could capture the seasonal feature.
0( ) = 0+ =1 sin
( )
2365 + cos( )
2365 , (2.6)where cx0, cxi, dxi are coefficients; M is the order of Fourier series; and according to Cao et al. (2004), the maximum order is set as 5. To simplify the model, we assume the transition probabilities stay constant between years,
px0 (t + 365) = px0 (t). (2.7)
which means we only discuss the behavior of precipitation during a year while the varying differences across years are not in our consideration.
For a time-homogeneous Markov chain, the distribution of steady states is an important property that we usually care about. According to Alexandridis and Zapranis (2013), the stationary probability π1 in a first-order, two-state Markov chain, which describes the unconditional probability of precipitation, is given by
1 = 01
1+ 01 11 = 01
01+ 10
. (2.8)
And the stationary probability π0 for state 0 is
0 =1 1 = 10
01+ 10
. (2.9)
For the steady states, both π0 and π1 should be constants. It implies that
0 1 = 10
01 10 = 0
1 01, (2.10)
which means p10 should be a multiple of p01, and seems to be not very reasonable. This strong condition will limit the efficiency of capturing seasonal features, since four transition probabilities will all share one truncated Fourier series. Therefore, we consider not keeping the steady states in this paper.
2.2 Magnitude Model
The magnitude model describes the amount of rain on the condition that it rains on day t, which is estimated by fitting an appropriate distribution to the historical data.
According to Table 2 in subsection 4.1.2, the descriptive statistics of daily precipitation magnitude show that the distribution should be asymmetric, positively skewed and contain only positive values.
In previous literature, many distributions with a nonnegative domain have been proposed, such as the gamma distribution, the exponential distribution, and the mixed exponential distribution. According to Wilks (2011), the gamma distribution is often used to describe the distribution of various atmospheric variables such as precipitation magnitudes and wind speeds. However, the estimators of the gamma distribution can yield poor results from small values of the shape parameter. A special case of the gamma distribution is the exponential distribution, in which α = 1.
In this paper, we prefer applying the mixed exponential distribution, which is also widely applied in hydrology and in weather derivatives pricing, and can better represent extreme events. This distribution is a weighted combination of two exponential distributions, where α is the weight of the first exponential distribution with parameter β1, and 1-α is the weight of the other exponential distribution with parameter β2.
( ) =
1
1 +1
2
2 , 0 1, 0< 1 < 2
( )
( )
.And since the magnitude of precipitation varies seasonally as described in aforementioned characteristics (3) and (4), the parameters of mixed-exponential distribution α, β1, and β2 are also approximated by the truncated Fourier series like (2.6).
2.3 Estimations
2.3.1 Maximum Likelihood Estimation
Since the parameters in occurrence model (p00 and p10) and magnitude model (α, β1, and β2) are all approximated by the truncated Fourier series, we express the parameter set of the precipitation model as (2.11):
θ(t) = { p00(t), p10(t), α(t), β1(t), β2(t)}. (2.11)
Each parameter is specified by truncated Fourier series as (2.12):
( ) = 0+ sin 2365 + cos
=1
( ) ( 2365 , (2.12))
where θj(t) = {p00, p10, α, β1, β2} with j = 0, 1, 2, …, 4 to distinguish different parameters.
The constraints for these parameters are 0 ≤ px0(t) ≤ 1, 0 ≤ α ≤ 1 and 0 < β1 < β2. We then determine the coefficients of all parameters with MLE and provide the log-likelihood functions of both models as follows:
(1) The log-likelihood functions U0 and U1 for the occurrence model:
We assign Φ1 to be a vector whose elements are the coefficients of the Fourier series describing p00(t) and p10(t). The maximum likelihood estimate Φ1 of Φ1 is found by maximizing the likelihood function L(Xt│Φ1 ), or logarithm: U = logL(Xt│Φ1 ).
We assume there are N years of daily precipitation magnitudes in the data. Then we calculate the historical transition probabilities using the following steps:
Step 1: The transition states for each day are marked as in Table 1. For example, if January 1 is rainy, then January 2 is not rainy. In this case, we mark January 2 as “10”.
Step 2: The marks of transition states are divided into 365 groups. Each group represents one calendar day which comprises N observations. Since February 29 is also included, we need to delete related transition states after having marked them.
Step 3: Starting from January 1, the historical counts are calculated as shown in Table 1.
For example, there are n01(t) days changing from state 0 to state 1 in N years on calendar day t. Finally, there are four series of historical counts and each series contains 365 values.
Table 1 Transition States, Marks and Counts
Transition No Rain to No Rain No Rain to Rain Rain to No Rain Rain to Rain
Marks 00 01 10 11
Counts n00(t) n01(t) n10(t) n11(t)
Step 4: The historical transition probabilities are calculated as (2.13) and (2.14). As mentioned above, p01(t) = 1 - p00(t) and p11(t) = 1 - p10(t).
00( ) = 00( )
00( ) + 01( ). (2.13)
10( ) = 10( )
10( ) + 11( ) . (2.14)
The logarithm of the likelihood function is:
U = U0 + U1 . (2.15)
0 = 00( ) log 00( ) + 01( ) log[1 00( )]
365
=1
. (2.16)
1 = 365 10( ) log 10( ) + 11( ) log[1 10( )]
=1
, (2.17)
since 0
10( ) = 0 and 1
00( ) = 0, U can be considered the sum of U0 and U1, which are log-likelihood functions of p00(t) and p10(t). Therefore, we could maximize U0 and U1 separately, but notice that the constraints are 0 ≤ px0(t) ≤ 1.
(2) The log-likelihood functions U2 for the magnitude model:
We assign Φ2 to be a vector whose elements are the coefficients of the Fourier series describing α, β1, and β2. The maximum likelihood estimate Φ2 of Φ2 is found by maximizing U2, the logarithm of likelihood function L(Φ2│Mt ):
2 { >0} log ( 2| )
365
=1
= 1 , (2.18)
where n = (t modulus 365) + 1, since the period is 1 year (365 days). There are N years of data. The function fn (Φ2│Mt ) is the probability density function (PDF) of the mixed- exponential distribution, in which α(t), β1(t), and β2(t) are truncated Fourier series as (2.12), specified by coefficients Φ2. The indicator function 1{Mt > 0} is defined as (2.19), depending on whether it rains or not.
0, Mt = 0, 1{Mt > 0} =⎧
⎨⎩ (2.19) 1, Mt > 0.
2.3.2 Order Determination of Fourier Series
When fitting models, it is possible that we increase the likelihood by adding parameters, but doing so may result in overfitting. The criteria for model selection are to find the best-fitted model with the minimum number of parameters by introducing a penalty term. For both models, the orders of the Fourier series in parameters θj are determined by the AIC and BIC:
AIC = 2k - 2ln(Lmax), (2.20)
BIC = ln(n)k - 2ln(Lmax), (2.21)
where k is the number of estimated coefficients, n is the number of data points (sample size), and Lmax is the maximum value of the likelihood function. The difference between the two criteria is the penalty term. For either the AIC or BIC, the model with the lowest criterion is preferred. According to Stowasser (2012) and Wilks (2011), the BIC is generally more preferable.
3. The Model of Pricing Precipitation Derivatives
3.1 Simulation Methods of Precipitation 3.1.1 Simulation of the Occurrence Model
The occurrence model could be simulated recursively with the help of a uniform
random variable u1(t) ~ U[0,1] and the initial value X0:
0, if u1(t) ≤ px0(t), Xt =⎧
⎨⎩ (3.1) 1, otherwise.
where px0(t) is the estimated transition probability function, and the use of either p00(t) or p10(t) is determined by Xt-1 = 0 or 1. The variable X0 equals to the initial value of X0. 3.1.2 Simulation of the Magnitude Model
Recalling the PDF of the mixed exponential distribution followed by the magnitude model
( ) = ( )
1( ) 1( ) +1 ( )
2( ) 2( ) , (3.2)
where 0 ≤ α(t) ≤ 1 and 0 < β1(t) < β2(t) for all t. The magnitude model could be simulated with two independent uniform random variables u2(t), u3(t) ~ U[0,1], which are also independent of u1(t). The simulation of the magnitude model is
Mt = Mmin - δtln[u3(t)] , (3.3)
where Mmin is the minimum magnitude of precipitation detected as rain (0.01 inch)4 , and δt is given by
β1(t), if u2(t) ≤ α(t), δt = ⎧
⎨⎩ (3.4) β2(t), if u2(t) > α(t).
The uniform random variable u2(t) decides which exponential distribution is chosen.
If simulates enough times, the distributions will greatly mixed. Since the cumulative distribution function (CDF) of the exponential distribution is
4 Rainfall measuring less than 0.01 inch is defined as the trace precipitation.
{ } = ( ) =1= 1 , >0, (3.5) the uniform random variable u3(t) can be viewed as 1 - F(m) with a specific magnitude m.
Therefore, we can reach the simulated magnitude Mt by (3.3).
The simulated rainfall index is calculated as the production of two simulations:
Rt = Xt ∙ Mt . (3.6)
Since the distribution of the rainfall index is usually asymmetric, skewed and heavy- tailed, normal distribution is obviously not suitable. Therefore, the distribution of the simulated monthly rainfall index should be fitted with non-normal distribution, like log- normal, exponential and normal-inverse Gaussian (NIG) distribution by MLE.
3.2 Esscher Transform
Besides the skewness and semi-heavy tails of the rainfall index distribution, the properties of the weather derivatives market also make it hard to price. Traditionally, derivatives are priced with the Black-Scholes model (1973), in which the risk-neutral probability measure ℚ equivalent to the physical measure ℙ is specified, such that all discounted tradable assets in the market are martingales under the ℚ measure. This unique risk-neutral measure ℚ is obtained by the Girsanov transform (also known as the Girsanov theorem), changing the drift term of the Brownian motion.
But the general consensus is that weather is not a tradable asset, which means the markets for precipitation derivatives, even for weather derivatives, are inherently incomplete. It is impossible to construct a riskless hedge portfolio containing such derivatives. In this case, we cannot find a unique risk-neutral measure ℚ equivalent to the physical measure ℙ. Instead, there are many equivalent martingale measures.
Hence, we choose to apply the Esscher transform to price the precipitation derivatives. According to Gerber and Shiu (1994), the Esscher transform provides a class of risk-neutral measures, which can be justified by assuming a representative investor who wants to maximize the expected utility with the market price of risk (MPR) parameterized by θ. The MPR is defined as a measure of the extra return or risk premium that investors demand to bear risk, which also reflects the risk attitude of market participants. The Esscher transform changes a probability density f(x) of a random variable X into a new
probability density f(x;θ) with MPR parameter θ.
( ; ) = ( )
( ) . (3.7)
The Esscher transform can also be extended to stochastic processes. We make Xt, 0 ≤ t ≤ T a Lévy process with the Lévy-Khintchine triplet (σ2, υ, γ) and assume the Lévy measure υ is such that | | 1 ( ) < . According to Cont and Tankov (2004), the Radon-Nikodym derivative corresponding to the Esscher transform is
|
| = [ ] = ( ) , (3.8)
where κ(θ) = ln E[eθX] is the logarithm of the moment generating function of X.
A big advantage of Esscher transform is the great applicability to any distribution’s probability density function f (x). In this paper, we choose the normal-inverse Gaussian distribution, whose probability density function NIG(μ, α, β, σ) remains in the original form after the Esscher transform, as NIG(μ, α, β + θ, σ).
3.3 Pricing Formulas
The subject for discussion in this paper is CME rainfall index binary contracts. And we define the underlying rainfall index (e.g., monthly rainfall index) I(τ1, τ2) as
( 1, 2) =
2
=1
, [ 1, 2] , (3.9)
which is the sum of the daily rainfall index Rt, and [τ1, τ2] is the particular accumulation period. The payoff of a rainfall index binary call with strike price K is
10,000, I(τ1, τ2) > K.
C(τ2, τ1, τ2)=⎧
⎨⎩ (3.10) 0, otherwise.
According to Table A.1 in Appendix A, the price of CME rainfall index binary
contracts is quoted in terms of the respective CME rainfall index, in which each point represents $100. Under the Esscher transform measure ℚθ, the expected rainfall index (ERI) at time t is a conditional risk-adjusted expectation adapted to the filtration ℱt
ERI(t, τ1, τ2) = Eℚθ ) [I(τ1, τ2) | ℱt]. (3.11)
the risk-neutral pricing formulas of binary call C(t, τ1, τ2) and put P(t, τ1, τ2) are
C(t, τ1, τ2) = 100Eℚθ [I(τ1, τ2) 1{I(τ1, τ2) > K}| ℱt] - K. (3.12) P(t, τ1, τ2) = K -100Eℚθ [I(τ1, τ2) 1{I(τ1, τ2) < K}| ℱt]. (3.13)
4. Empirical Analysis
4.1 Data
4.1.1 Precipitation Derivatives Market Data
Table A.1 lists the specifications of the CME Rainfall Index Binary Contracts. We use the daily prices of these contracts (the rainfall binary monthly options), whose market data is obtained from the CME. The contracts are for March 2005 and 2006, which were traded in the OTC market and cleared by CME ClearPort. And the reference station is Des Moines International Airport, U.S. The details are shown in Table A.2 in Appendix A.
4.1.2 Precipitation Data
The rainfall data we use is the daily rainfall (in inches) at the Des Moines International Airport, U.S., from January 1, 1946 to December 31, 2004 (59 years in total), provided by the National Oceanic and Atmospheric Administration (NOAA). Table 2 presents the descriptive statistics, in which the kurtosis and skewness of non-zero data show that the daily rainfall is positively skewed and has a fat tail.
Table 2 Descriptive Statistics of Historical Daily Rainfall Data
Data Count Mean Std Kurtosis Skewness Min Max
Zero 15,164
Non-zero 6,386 0.3003 0.4446 16.9431 3.2821 0.01 6.18
In order to justify the rationality of the seasonal characteristic assumption in transition probability, we divide 59 years of data into three groups with each consisting of two decades (except the last group which consists of 19 years of data), and by taking p00 as an example, we calculate the historical daily transition probability in a year (365 days). After comparing the annual trends in Figure 2, we think the seasonal assumption of transition probability is reasonable.
Besides, we divide the data into six groups with each consisting of a decade (except the last group which consists of 9 years of data), and calculate the mean and median magnitude of monthly precipitation for each group. According to the trend observed in Figure 2, it is reasonable to make the seasonal assumption of the magnitude parameters.
Figure 2 The Test of the Seasonal Characteristic Assumptions
4.2 Estimations of Precipitation Models 4.2.1 Estimation of Occurrence Model
We fit the transition probabilities of the occurrence model with the truncated Fourier series from first-order (M = 1) to fifth-order (M = 5), and estimate the coefficients by MLE.
Since MLE is a numerical optimization problem, we should notice that the coefficients need to satisfy the constraint (4.1) for all transition probabilities.
0 ( ) = 0+ 2
365 +
2
365 1. (4.1)
When we apply MLE in numerical computation, the initial values of coefficients are obtained by the least-squares approach. However, according to Woolhiser and Pegram (1979), the data points (transition probabilities) are estimates and have unequal variances because of varying properties of the distribution. Due to this major deficiency, we do not directly adopt the least-squares approach to estimate the coefficients.
Table 3 shows the results estimated by MLE and the values of the AIC and BIC, which are used to decide the order of each transition probability. We find that the criteria of p00 are close between the 1st order and 2nd order in the AIC but in the BIC, these two orders are obviously different. According to the lowest criterion, we choose the 1st order for p00 and the 2nd order for p10. Figure 3 shows the estimated truncated Fourier series of p00 and p10, and all orders are shown in Figure B.2 of Appendix B.
4.2.2 Estimation of Magnitude Model
We fit the daily rainfall data with the mixed exponential distribution, where the parameters α, β1, and β2 are described with the truncated Fourier series from first-order (M = 1) to fifth-order (M = 5). The coefficients of each parameter are estimated by MLE with constraints 0 ≤ α ≤ 1, 0 < β1 < β2. According to the estimated results and the values of the AIC and BIC in Table 4, we choose the 3rd order for α, β1, and β2. Figure 4 shows the estimated truncated Fourier series, and all orders are shown in Figure B.1 of Appendix B.
Table 3 The Coefficients of the Truncated Fourier Series in the Occurrence Model Estimated by MLE and the Values of the AIC and BIC Orderθjcj0cj1dj1cj2dj2cj3dj3cj4dj4cj5dj5AICBICU 1p000.7613 0.7613 -0.0314 0.0442 16451.316463.0-8222 p100.5644 -0.0302 0.0185 8741.18752.8-4367 2p000.7612 -0.0320 0.0452 0.0094 -0.0027 16451.216470.7-8220 p100.5652 -0.0274 0.0184 0.0254 0.0175 8732.98752.4-4361 3p000.7612 -0.0322 0.0451 0.0098 -0.0020 0.0007 -0.0067 16453.216480.5-8219 p100.5651 -0.0273 0.0184 0.0255 0.0184 -0.0047 0.0076 8735.88763.1-4360 4p000.7611 -0.0322 0.0451 0.0099 -0.0021 0.0000 -0.0065 0.0052 0.0019 16455.816490.9-8218 p100.5651 -0.0274 0.0183 0.0253 0.0183 -0.0036 0.0079 0.0036 0.0101 8735.88763.1-4360 5p000.7610 -0.0320 0.0451 0.0097 -0.0020 0.0001 -0.0063 0.0056 0.0011 -0.0061 0.0033 16457.616500.5-8217 p100.5651 -0.0274 0.0183 0.0253 0.0184 -0.0036 0.0080 0.0036 0.0101 0.0002 -0.0004 8742.38785.2-4360 Figure 3 Estimated Fourier Series of Transition Probabilities
Table 4 The Coefficients of the Truncated Fourier Series in the Magnitude Model Estimated by MLE and the Values of the AIC and BIC Orderθjcj0cj1dj1cj2dj2cj3dj3cj4dj4cj5dj5AICBICU 1α0.7485 0.2511 0.0059 β10.2816 -0.0953 -0.1462 β20.6861 0.5602 -0.3961 -4371.2-4299.42195 2α0.3479 -0.1130 -0.2141 0.0737 -0.0660 β10.3108 0.3121 0.1723 0.1926 -0.0813 β20.4140-0.2589 -0.3418 0.0774 0.0331 -5225.8-5106.12628 3α0.5814-0.0846 -0.0532 -0.0248 0.0187 -0.1370-0.0355 β10.4851 -0.0428 -0.1844 -0.0339 -0.0692 0.1468 0.0182 β20.0486 0.0133 -0.0253 -0.0072 -0.0041 0.0229 0.0131 -5959.6-5792.73001 4α0.4104 0.0373 -0.0140 -0.0157 -0.0358 0.0370 -0.0433 -0.0424 -0.0827 β10.1624 -0.1139 -0.2078 0.1468 0.0645 -0.1004 0.0277 0.0295 -0.0468 β20.4023 0.0156 -0.1002 -0.1459 -0.0755 0.1235 -0.0518 -0.0296 -0.0159 -5506.1-5290.72780 5α0.4089 -0.0581 0.0165 0.1007 0.0445 0.0460 0.0199 -0.0049 -0.0319 0.0288 -0.0095 β10.2685 -0.0843 -0.1311 0.1333 0.0579 -0.1079 0.0795 -0.0316 -0.0965 -0.1221 0.0671 β20.2869 0.0042 -0.0913 -0.1114 -0.0515 0.1045 -0.0363 0.00910.0430 0.0495 -0.0333 -4928.2-4665.02497 Figure 4 Estimated Truncated Fourier Series of Parameters α, β1 and β2
4.3 Simulation and Estimation of Index Distribution
Since the rainfall of December 31, 2004 and December 31, 2005 is both 0 inches, the path of the rainfall index is simulated for 10,000 times with the same initial value (0 inches), which implies that the distribution of these two years are the same. Simulated monthly rainfall indexes from March to October are calculated and fitted with normal, exponential, and NIG distribution by MLE. The histograms and fitted PDF are presented in Figure C.1 of Appendix C, which shows that the NIG distribution fits best. We draw the QQ plot with the NIG distribution for each month in Figure C.2 of Appendix C. We also use the one-sample Kolmogorov-Smirnov test (K-S test) to check the equality between the simulated sample and the fitted NIG distribution. The result is shown in Table 5.
Table 5 Result of One-Sample Kolmogorov-Smirnov Test
Month Mar. Apr. May Jun. Jul. Aug. Sept. Oct.
D 0.0122 0.0118 0.0129 0.0143 0.0128 0.0166 0.0190 0.0184
P-value 0.3328 0.3748 0.2713 0.1731 0.2776 0.0738 0.0263 0.0345
According to the definition, there is no significant difference between the sample and the fitted distribution when the p-value of the K-S test is greater than 0.05. In addition, it has a better performance from March to August, but does not fit so well in September and October. The fitted parameters of the NIG distribution are shown in Table 6.
Table 6 Fitted Parameters of NIG Distribution (Mar. to Oct.)
Param Mar. Apr. May Jun. Jul. Aug. Sept. Oct.
α 79.5702 49.3685 37.6038 33.8541 16.5078 43.8524 41.1671 42.1280 β 78.2652 48.4672 37.0062 33.2425 15.6259 43.0413 40.6162 41.5929 μ -1.2892 -1.7029 -2.0678 -2.2737 -2.1244 -2.2339 -1.8373 -1.4369
δ 0.6200 0.9637 1.1913 1.3258 2.0743 1.1489 0.8826 0.7375
4.4 Theoretical Prices and Risk Management 4.4.1 Market Price of Risk
Since we fit the simulated monthly rainfall index with the NIG distribution, by taking the rainfall index binary call option as an example, the market price of risk (MPR) θ under the Esscher transform can be derived as (4.2), which is proved in Appendix C.
=
+100
2+ +
100
2 . (4.2)
With the parameters estimated in Table 6, we calculate the theoretical expected rainfall index for different θ from 0.5 to -0.5, which is presented in Table 7.
4.4.2 Burn Analysis and Model Performance
According to Dorfleitner and Wimmer (2010) and Cabrera et al. (2013), as an actuarial method, burn analysis (BA) is widely used in pricing under physical measure ℙ, which values a derivative by averaging all payoffs that have been realized in the past.
Since weather derivatives are considered a substitute for weather insurances, the BA is widely applied as a benchmark approach in practice.
We calculate the BA with the historical monthly rainfall index from January 1946 to December 2004. It is noteworthy that the BA cannot be considered equivalent to θ = 0 since BA is under physical measure ℙ and θ = 0 is under the risk-neutral measure ℚθ. The pricing results (θ = 0) are compared with the BA and the real monthly rainfall index of Des Moines, IA, USA in 2005, which are summarized in Table 7.
Table 7 Burn Analysis, Real Rainfall Index in 2005 and Theoretical MPR
MPR Rainfall Index
θ Mar. Apr. May Jun. Jul. Aug. Sept. Oct.
2005 1.6100 4.7300 5.1600 4.6800 3.2800 1.3300 1.4100 0.9200
BA 2.2347 3.2903 4.2604 4.4301 3.7308 3.9024 2.9459 2.3934
(0.6247) (1.4397) (0.8996) (0.2499) (0.4508) (2.5724) (1.5359) (1.4734) -0.50 1.5719 2.2551 2.7541 2.7720 2.6212 2.3586 1.9935 1.8283 -0.40 1.6575 2.4108 3.0005 3.0277 2.8275 2.5529 2.1975 2.0048 -0.30 1.7508 2.5853 3.2863 3.3237 3.0585 2.7725 2.4360 2.2117 -0.20 1.8532 2.7828 3.6237 3.6721 3.3195 3.0235 2.7201 2.4590 -0.10 1.9663 3.0088 4.0304 4.0906 3.6181 3.3142 3.0668 2.7622 0.00 2.0920 3.2710 4.5343 4.6069 3.9645 3.6566 3.5033 3.1460 (0.4820) (1.4590) (0.6257) (0.0731) (0.6845) (2.3266) (2.0933) (2.2260) 0.10 2.2329 3.5804 5.1820 5.2664 4.3734 4.0681 4.0768 3.6547 0.20 2.3923 3.9533 6.0589 6.1516 4.8670 4.5761 4.8791 4.3753 0.30 2.5748 4.4149 7.3444 7.4315 5.4801 5.2255 6.1201 5.5147 0.40 2.7866 5.0076 9.5059 9.5303 6.2715 6.0979 8.4422 7.7498 0.50 3.0365 5.8087 14.4317 14.0177 7.3513 7.3604 15.8938 16.6189
The parentheses contain the absolute difference between each approach and the real rainfall index in 2005. The results of BA are around the expected rainfall index of θ = 0, ranging from θ = 0.1 to θ = -0.1. Compared with the BA, the benchmark approach, the model we use has a better performance in the first half of the year.
4.4.3 Sensitivity Analysis
Since the quoted theoretical price can be expressed in a linear form with expected rainfall index Eℚθ[I(τ1, τ2) | ℱt] and strike price K as shown in (3.11) and (3.12), we then discuss the relationship between the MPR (θ) and ERI from two aspects.
(1) The Influence of the NIG Distribution
Comparing the expectation before and after the Esscher transform, the bigger θ means a larger β, which positively measures the asymmetry degree of the NIG distribution. It also influences the expectation of the distribution, since a bigger θ means that the investors have a larger estimation of the monthly rainfall index.
(2) The Influence of the Utility Function
In the Esscher transform, we assume that all market participants are represented by a representative agent with exponential utility function u(x) = (1- e-θx) ⁄ θ, in which the parameter θ is the MPR. The Arrow-Pratt measure of Absolute Risk Aversion (ARA) represents the risk attitude of market participants as (4.3),
( ) = ( )
( ) = = , (4.3)
which means it is a Constant Absolute Risk Aversion (CARA) under the exponential utility function. When θ > 0, market participants are more risk averse as opposed to risk loving when θ < 0.
By taking the data of March as an example, Figure 5 shows the monthly ERI is positively related to the MPR θ, which means market participants are more risk averse with greater ERI.
Figure 5 Monthly ERI with Varieties of the MPR in Mar.
4.4.4 Market Data Analysis
For contracts in Table A.2 of Appendix A, we calculate the MPR from January 3 to March 31 in 2005 and 2006 with market data from the CME and analyze the temporal behavior.
Figure 6 The MPR of Binary Calls and Binary Puts with Different Strike Prices in 2005 and 2006
