# Modeling Drought Option Contracts

Authors: Jielin Zhu; Marco Pollanen; Kenzu Abdella; Bruce Cater
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

We introduce a new financial weather derivative—a drought option contract—designed
to protect agricultural producers from potential income loss due to agricultural drought. The contract
is based on an index that reflects the severity of drought over a long period. By modeling temperature
and precipitation, we price a hypothetical drought contract based on data from the Jinan climate station
located in a dry region of China.


---

## Body

## 1. Introduction

Since the dawn of agriculture, farmers have been vulnerable to the effects of drought. These effects include not only any immediate financial losses that result from reduced crop production, but also the adverse effects that those losses may have on their operational viability going forward.Historically, the only financial instruments available to mitigate the risks associated with drought were crop insurance contracts. Those contracts have the obvious advantage of minimizing the financial losses incurred by farmers, thus stabilizing their operations. But because the marginal benefits of farmers’ crop damage-mitigating strategies are offset by reductions in the expected value of the insurance claim, the private benefits they realize from those efforts will fall short of the social benefits. From a social perspective, the resulting moral hazard then leads to an inefficiently low level of damage-mitigating effort and an inefficiently high level of realized crop loss. Moreover, because of drought’s high positive spatial correlation, regional insurers may themselves have difficulty meeting their obligations in the event of drought, thus undermining the ability of insurance contracts to reduce the drought risks farmers face.By the mid-1990s, however, new financial contracts known as weather derivatives had emerged as a way of reducing particular weather-related risks, including those resulting from uncertain temperature, precipitation, and frost. Because the payoffs of these derivatives depend only on weather variables over which the farmer has no control, these instruments have the advantage of sustaining the damage-mitigation incentives that crop insurance contracts may undermine, thus mitigating real crop losses [1]. Moreover, relative to an insurance market, a robust weather derivative market may better distribute correlated risks.Weather derivatives are very similar to traditional derivatives on financial assets. But because the weather indices they depend on have no underlying measurable value [2–5], there is no underlying instrument with which to hedge, and the Black-Scholes sort of risk-neutral pricing models cannot be applied. As such, much of the literature on weather derivatives—including Garman et al. [6], Brody et al. [7], Richards et al. [8], Cao and Wei [4], Geman and Leonardi [9], and Taylor and Buizza [10]—has focused on the development of alternative models to price weather risk. Because the weather variables on which those models are based are likely to be imperfectly correlated with their impact on the value of the crop yield, however, a basis risk remains.In this paper, we propose, model, and price a weather derivative, which we call a “drought option,” that is based not only on multiple weather factors (i.e., precipitation and temperature), but also on an agricultural drought index that depends on the particular crop growth being hedged as a way of reducing this basis risk.The balance of the paper is organized as follows. In Section2, we describe the structure of the drought option contract as well as the drought index upon which the contract depends. In Section 3, we introduce three alternative approaches to pricing the option. In Section 4, we use data from Jinan climate station in eastern China to simulate the option prices for each of our three approaches. Section 5 concludes with comparisons of the three approaches as well as a discussion of possible directions for future work.


## 2. Background

### 2.1. Drought Options

We model our derivative contract as a European put option on an agricultural drought index,I, with strike level K, and tip S, which reflects the relationship between the drought index and the loss of farmers’ income. At maturity, the payoff function is(2.1)f(I)=S⋅max⁡(K-I,0).
If I—an index that takes on smaller values when a drought is more severe—is less than K, the holder will receive S(K-I); if I is larger than K, the holder would choose not to exercise the contract and the profit would be 0.Because no-arbitrage pricing methods are unsuitable, we use the equilibrium approach and price the option as the present value of the expectation of its payoff [4, 11]. That expectation can be written as(2.2)V=EQ[e-∫0τr(t)dtf(I)],
where τ is the contract’s duration, r(t) is the risk-free interest rate, I is the value of drought index, and ℚ is the probability density over all possible values of I.To simplify the problem, we will assume thatS=1 and that r(t) is constant over the length of the contract, τ. Equation (2.2) can then be rewritten as(2.3)V=e-rτEQ[max⁡(K-I,0)].
Our problem is thus reduced to choosing an appropriate drought index, setting the strike level, and estimating the expectation.


### 2.2. Drought Index

Because our interest lies in developing a tool to help farmers hedge against the risks of reduced agricultural yield, our approach will utilize the Reconnaissance Drought Index (RDI) [12, 13], which measures the severity of drought in agricultural terms. For any period within year i, the RDI can be calculated as(2.4)RDIi=∑j=1mPij∑j=1mPETij,i=1,…,N,j=1,…,m,
where Pij and PETij are the precipitation and potential evapotranspiration, respectively, for the jth month of the ith year, m is the number of observed months in year i, and N is the number of observed years (Evapotranspiration (ET) is a measure of water movement from the Earth’s surface to the atmosphere; it is the sum of evaporation, such as that from soil and canopy, and plant transpiration. Potential Evapotranspiration (PET) is the evapotranspiration that could occur where there is sufficient water supply. Given the relationship between P and PET, the RDI takes on smaller values when the drought is more severe, making our choice of a put option appropriate).


### 2.3. Potential Evapotranspiration and Actual Evapotranspiration

Potential evapotranspiration is difficult to measure because it depends on the type of plants, the type of soil, and the climatic conditions. For practical purposes, we will therefore use theactual rather than the potential evapotranspiration in our estimation of the RDI.This gives us a newAdjusted-RDI, that can be written as(2.5)RDIadi=∑j=1mPij∑j=1mETij,i=1,…,N,j=1,…,m,
where ETij is the actual evapotranspiration in the jth month of the ith year.To estimate theET, we will use the temperature-based methods of Blaney-Criddle [14] (Methods of estimating evapotranspiration can be divided into three categories [15]: (1) mass-transfer-based methods; (2) radiation-based methods; (3) temperature-based methods. Our use of a temperature-based method is necessitated by data limitations). The formula for the SCS Blaney-Criddle method is(2.6)ET=kcktf,
where(2.7)kt=0.0311T+0.24,f=d(1.8T+32)100,
and where the monthly ET measured in inches, kc is the consumptive crop coefficient for the SCS version, T is monthly mean temperature, and d is the percentage of daylight hours.


## 2.1. Drought Options

We model our derivative contract as a European put option on an agricultural drought index,I, with strike level K, and tip S, which reflects the relationship between the drought index and the loss of farmers’ income. At maturity, the payoff function is(2.1)f(I)=S⋅max⁡(K-I,0).
If I—an index that takes on smaller values when a drought is more severe—is less than K, the holder will receive S(K-I); if I is larger than K, the holder would choose not to exercise the contract and the profit would be 0.Because no-arbitrage pricing methods are unsuitable, we use the equilibrium approach and price the option as the present value of the expectation of its payoff [4, 11]. That expectation can be written as(2.2)V=EQ[e-∫0τr(t)dtf(I)],
where τ is the contract’s duration, r(t) is the risk-free interest rate, I is the value of drought index, and ℚ is the probability density over all possible values of I.To simplify the problem, we will assume thatS=1 and that r(t) is constant over the length of the contract, τ. Equation (2.2) can then be rewritten as(2.3)V=e-rτEQ[max⁡(K-I,0)].
Our problem is thus reduced to choosing an appropriate drought index, setting the strike level, and estimating the expectation.


## 2.2. Drought Index

Because our interest lies in developing a tool to help farmers hedge against the risks of reduced agricultural yield, our approach will utilize the Reconnaissance Drought Index (RDI) [12, 13], which measures the severity of drought in agricultural terms. For any period within year i, the RDI can be calculated as(2.4)RDIi=∑j=1mPij∑j=1mPETij,i=1,…,N,j=1,…,m,
where Pij and PETij are the precipitation and potential evapotranspiration, respectively, for the jth month of the ith year, m is the number of observed months in year i, and N is the number of observed years (Evapotranspiration (ET) is a measure of water movement from the Earth’s surface to the atmosphere; it is the sum of evaporation, such as that from soil and canopy, and plant transpiration. Potential Evapotranspiration (PET) is the evapotranspiration that could occur where there is sufficient water supply. Given the relationship between P and PET, the RDI takes on smaller values when the drought is more severe, making our choice of a put option appropriate).


## 2.3. Potential Evapotranspiration and Actual Evapotranspiration

Potential evapotranspiration is difficult to measure because it depends on the type of plants, the type of soil, and the climatic conditions. For practical purposes, we will therefore use theactual rather than the potential evapotranspiration in our estimation of the RDI.This gives us a newAdjusted-RDI, that can be written as(2.5)RDIadi=∑j=1mPij∑j=1mETij,i=1,…,N,j=1,…,m,
where ETij is the actual evapotranspiration in the jth month of the ith year.To estimate theET, we will use the temperature-based methods of Blaney-Criddle [14] (Methods of estimating evapotranspiration can be divided into three categories [15]: (1) mass-transfer-based methods; (2) radiation-based methods; (3) temperature-based methods. Our use of a temperature-based method is necessitated by data limitations). The formula for the SCS Blaney-Criddle method is(2.6)ET=kcktf,
where(2.7)kt=0.0311T+0.24,f=d(1.8T+32)100,
and where the monthly ET measured in inches, kc is the consumptive crop coefficient for the SCS version, T is monthly mean temperature, and d is the percentage of daylight hours.


## 3. Drought Option Pricing Models

We now wish to calculate possible values of ourAdjusted-RDI in order to estimate the option price. Common techniques to calculate index values include Historical Burn Analysis, Index Value Simulation, and Stochastic Simulation [11].


### 3.1. Historical Burn Analysis and Index Value Simulation

For Historical Burn Analysis, the basic assumption is that the historical data reasonably approximate future scenarios, allowing the price of the option to be calculated from data we already have. From (2.3), we can see that the price of a put option for drought is given by(3.1)V=1Ne-rτ(∑i=1Nmax⁡(K-Ii,0)),
where Ii is the value of the drought index for year i and N is the number of years of historical data.Our Index Value Simulation approach involves finding the best fit distribution for index values calculated using the Historical Burn Analysis approach and then sampling this distribution to produce possible future values of the index. The option price will then be calculated using (3.1).


### 3.2. Stochastic Simulation

Mean reversion stochastic processes are used to simulate the daily mean temperature and the speed of monthly rainfall over the entire year. The period’s rainfall and evapotranspiration are then obtained by summing the related simulated values. Finally, the option payoff is computed based on the simulation where the average discounted payoff is used as the option price.


#### 3.2.1. Simulation of Mean Monthly Temperature

Although based on (2.6) and (2.7), only the mean monthly temperatures are needed to calculate the evapotranspiration, and those mean temperatures can be simulated directly.Mean Reversion Process
Average daily temperature is, of course, subject to seasonal changes. For one specific day in each year, however, we will assume that the temperature fluctuates around some mean value, which means we can choose a mean reversion stochastic process to simulate the behavior of daily temperature [2, 16].
We use the following Stochastic Differential Equation (SDE) to model daily temperature:(3.2)dTt=dTtm+a(Ttm-Tt)dt+σtdWt.
The solution is
(3.3)Tt=(Ts-Tsm)e-a(t-s)+Ttm+∫ste-a(t-τ)στdWτ,
where Tt is the temperature at time t, Ttm is the mean daily temperature at time t, a reflects the speed of the mean reversion, σt reflects the variation from the mean value at time t, and Wτ is a Wiener process.Parameter Estimation
Ttm in (3.2) represents simulated mean daily temperature. If the mean temperature is assumed to be a continuous variable then it is possible to simulate this curve by a sinusoidal function of the form:
(3.4)sin⁡(ωt+φ),
where t denotes the time, measured in days.
To account for local warming trends observed in the data,Ttm is modeled as having the form:(3.5)Ttm=A+Bt+Csin⁡(ωt+φ).
Expanding (3.5) gives us
(3.6)Ttm=a1+a2t+a3sin⁡(ωt)+a4cos⁡(ωt).
We then estimate the parameters in (3.6) by using the Gauss-Newton algorithm to solve(3.7)min⁡ξ⟦Tt-X2⟧,
where X is the data vector and ξ is the parameter vector (a1,a2,a3,a4).
Estimates ofξ then also allows us to calculate the parameters in (3.5) as follows:(3.8)A=a1,B=a2,C=a32+a42,φ=arctan(a4a3)-π.
The termσtdWt in (3.2) reflects the discrepancy of real temperature from the simulated mean temperature. We simplify the function σt as a piecewise constant function. Here we assume the value of σt is a constant number for each month.
As we see,σt reflects the variation of daily temperature, the first estimator is based on the quadratic variation of Tt:
(3.9)σμ2=1Nμ∑j=0Nμ-1(Tj+1-Tj)2,
where Nμ is the number of days in month μ.
If we discretize (3.2), for a given month μ, we have(3.10)Tj=Tjm-Tj-1m+aTj-1m+(1-a)Tj-1+σμϵj-1,
with j=1,…,Nμ, where {ϵj}j=1Nμ-1 follow the standard normal distribution. If we let T̅j be Tj-(Tjm-Tj-1m), then (3.10) becomes
(3.11)T̅j=aTj-1m+(1-a)Tj-1+σμϵj-1,
which we can see as a regression of today’s temperature on yesterday’s temperature. Therefore, the second estimator of σμ is
(3.12)σ̂μ2=1Nμ-2∑j=2Nμ(T̅j-âTj-1m-(1-â)Tj-1)2,
where â is an estimator of a.
To estimate the mean reversion parametera, we use the martingale estimation functions method. After collecting observations during n days, an efficient estimator ân of a can be obtained by solving the equation(3.13)Gn(ân)=0,
where
(3.14)Gn(a)=∑i=1nḃ(Ti-1;a)σi-12(Ti-E[Ti∣Ti-1]),
and where ḃ(Ti-1;a) denotes the derivative with respect to a of the term
(3.15)b(Tt;a)=dTtmdt+a(Ttm-Tt).
Therefore, we have
(3.16)ḃ(Ti-1;a)=Ttm-Tt.
To then solve (3.14), all we need is the formula Ti-E[Ti∣Ti-1]. As we know the solution of (3.2) is (3.3), and the expectation of the integration part is 0, when t=i and s=i-1, we have(3.17)E[Ti∣Ti-1]=(Ti-1-Ti-1m)e-a+Tim.
So we have
(3.18)Gn(a)=∑i=1nTi-1m-Ti-1σi-12[Ti-(Ti-1-Ti-1m)e-a-Tim].
It is easy to estimate the parameter a from
(3.19)ân=-log⁡(∑i=1n-1Yi-1(Ti-Tim)∑i=1nYi-1(Ti-1-Ti-1m)),
where
(3.20)Yi-1≡Ti-1m-Ti-1σi-12,i=1,2,…n.
Based on (3.15), we can see that ân is the unique estimation.Numerical Solution
Because it is difficult to use the integral solution in (3.3), we employ the implicit Milstein numerical method [17, 18]. The problem is then reduced to the following iterative procedure:(3.21)Tn+1(0)=Tn+(Tnm′+a(Tnm-Tn))Δt+σkΔW,Tn+1(1)=Tn+(Tn+1m′+a(Tn+1m-Tn+1(0)))Δt+σkΔW,
where Tn+1(1) is the simulated point we get, σk is the value of diffusion parameter depending on the nth point located in month k.


#### 3.2.2. Simulation of the Speed of Monthly Rain

Transforming Monthly Rain to a Continuous Process
While precipitation does not behave continuously, the speed of precipitation can be assumed to change continuously in time. When rain starts, the speed of precipitation increases continuously from zero to its peak and decrease continuously to zero where it remains at zero until it rains again. In this formulation, the total amount of precipitation can be obtained by integrating the speed of precipitation over the period of interest, as given by the formula:(3.22)Pt1,t2=∫t1t2x(t)dt,
where x(t) represents the speed of rainfall at time t and Pt1,t2 represents the amount of precipitation during the period from time t1 to time t2.
The monthly rainfall, then, is essentially the monthly mean speed of precipitation with units (mm/m2·month). In order to simulate monthly rainfall, it will be sufficient to simulate the speed of rainfall using a mean reversion process.Mean Reversion Process
In the previous model of monthly temperature, diffusion was modeled as a piecewise functionσt, where σt was assumed constant for each month. In order to make the simulation of the speed of monthly rain more realistic, the diffusion can be modeled as a function of Xt according to the following SDE:(3.23)dXt=dθ(t)+a(θ(t)-Xt)dt+σXtpdWt,
with t≥0, a≥0, and θ(t)>0. As described in (3.2), θ(t) is the simulated mean speed of monthly rain at time t over a year and a is the mean reversion parameter. Here σ and p together reflect the volatility which is now dependent on previous precipitation.Parameter Estimation
The data we will utilize are for the Jinan climate station in eastern China, obtained from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/) for the period from January 1st, 1951 to December 31st, 2006.
Figure1 depicts the mean monthly rain speed over the 56-year period. As in the temperature model, the real mean speed of monthly precipitation will be assumed to fit a sinusoidal function of the form [19]:(3.24)θ(t)=m+αsin⁡(2π(t-υ)12),
where m is the mean of the sine curve, α determines the oscillation, and υ is the horizontal shift. To make the simulated curve closer to the historical experience, the simulation can be further improved by expanding θ in terms of a Fourier series [19]:
(3.25)θ(t)=m+∑i=0nαisin⁡((2i+1)π(t-υ)6).
The parameters inθ(t) can then be estimated using the Gauss-Newton method to solve the least squares problem. To do this, we choose n=1.
We now wish to find an unbiased estimator of the mean reversion parametera. The original estimation of a is [19]:
(3.26)â=1n∑i=0n-1Xi+1-Xi-θ(i+1)+θ(i)(θ(i)-Xi)Δ,
where Δ is the time increment. Because this estimation would change due to different lengths of data, we have the modification:
(3.27)ab̂=1|Ib|∑i∈IbXi+1-Xi-θ(i+1)+θ(i)(θ(i)-Xi)Δ,
with Ib={i=1,…,n:|θ(i)-Xi|>b}, b a positive real number. As b gets larger, ab̂ is convergent. Here we choose b=50 for a suitable estimation.
The diffusion parameters can be estimated by first squaring equation (3.23) to get(3.28)(dXt)2=a2(θ-Xt)2⋅dt⋅dt+a⋅(θ-Xt)⋅σ⋅Xtp⋅dt⋅dWt.
Then, based on Itô’s integration rule, we have
(3.29)ln⁡(dXt)2=2ln⁡(σ)+ln⁡(dt)+2pln⁡(Xt),
such that ln(dXt)2 and ln(dXt) have a linear relationship. The data can now be substituted into (3.29). We use linear regression method to estimate σ and p.Figure 1
Mean speed of monthly rain over 56 years.Numerical Solution
After the parameter estimation, based on the implicit Milstein method, the monthly precipitation can be simulated by the following:(3.30)Xn+1(0)=Xn+(θ′(n)+a(θ(n)-Xn(0)))Δtn+σXn(1)pΔWn+12σ2pXn(1)2p-1(ΔWn2-Δtn),Xn+1(1)=Xn+(θ′(n+1)+a(θ(n+1)-Xn+1(0)))Δtn+σXn(1)pΔWn+12σ2pXn(1)2p-1(ΔWn2-Δtn),
where Xn+1(1) is the simulated point at the (n+1)th step, n=0,…,N and N is the step size for one simulation process.


### 3.3. Correction of the Model

#### 3.3.1. Correction for the Mean Curve

For both temperature and precipitation models, the distance between the average of historical records and the simulated mean curve can also be simulated with another sinusoidal function like (3.25), with the exception that the period of this function may be different depending on the climate conditions of the location.Figures2 and 3 compare the simulated mean curve with and without the simulation of distance. It is clear that the inclusion of simulated distance provides results that are much closer to the historical record than the original one.Figure 2
Simulated daily temperature comparison.Figure 3
Simulated mean speed of rain comparison.


#### 3.3.2. Maintaining the Positivity of the Rain Speed

SinceXt represents the speed function, it must be nonnegative to have a meaningful interpretation. Moreover, according to (3.23), negative Xt leads to no value for Xt2p-1. However, there is no guarantee that the simulation of (3.23) yields positive values, especially when Xt is close to 0 at the beginning of the year [17]. If we examine Figure 3, we see that the mean speed of rain is under the time axis at the beginning of the year.IfXt represents a stochastic process with(3.31)Prob({Xt>0∀t})=1,
the stochastic integration scheme possesses an eternal life time if(3.32)Prob({Xn+1>0∣Xn>0})=1.
Otherwise, we say it has a finite lifetime.If we can find a numerical method to solve this model which has an eternal life, as long as the initial value of historical precipitation is positive, we can make sure that all the following simulated points are positive. This is another reason to choose the implicit Milstein method to do the simulation [19].For the mean reversion process given by(3.33)dXt=(α-βXt)dt+σXtpdW,
with α,β,σ,p∈R+ and p>1/2, the Milstein method provides numerical positivity with the following restriction:(3.34)Δt<1σ2.Compared to (3.23), (3.33) will guarantee positivity with(3.35)α=a⋅θ(t)+θ′(t)∈R+.We now use the daily data from Jinan station withn=1. As Figure 4 shows, the graph of θ(t) for the whole year based on (3.35), the simulation of the mean monthly precipitation does not stay positive due to the value α which is not always positive. Moreover, for the model whose diffusion part has term Xt, once a negative Xt value is obtained, according to (3.30), the simulated value for the next step could never be calculated since p<1/2, which would lead to 2p-1<0 and therefore no value for Xn2p-1.Figure 4
Simulated speed of monthly rain in Jinan.Therefore, in order to ensure positivity,θ(t) must be modeled appropriately. Since θ(t) is a continuous sinusoidal function and the speed of mean reversion a is assumed constant, there should be a boundary for the value of α in (3.35) that ensures α>0. Thus we require,(3.36)θ(t)>-θ′(t)a,t=t0,…,tN.
If we only add a constant on the right hand of (3.25) to make sure the new θ(t) satisfies (3.35), we still have the same θ′(t). Therefore we must find maxt=t0,…,tN(-θ′(t)/a), denoted by M, so that the adjusted simulated mean speed of precipitation can be expressed as θad(t)=θ(t)+M+1, and θad′(t)=θ′(t).Continuing with our use of the Jinan data as an example, the results are shown in Figure5 and the corrector for θ(t) is −40.977.Figure 5
Adjusted simulated speed of monthly Rain in Jinan.To use the newθ in the algorithm, note that all the simulated points, Xt, have mean value θad(t). This change could also affect the diffusion part. In order to keep the process similar to the original one, we modify the algorithm as follows.(1)
E=max(Xn(1)-M-1,v), where v is a uniform random variable in (0, 0.1). If we use 0 to replace v, the diffusion part could be 0 if 2p-1>0 and infinity if 2p-1<0.(2)
Change (3.30) to
(3.37)Xn+1(0)=Xn+(θ′(n)+a(θ(n)-Xn(1)))Δtn+σEpΔWn+12σ2pE2p-1(ΔWn2-Δtn),Xn+1(1)=Xn+(θ′(n+1)+a(θ(n+1)-Xn+1(0)))Δtn+σEpΔWn+12σ2pE2p-1(ΔWn2-Δtn).(3)
After calculating all the simulated points, we need to change the whole list back down to the original level:(3.38)Xn=Xn(1)-M-1.
Note that if p<1/2, then positivity is not yet guaranteed for a large size simulation.


## 3.1. Historical Burn Analysis and Index Value Simulation

For Historical Burn Analysis, the basic assumption is that the historical data reasonably approximate future scenarios, allowing the price of the option to be calculated from data we already have. From (2.3), we can see that the price of a put option for drought is given by(3.1)V=1Ne-rτ(∑i=1Nmax⁡(K-Ii,0)),
where Ii is the value of the drought index for year i and N is the number of years of historical data.Our Index Value Simulation approach involves finding the best fit distribution for index values calculated using the Historical Burn Analysis approach and then sampling this distribution to produce possible future values of the index. The option price will then be calculated using (3.1).


## 3.2. Stochastic Simulation

Mean reversion stochastic processes are used to simulate the daily mean temperature and the speed of monthly rainfall over the entire year. The period’s rainfall and evapotranspiration are then obtained by summing the related simulated values. Finally, the option payoff is computed based on the simulation where the average discounted payoff is used as the option price.


### 3.2.1. Simulation of Mean Monthly Temperature

Although based on (2.6) and (2.7), only the mean monthly temperatures are needed to calculate the evapotranspiration, and those mean temperatures can be simulated directly.Mean Reversion Process
Average daily temperature is, of course, subject to seasonal changes. For one specific day in each year, however, we will assume that the temperature fluctuates around some mean value, which means we can choose a mean reversion stochastic process to simulate the behavior of daily temperature [2, 16].
We use the following Stochastic Differential Equation (SDE) to model daily temperature:(3.2)dTt=dTtm+a(Ttm-Tt)dt+σtdWt.
The solution is
(3.3)Tt=(Ts-Tsm)e-a(t-s)+Ttm+∫ste-a(t-τ)στdWτ,
where Tt is the temperature at time t, Ttm is the mean daily temperature at time t, a reflects the speed of the mean reversion, σt reflects the variation from the mean value at time t, and Wτ is a Wiener process.Parameter Estimation
Ttm in (3.2) represents simulated mean daily temperature. If the mean temperature is assumed to be a continuous variable then it is possible to simulate this curve by a sinusoidal function of the form:
(3.4)sin⁡(ωt+φ),
where t denotes the time, measured in days.
To account for local warming trends observed in the data,Ttm is modeled as having the form:(3.5)Ttm=A+Bt+Csin⁡(ωt+φ).
Expanding (3.5) gives us
(3.6)Ttm=a1+a2t+a3sin⁡(ωt)+a4cos⁡(ωt).
We then estimate the parameters in (3.6) by using the Gauss-Newton algorithm to solve(3.7)min⁡ξ⟦Tt-X2⟧,
where X is the data vector and ξ is the parameter vector (a1,a2,a3,a4).
Estimates ofξ then also allows us to calculate the parameters in (3.5) as follows:(3.8)A=a1,B=a2,C=a32+a42,φ=arctan(a4a3)-π.
The termσtdWt in (3.2) reflects the discrepancy of real temperature from the simulated mean temperature. We simplify the function σt as a piecewise constant function. Here we assume the value of σt is a constant number for each month.
As we see,σt reflects the variation of daily temperature, the first estimator is based on the quadratic variation of Tt:
(3.9)σμ2=1Nμ∑j=0Nμ-1(Tj+1-Tj)2,
where Nμ is the number of days in month μ.
If we discretize (3.2), for a given month μ, we have(3.10)Tj=Tjm-Tj-1m+aTj-1m+(1-a)Tj-1+σμϵj-1,
with j=1,…,Nμ, where {ϵj}j=1Nμ-1 follow the standard normal distribution. If we let T̅j be Tj-(Tjm-Tj-1m), then (3.10) becomes
(3.11)T̅j=aTj-1m+(1-a)Tj-1+σμϵj-1,
which we can see as a regression of today’s temperature on yesterday’s temperature. Therefore, the second estimator of σμ is
(3.12)σ̂μ2=1Nμ-2∑j=2Nμ(T̅j-âTj-1m-(1-â)Tj-1)2,
where â is an estimator of a.
To estimate the mean reversion parametera, we use the martingale estimation functions method. After collecting observations during n days, an efficient estimator ân of a can be obtained by solving the equation(3.13)Gn(ân)=0,
where
(3.14)Gn(a)=∑i=1nḃ(Ti-1;a)σi-12(Ti-E[Ti∣Ti-1]),
and where ḃ(Ti-1;a) denotes the derivative with respect to a of the term
(3.15)b(Tt;a)=dTtmdt+a(Ttm-Tt).
Therefore, we have
(3.16)ḃ(Ti-1;a)=Ttm-Tt.
To then solve (3.14), all we need is the formula Ti-E[Ti∣Ti-1]. As we know the solution of (3.2) is (3.3), and the expectation of the integration part is 0, when t=i and s=i-1, we have(3.17)E[Ti∣Ti-1]=(Ti-1-Ti-1m)e-a+Tim.
So we have
(3.18)Gn(a)=∑i=1nTi-1m-Ti-1σi-12[Ti-(Ti-1-Ti-1m)e-a-Tim].
It is easy to estimate the parameter a from
(3.19)ân=-log⁡(∑i=1n-1Yi-1(Ti-Tim)∑i=1nYi-1(Ti-1-Ti-1m)),
where
(3.20)Yi-1≡Ti-1m-Ti-1σi-12,i=1,2,…n.
Based on (3.15), we can see that ân is the unique estimation.Numerical Solution
Because it is difficult to use the integral solution in (3.3), we employ the implicit Milstein numerical method [17, 18]. The problem is then reduced to the following iterative procedure:(3.21)Tn+1(0)=Tn+(Tnm′+a(Tnm-Tn))Δt+σkΔW,Tn+1(1)=Tn+(Tn+1m′+a(Tn+1m-Tn+1(0)))Δt+σkΔW,
where Tn+1(1) is the simulated point we get, σk is the value of diffusion parameter depending on the nth point located in month k.


### 3.2.2. Simulation of the Speed of Monthly Rain

Transforming Monthly Rain to a Continuous Process
While precipitation does not behave continuously, the speed of precipitation can be assumed to change continuously in time. When rain starts, the speed of precipitation increases continuously from zero to its peak and decrease continuously to zero where it remains at zero until it rains again. In this formulation, the total amount of precipitation can be obtained by integrating the speed of precipitation over the period of interest, as given by the formula:(3.22)Pt1,t2=∫t1t2x(t)dt,
where x(t) represents the speed of rainfall at time t and Pt1,t2 represents the amount of precipitation during the period from time t1 to time t2.
The monthly rainfall, then, is essentially the monthly mean speed of precipitation with units (mm/m2·month). In order to simulate monthly rainfall, it will be sufficient to simulate the speed of rainfall using a mean reversion process.Mean Reversion Process
In the previous model of monthly temperature, diffusion was modeled as a piecewise functionσt, where σt was assumed constant for each month. In order to make the simulation of the speed of monthly rain more realistic, the diffusion can be modeled as a function of Xt according to the following SDE:(3.23)dXt=dθ(t)+a(θ(t)-Xt)dt+σXtpdWt,
with t≥0, a≥0, and θ(t)>0. As described in (3.2), θ(t) is the simulated mean speed of monthly rain at time t over a year and a is the mean reversion parameter. Here σ and p together reflect the volatility which is now dependent on previous precipitation.Parameter Estimation
The data we will utilize are for the Jinan climate station in eastern China, obtained from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/) for the period from January 1st, 1951 to December 31st, 2006.
Figure1 depicts the mean monthly rain speed over the 56-year period. As in the temperature model, the real mean speed of monthly precipitation will be assumed to fit a sinusoidal function of the form [19]:(3.24)θ(t)=m+αsin⁡(2π(t-υ)12),
where m is the mean of the sine curve, α determines the oscillation, and υ is the horizontal shift. To make the simulated curve closer to the historical experience, the simulation can be further improved by expanding θ in terms of a Fourier series [19]:
(3.25)θ(t)=m+∑i=0nαisin⁡((2i+1)π(t-υ)6).
The parameters inθ(t) can then be estimated using the Gauss-Newton method to solve the least squares problem. To do this, we choose n=1.
We now wish to find an unbiased estimator of the mean reversion parametera. The original estimation of a is [19]:
(3.26)â=1n∑i=0n-1Xi+1-Xi-θ(i+1)+θ(i)(θ(i)-Xi)Δ,
where Δ is the time increment. Because this estimation would change due to different lengths of data, we have the modification:
(3.27)ab̂=1|Ib|∑i∈IbXi+1-Xi-θ(i+1)+θ(i)(θ(i)-Xi)Δ,
with Ib={i=1,…,n:|θ(i)-Xi|>b}, b a positive real number. As b gets larger, ab̂ is convergent. Here we choose b=50 for a suitable estimation.
The diffusion parameters can be estimated by first squaring equation (3.23) to get(3.28)(dXt)2=a2(θ-Xt)2⋅dt⋅dt+a⋅(θ-Xt)⋅σ⋅Xtp⋅dt⋅dWt.
Then, based on Itô’s integration rule, we have
(3.29)ln⁡(dXt)2=2ln⁡(σ)+ln⁡(dt)+2pln⁡(Xt),
such that ln(dXt)2 and ln(dXt) have a linear relationship. The data can now be substituted into (3.29). We use linear regression method to estimate σ and p.Figure 1
Mean speed of monthly rain over 56 years.Numerical Solution
After the parameter estimation, based on the implicit Milstein method, the monthly precipitation can be simulated by the following:(3.30)Xn+1(0)=Xn+(θ′(n)+a(θ(n)-Xn(0)))Δtn+σXn(1)pΔWn+12σ2pXn(1)2p-1(ΔWn2-Δtn),Xn+1(1)=Xn+(θ′(n+1)+a(θ(n+1)-Xn+1(0)))Δtn+σXn(1)pΔWn+12σ2pXn(1)2p-1(ΔWn2-Δtn),
where Xn+1(1) is the simulated point at the (n+1)th step, n=0,…,N and N is the step size for one simulation process.


## 3.2.1. Simulation of Mean Monthly Temperature

Although based on (2.6) and (2.7), only the mean monthly temperatures are needed to calculate the evapotranspiration, and those mean temperatures can be simulated directly.Mean Reversion Process
Average daily temperature is, of course, subject to seasonal changes. For one specific day in each year, however, we will assume that the temperature fluctuates around some mean value, which means we can choose a mean reversion stochastic process to simulate the behavior of daily temperature [2, 16].
We use the following Stochastic Differential Equation (SDE) to model daily temperature:(3.2)dTt=dTtm+a(Ttm-Tt)dt+σtdWt.
The solution is
(3.3)Tt=(Ts-Tsm)e-a(t-s)+Ttm+∫ste-a(t-τ)στdWτ,
where Tt is the temperature at time t, Ttm is the mean daily temperature at time t, a reflects the speed of the mean reversion, σt reflects the variation from the mean value at time t, and Wτ is a Wiener process.Parameter Estimation
Ttm in (3.2) represents simulated mean daily temperature. If the mean temperature is assumed to be a continuous variable then it is possible to simulate this curve by a sinusoidal function of the form:
(3.4)sin⁡(ωt+φ),
where t denotes the time, measured in days.
To account for local warming trends observed in the data,Ttm is modeled as having the form:(3.5)Ttm=A+Bt+Csin⁡(ωt+φ).
Expanding (3.5) gives us
(3.6)Ttm=a1+a2t+a3sin⁡(ωt)+a4cos⁡(ωt).
We then estimate the parameters in (3.6) by using the Gauss-Newton algorithm to solve(3.7)min⁡ξ⟦Tt-X2⟧,
where X is the data vector and ξ is the parameter vector (a1,a2,a3,a4).
Estimates ofξ then also allows us to calculate the parameters in (3.5) as follows:(3.8)A=a1,B=a2,C=a32+a42,φ=arctan(a4a3)-π.
The termσtdWt in (3.2) reflects the discrepancy of real temperature from the simulated mean temperature. We simplify the function σt as a piecewise constant function. Here we assume the value of σt is a constant number for each month.
As we see,σt reflects the variation of daily temperature, the first estimator is based on the quadratic variation of Tt:
(3.9)σμ2=1Nμ∑j=0Nμ-1(Tj+1-Tj)2,
where Nμ is the number of days in month μ.
If we discretize (3.2), for a given month μ, we have(3.10)Tj=Tjm-Tj-1m+aTj-1m+(1-a)Tj-1+σμϵj-1,
with j=1,…,Nμ, where {ϵj}j=1Nμ-1 follow the standard normal distribution. If we let T̅j be Tj-(Tjm-Tj-1m), then (3.10) becomes
(3.11)T̅j=aTj-1m+(1-a)Tj-1+σμϵj-1,
which we can see as a regression of today’s temperature on yesterday’s temperature. Therefore, the second estimator of σμ is
(3.12)σ̂μ2=1Nμ-2∑j=2Nμ(T̅j-âTj-1m-(1-â)Tj-1)2,
where â is an estimator of a.
To estimate the mean reversion parametera, we use the martingale estimation functions method. After collecting observations during n days, an efficient estimator ân of a can be obtained by solving the equation(3.13)Gn(ân)=0,
where
(3.14)Gn(a)=∑i=1nḃ(Ti-1;a)σi-12(Ti-E[Ti∣Ti-1]),
and where ḃ(Ti-1;a) denotes the derivative with respect to a of the term
(3.15)b(Tt;a)=dTtmdt+a(Ttm-Tt).
Therefore, we have
(3.16)ḃ(Ti-1;a)=Ttm-Tt.
To then solve (3.14), all we need is the formula Ti-E[Ti∣Ti-1]. As we know the solution of (3.2) is (3.3), and the expectation of the integration part is 0, when t=i and s=i-1, we have(3.17)E[Ti∣Ti-1]=(Ti-1-Ti-1m)e-a+Tim.
So we have
(3.18)Gn(a)=∑i=1nTi-1m-Ti-1σi-12[Ti-(Ti-1-Ti-1m)e-a-Tim].
It is easy to estimate the parameter a from
(3.19)ân=-log⁡(∑i=1n-1Yi-1(Ti-Tim)∑i=1nYi-1(Ti-1-Ti-1m)),
where
(3.20)Yi-1≡Ti-1m-Ti-1σi-12,i=1,2,…n.
Based on (3.15), we can see that ân is the unique estimation.Numerical Solution
Because it is difficult to use the integral solution in (3.3), we employ the implicit Milstein numerical method [17, 18]. The problem is then reduced to the following iterative procedure:(3.21)Tn+1(0)=Tn+(Tnm′+a(Tnm-Tn))Δt+σkΔW,Tn+1(1)=Tn+(Tn+1m′+a(Tn+1m-Tn+1(0)))Δt+σkΔW,
where Tn+1(1) is the simulated point we get, σk is the value of diffusion parameter depending on the nth point located in month k.


## 3.2.2. Simulation of the Speed of Monthly Rain

Transforming Monthly Rain to a Continuous Process
While precipitation does not behave continuously, the speed of precipitation can be assumed to change continuously in time. When rain starts, the speed of precipitation increases continuously from zero to its peak and decrease continuously to zero where it remains at zero until it rains again. In this formulation, the total amount of precipitation can be obtained by integrating the speed of precipitation over the period of interest, as given by the formula:(3.22)Pt1,t2=∫t1t2x(t)dt,
where x(t) represents the speed of rainfall at time t and Pt1,t2 represents the amount of precipitation during the period from time t1 to time t2.
The monthly rainfall, then, is essentially the monthly mean speed of precipitation with units (mm/m2·month). In order to simulate monthly rainfall, it will be sufficient to simulate the speed of rainfall using a mean reversion process.Mean Reversion Process
In the previous model of monthly temperature, diffusion was modeled as a piecewise functionσt, where σt was assumed constant for each month. In order to make the simulation of the speed of monthly rain more realistic, the diffusion can be modeled as a function of Xt according to the following SDE:(3.23)dXt=dθ(t)+a(θ(t)-Xt)dt+σXtpdWt,
with t≥0, a≥0, and θ(t)>0. As described in (3.2), θ(t) is the simulated mean speed of monthly rain at time t over a year and a is the mean reversion parameter. Here σ and p together reflect the volatility which is now dependent on previous precipitation.Parameter Estimation
The data we will utilize are for the Jinan climate station in eastern China, obtained from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/) for the period from January 1st, 1951 to December 31st, 2006.
Figure1 depicts the mean monthly rain speed over the 56-year period. As in the temperature model, the real mean speed of monthly precipitation will be assumed to fit a sinusoidal function of the form [19]:(3.24)θ(t)=m+αsin⁡(2π(t-υ)12),
where m is the mean of the sine curve, α determines the oscillation, and υ is the horizontal shift. To make the simulated curve closer to the historical experience, the simulation can be further improved by expanding θ in terms of a Fourier series [19]:
(3.25)θ(t)=m+∑i=0nαisin⁡((2i+1)π(t-υ)6).
The parameters inθ(t) can then be estimated using the Gauss-Newton method to solve the least squares problem. To do this, we choose n=1.
We now wish to find an unbiased estimator of the mean reversion parametera. The original estimation of a is [19]:
(3.26)â=1n∑i=0n-1Xi+1-Xi-θ(i+1)+θ(i)(θ(i)-Xi)Δ,
where Δ is the time increment. Because this estimation would change due to different lengths of data, we have the modification:
(3.27)ab̂=1|Ib|∑i∈IbXi+1-Xi-θ(i+1)+θ(i)(θ(i)-Xi)Δ,
with Ib={i=1,…,n:|θ(i)-Xi|>b}, b a positive real number. As b gets larger, ab̂ is convergent. Here we choose b=50 for a suitable estimation.
The diffusion parameters can be estimated by first squaring equation (3.23) to get(3.28)(dXt)2=a2(θ-Xt)2⋅dt⋅dt+a⋅(θ-Xt)⋅σ⋅Xtp⋅dt⋅dWt.
Then, based on Itô’s integration rule, we have
(3.29)ln⁡(dXt)2=2ln⁡(σ)+ln⁡(dt)+2pln⁡(Xt),
such that ln(dXt)2 and ln(dXt) have a linear relationship. The data can now be substituted into (3.29). We use linear regression method to estimate σ and p.Figure 1
Mean speed of monthly rain over 56 years.Numerical Solution
After the parameter estimation, based on the implicit Milstein method, the monthly precipitation can be simulated by the following:(3.30)Xn+1(0)=Xn+(θ′(n)+a(θ(n)-Xn(0)))Δtn+σXn(1)pΔWn+12σ2pXn(1)2p-1(ΔWn2-Δtn),Xn+1(1)=Xn+(θ′(n+1)+a(θ(n+1)-Xn+1(0)))Δtn+σXn(1)pΔWn+12σ2pXn(1)2p-1(ΔWn2-Δtn),
where Xn+1(1) is the simulated point at the (n+1)th step, n=0,…,N and N is the step size for one simulation process.


## 3.3. Correction of the Model

### 3.3.1. Correction for the Mean Curve

For both temperature and precipitation models, the distance between the average of historical records and the simulated mean curve can also be simulated with another sinusoidal function like (3.25), with the exception that the period of this function may be different depending on the climate conditions of the location.Figures2 and 3 compare the simulated mean curve with and without the simulation of distance. It is clear that the inclusion of simulated distance provides results that are much closer to the historical record than the original one.Figure 2
Simulated daily temperature comparison.Figure 3
Simulated mean speed of rain comparison.


### 3.3.2. Maintaining the Positivity of the Rain Speed

SinceXt represents the speed function, it must be nonnegative to have a meaningful interpretation. Moreover, according to (3.23), negative Xt leads to no value for Xt2p-1. However, there is no guarantee that the simulation of (3.23) yields positive values, especially when Xt is close to 0 at the beginning of the year [17]. If we examine Figure 3, we see that the mean speed of rain is under the time axis at the beginning of the year.IfXt represents a stochastic process with(3.31)Prob({Xt>0∀t})=1,
the stochastic integration scheme possesses an eternal life time if(3.32)Prob({Xn+1>0∣Xn>0})=1.
Otherwise, we say it has a finite lifetime.If we can find a numerical method to solve this model which has an eternal life, as long as the initial value of historical precipitation is positive, we can make sure that all the following simulated points are positive. This is another reason to choose the implicit Milstein method to do the simulation [19].For the mean reversion process given by(3.33)dXt=(α-βXt)dt+σXtpdW,
with α,β,σ,p∈R+ and p>1/2, the Milstein method provides numerical positivity with the following restriction:(3.34)Δt<1σ2.Compared to (3.23), (3.33) will guarantee positivity with(3.35)α=a⋅θ(t)+θ′(t)∈R+.We now use the daily data from Jinan station withn=1. As Figure 4 shows, the graph of θ(t) for the whole year based on (3.35), the simulation of the mean monthly precipitation does not stay positive due to the value α which is not always positive. Moreover, for the model whose diffusion part has term Xt, once a negative Xt value is obtained, according to (3.30), the simulated value for the next step could never be calculated since p<1/2, which would lead to 2p-1<0 and therefore no value for Xn2p-1.Figure 4
Simulated speed of monthly rain in Jinan.Therefore, in order to ensure positivity,θ(t) must be modeled appropriately. Since θ(t) is a continuous sinusoidal function and the speed of mean reversion a is assumed constant, there should be a boundary for the value of α in (3.35) that ensures α>0. Thus we require,(3.36)θ(t)>-θ′(t)a,t=t0,…,tN.
If we only add a constant on the right hand of (3.25) to make sure the new θ(t) satisfies (3.35), we still have the same θ′(t). Therefore we must find maxt=t0,…,tN(-θ′(t)/a), denoted by M, so that the adjusted simulated mean speed of precipitation can be expressed as θad(t)=θ(t)+M+1, and θad′(t)=θ′(t).Continuing with our use of the Jinan data as an example, the results are shown in Figure5 and the corrector for θ(t) is −40.977.Figure 5
Adjusted simulated speed of monthly Rain in Jinan.To use the newθ in the algorithm, note that all the simulated points, Xt, have mean value θad(t). This change could also affect the diffusion part. In order to keep the process similar to the original one, we modify the algorithm as follows.(1)
E=max(Xn(1)-M-1,v), where v is a uniform random variable in (0, 0.1). If we use 0 to replace v, the diffusion part could be 0 if 2p-1>0 and infinity if 2p-1<0.(2)
Change (3.30) to
(3.37)Xn+1(0)=Xn+(θ′(n)+a(θ(n)-Xn(1)))Δtn+σEpΔWn+12σ2pE2p-1(ΔWn2-Δtn),Xn+1(1)=Xn+(θ′(n+1)+a(θ(n+1)-Xn+1(0)))Δtn+σEpΔWn+12σ2pE2p-1(ΔWn2-Δtn).(3)
After calculating all the simulated points, we need to change the whole list back down to the original level:(3.38)Xn=Xn(1)-M-1.
Note that if p<1/2, then positivity is not yet guaranteed for a large size simulation.


## 3.3.1. Correction for the Mean Curve

For both temperature and precipitation models, the distance between the average of historical records and the simulated mean curve can also be simulated with another sinusoidal function like (3.25), with the exception that the period of this function may be different depending on the climate conditions of the location.Figures2 and 3 compare the simulated mean curve with and without the simulation of distance. It is clear that the inclusion of simulated distance provides results that are much closer to the historical record than the original one.Figure 2
Simulated daily temperature comparison.Figure 3
Simulated mean speed of rain comparison.


## 3.3.2. Maintaining the Positivity of the Rain Speed

SinceXt represents the speed function, it must be nonnegative to have a meaningful interpretation. Moreover, according to (3.23), negative Xt leads to no value for Xt2p-1. However, there is no guarantee that the simulation of (3.23) yields positive values, especially when Xt is close to 0 at the beginning of the year [17]. If we examine Figure 3, we see that the mean speed of rain is under the time axis at the beginning of the year.IfXt represents a stochastic process with(3.31)Prob({Xt>0∀t})=1,
the stochastic integration scheme possesses an eternal life time if(3.32)Prob({Xn+1>0∣Xn>0})=1.
Otherwise, we say it has a finite lifetime.If we can find a numerical method to solve this model which has an eternal life, as long as the initial value of historical precipitation is positive, we can make sure that all the following simulated points are positive. This is another reason to choose the implicit Milstein method to do the simulation [19].For the mean reversion process given by(3.33)dXt=(α-βXt)dt+σXtpdW,
with α,β,σ,p∈R+ and p>1/2, the Milstein method provides numerical positivity with the following restriction:(3.34)Δt<1σ2.Compared to (3.23), (3.33) will guarantee positivity with(3.35)α=a⋅θ(t)+θ′(t)∈R+.We now use the daily data from Jinan station withn=1. As Figure 4 shows, the graph of θ(t) for the whole year based on (3.35), the simulation of the mean monthly precipitation does not stay positive due to the value α which is not always positive. Moreover, for the model whose diffusion part has term Xt, once a negative Xt value is obtained, according to (3.30), the simulated value for the next step could never be calculated since p<1/2, which would lead to 2p-1<0 and therefore no value for Xn2p-1.Figure 4
Simulated speed of monthly rain in Jinan.Therefore, in order to ensure positivity,θ(t) must be modeled appropriately. Since θ(t) is a continuous sinusoidal function and the speed of mean reversion a is assumed constant, there should be a boundary for the value of α in (3.35) that ensures α>0. Thus we require,(3.36)θ(t)>-θ′(t)a,t=t0,…,tN.
If we only add a constant on the right hand of (3.25) to make sure the new θ(t) satisfies (3.35), we still have the same θ′(t). Therefore we must find maxt=t0,…,tN(-θ′(t)/a), denoted by M, so that the adjusted simulated mean speed of precipitation can be expressed as θad(t)=θ(t)+M+1, and θad′(t)=θ′(t).Continuing with our use of the Jinan data as an example, the results are shown in Figure5 and the corrector for θ(t) is −40.977.Figure 5
Adjusted simulated speed of monthly Rain in Jinan.To use the newθ in the algorithm, note that all the simulated points, Xt, have mean value θad(t). This change could also affect the diffusion part. In order to keep the process similar to the original one, we modify the algorithm as follows.(1)
E=max(Xn(1)-M-1,v), where v is a uniform random variable in (0, 0.1). If we use 0 to replace v, the diffusion part could be 0 if 2p-1>0 and infinity if 2p-1<0.(2)
Change (3.30) to
(3.37)Xn+1(0)=Xn+(θ′(n)+a(θ(n)-Xn(1)))Δtn+σEpΔWn+12σ2pE2p-1(ΔWn2-Δtn),Xn+1(1)=Xn+(θ′(n+1)+a(θ(n+1)-Xn+1(0)))Δtn+σEpΔWn+12σ2pE2p-1(ΔWn2-Δtn).(3)
After calculating all the simulated points, we need to change the whole list back down to the original level:(3.38)Xn=Xn(1)-M-1.
Note that if p<1/2, then positivity is not yet guaranteed for a large size simulation.


## 4. Analysis of Results

Using temperature and precipitation data from the Jinan climate station, and the curve of consumptive crop coefficient based on McGuinness and Bordne [14], we price a drought put option contract based on three different within-year data periods: (1) January to December (Table 1); April to August (Table 2); April to June (Table 3). In each case, we set the strike level at K=0.7 and use each of our three methods. Each of Tables 1 through 3 shows the option prices based on data from four different historical spans: 1951–2006, 1967–2006, 1977–2006, and 1987–2006.Table 1
Option price comparison for months 1–12.
1951–20061967–20061977–20061987–2006Historical Burn Analysis0.17640.18130.17990.1561Index Value—Weib. Distribution0.17690.18160.18050.1558Stochastic Simulation0.16660.17820.17390.1502Table 2
Option price comparison for months 4–8.
1951–20061967–20061977–20061987–2006Historical Burn Analysis0.18720.18770.17970.1453Index Value—Ext. Val. Distribution0.18970.19310.18430.1535Stochastic Simulation0.18460.19880.18730.1584Table 3
Option price comparison for months 4–6.
1951–20061967–20061977–20061987–2006Historical Burn Analysis0.19630.19270.16860.1505Index Value—Ext. Val. Distribution0.19280.18970.16470.1522Stochastic Simulation0.16160.16760.15200.1447Comparing the option prices for the full 56-year span, we see that the results from the Historical Burn Analysis and the Index Value approaches are very similar to one another, while the Stochastic Simulation method yields different results. This difference is attributable to the inaccuracy of the simulation of the speed of monthly rain. The Stochastic Simulation method does, however, yield results that have the lowest variance across the four columns.Looking at the 1987–2006 results in the last column in each of Tables1–3, we see that, for each method, the results are quite different from those obtained using the full 56 years of data. This points to the overall sensitivity of pricing to the historical span of time upon which the pricing is based.


## 5. Future Work

To simplify our analysis, we have ignored a number of financial factors. Future work could incorporate the relationship between a drought index and farmers’ profits, the possibility of changing interest rates, and the potential market price of risk [20].In addition, a number of improvements could be made to the weather models we have used to price drought contracts. In particular, the price of the drought option clearly depends on the joint distribution of the temperatures and the precipitation at maturity. To simplify our analysis, we have simulated these two variables independently. In future work, a joint model could be developed. We have also not considered the location limitations of the climate models. In particular, the climate stations which collect precipitation data are typically located in large urban areas, far from the rural areas where farming is concentrated. Because the levels of precipitation might be quite different between even neighboring urban and rural areas, future work should incorporate this spatial basis risk.



---

*Source: 251835-2012-04-11.xml*