Modeling Drought Option Contracts

In this paper we introduce a new financial weather derivative – a drought option contract – to protect agricultural producers from 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 estimate the value of a hypothetical drought contract based on data from the Jinan weather station located in a dry region of China.


Introduction
Since the dawn of agriculture, farmers have been vulnerable to the effects of drought.Indeed, in regions that rely primarily on rain-fed agriculture, drought continues to be one of the greatest risks farmers face.Because drought is a highly positively spatially-correlated event, serious disruptions to regional food supplies may result.And in addition to the immediate crop losses themselves, the resulting economic impact of those losses to a region's agricultural sector may undermine its stability going forward.
Historically, the only financial instruments available to mitigate the risks associated with drought were crop insurance contracts.These contracts have the obvious advantage of minimizing the financial losses incurred by farmers, thus stabilizing their operations.But because the marginal benefits of crop damage-mitigating strategies are offset by reductions in the expected value of the insurance claim, the private benefits farmers 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 undermine, thus mitigating real crop losses (Gronberg, 2007).Moreover, relative to insurance markets, 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 (Alaton, 2002;Mraoua, 2005;Yoo, 2003;Cao, 2004), 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 (2000), Brody et al (2002), Richards et al (2004), Cao and Wei (2004), Geman et al (2005), and Taylor and Buizza (2006) -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 the value of the crop yield they affect, however, a basis risk remains.
In this paper, we propose and model 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 in order to reduce basis risk.
The balance of the paper is organized as follows.In Section 2, we describe the structure of the drought option contract, the drought index, and all the parameters which can affect the value of this drought index.In Section 3, we introduce three different models to estimate the option prices, while section 4 uses data from Jinan station to simulate the option prices based on those three models.In Section 5, we compare the three models and discuss possible improvements and future work.

Drought Options
We model our derivative contract as a European put option on an agricultural drought index, I, with strike level K, contract period τ , and the tip, S, which reflects the relationship between the drought index and the loss of farmers' income.At maturity, the payoff function is: If I 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.
As noted above, because no-arbitrage pricing methods are unsuitable for pricing weather options, we will use the equilibrium approach, where the option is priced as the present value of the expectation of its payoff (Cao, 2004;Mußhoff, 2006).
The equilibrium approach finds the value V of an option according to the following expectation: where τ is the contract expiry, f (I) is payoff of the contract, r(t) is the risk-free interest rate, I is the value of drought index, and Q is the probability measure for the expectation which includes all possible values of I.
To simplify the problem, we will assume that S = 1 and that r(t) is constant over the length of the contract, τ .Equation ( 2) can then be rewritten as Our problem is thus reduced to choosing an appropriate drought index, setting the strike level, and estimating the expectation.

Drought Index
Because our interest lies in developing a tool to help farmers hedge against the risks of reduced agricultural yield, our analysis will be based on the Reconnaissance Drought Index (RDI) (Tsakiris, 2005;Tsakiris, 2007), which measures the severity of drought in agricultural terms. 1 2 1 The American Meteorological Society divides drought into four categories: meteorological/climatological, agricultural, hydrological, and socioeconomic (Richard, 2002).Scientists have designed many different drought indices, each satisfying a different need (Richard, 2002;Keyantash, 2002).The most widely-referenced indices are the Palmer's Drought Index (PDI) and the Standardized Precipitation Index (SPI), both of which focus on meteorological or hydrological drought (Alley, 1984;Guttman, 1997;Guttman, 1999).Agricultural drought is concerned with climate factors including evapotranspiration, temperature, precipitation, and planting dates (Boken, 2005;Hanson, 1991).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 (P ET ) is the evapotranspiration that could occur were there sufficient water supply.
2 There are a number of advantages to the RDI (Tsakiris, 2006): (1) although equation (4) calculates the value for a year, the RDI can be calculated for any period of time; (2) because it considers both water supply (precipitation) and the water demand of the crop (evapotranspiration), it an well-suited for specifically measuring agricultural drought; We will use the initial value of the RDI for the option underlying, which can be calculated for any period in one year as: (4) where P ij and P ET ij are the precipitation and potential evapotranspiration, respectively, for the jth month of the ith year, m is the number of observed months in a given year, and N is the number of observed years.

Potential Evapotranspiration and Actual Evapotranspiration
Potential Evapotranspiration is difficult to measure and use in practice because it depends on the type of plants, the type of soil, and the the climatic conditions.For practical purposes, we will therefore use the actual, rather than the potential evapotranspiration in our estimation of the RDI.
This gives us a new Adjusted − RDI that can be written as: , i = 1 . . .N and j = 1 . . .m (5) where ET ij is the actual evapotranspiration in the jth month of the ith year.
To estimate the ET , we will use the temperaturebased methods of Blaney-Criddle, considered to be more accurate than the approach of Thornthwaite (McGuinness, 1972; USDA, 1964). 3  The formula for the SCS Blaney-Criddle method is: where k t = 0.0311T + 0.24 ( 7) and (3) because it incorporates more climate factors than PDI or SPI (see note 1), including day hour and crop coefficients, it may lead to a more accurate measure.
3 Scientists have used a number of methods to estimate evapotranspiration (Jensen, 1990).These can be divided into three categories (Xu, 2002): (1) mass-transfer based methods; (2) radiation-based methods; and (3) temperature-based methods.Our use of a temperaturebased method is necessitated by data limitations.and and where the monthly ET measured in inches, k c is the consumptive crop coefficient for the SCS version, T is monthly mean temperature, and d is the percentage of daylight hours.Note that when air temperature is no more than 1.67 • C, the value of k t is a constant value of 0.3.

Drought Option Pricing Models
We now wish to calculate possible values of the Adjusted − RDI in order to estimate the option price.Common techniques to calculate index values include Historical Burn Analysis, Index Value Simulation, and Stochastic Simulation (Mußhoff, 2006).

Historical Burn Analysis and Index
Value Simulation Model

Historical Burn Analysis
The basic assumption of this method is that the historical data reasonably approximates future scenarios, allowing the price of the option to be calculated from data we already have.
From Equation (3), we can see that the price of a put option for drought is given by where I i is the value of the drought index for the ith year during that specific period, and N is the number of years of data.
We choose a put option contract because of the characteristics of the RDI; based on the relationship of precipitation and evapotranspiration, the smaller the RDI value is, the more serious the drought we have.

Index Value Simulation
In this analysis, the index is considered to be the only variable for determining the drought option.
With the historical burn analysis approach, only observed scenarios captured in the historical data are used to calculate the option value.As a result, extreme events may not be adequately represented.By contrast, index value simulation considers not only the current and past values, but also the possible future values of the index.This is achieved by finding the best fit distribution for index values calculated from burn analysis and then sampling this distribution to produce new possible index values.The option price is then calculated again using Equation ( 9).4

Stochastic Simulation
For the purposes of our analysis, 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.

Simulation of Mean Monthly Temperature
Although based on the equations ( 6), (7), and (8), only the mean monthly temperatures are needed to calculate the evapotranspiration, and it is relatively simple to simulate those mean temperatures 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 (Alaton, 2002;Benth, 2007).
We will use the following Stochastic Differential Equation (SDE) to model daily temperature: The solution is where T t is the temperature at time t, T m t 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
T m t in equation ( 10) 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: sin(ωt + ϕ), where t denotes the time, measured in days.In order to account for trends observed in the data, due to local warming, such as urbanization, or global warming, the mean temperature is increased slightly each year through a linear relationship with respect to time t.Therefore, the mean temperature is modeled as T m t , which has the form Expanding equation ( 13), we have We can use the Gauss-Newton method to make the data fit this function, by solving min ξ where X is the data vector, ξ is the parameter vector (a 1 , a 2 , a 3 , a 4 ).
After estimating ξ, we also have the parameters in Equation ( 13) by changing their form: The term σ t dW t in equation ( 13) reflects the discrepancy of real temperature from the simulated mean temperature.To make the problem easier, 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 T t : where N µ is the number of days in month µ.
If we discretize equation ( 10), for a given month µ, we have which we can see as a regression of today's temperature on yesterday's temperature.Therefore, the second estimator of σ µ is where â is an estimator of a, which is described as follows.
To estimate the mean reversion parameter a, we can use the martingale estimation functions method.After we collect observations during n days, an efficient estimator ân of a can be obtained by solving the equation where and where ḃ(T i−1 ; a) denotes the derivative with respect to a of the term Therefore, we have To then solve equation ( 25), the only thing we need to know is the formula of As we know the solution of Equation ( 10) is Equation ( 11), and the expectation of the integration part is 0, when t = i and s = i − 1, we have So we have (29) It is easy to estimate the parameter a from ân = − log( where Based on Equation (26), we can see that ân is the unique estimation.

Numerical Solution
Because it is difficult to use the integral solution in equation ( 11), we employ the implicit Milstein numerical method (Kahl, 2004;Han, 2005).
The problem is then reduced to the following iterative procedure: T (1) (1) n+1 is the simulated point we get, σ k is the value of diffusion parameter depending on the nth point located in month k.

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: where x(t) represents the speed of rainfall at time t and P t 1 ,t 2 represents the amount of precipitation during the period from time t 1 to time t 2 .
The monthly rainfall, then, is essentially the monthly mean speed of precipitation with units mm/(m 2 • 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 X t according to the following SDE: with t ≥ 0, a ≥ 0, and θ(t) > 0. As described in equation ( 10), θ(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
Figure 1 depicts the mean monthly rain speed over 56 years using the daily data for the Jinan weather station.As in the temperature model, the real mean speed of monthly precipitation will be assumed to fit a sinusoidal function of the form (van Emmerich, 2005): 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 (van Emmerich, 2005):

Historical Mean Speed for each day
For the parameters in θ(t), one can use the Gauss-Newton method to solve the least squares problem.However, according to equation ( 37), we must first choose how many α i to use for the specific problem.
We now wish to find an unbiased estimator of the mean reversion parameter a.The original estimation of a is (van Emmerich, 2005): where ∆ is the time increment.Because this estimation would change due to different lengths of data, we have the modification: with Here we choose b = 50 for a suitable estimation.
The diffusion parameters can be estimated by first squaring equation ( 35) to get (40) Then, based on Itô's integration rule, we have ln(dX t ) 2 = 2 ln(σ) + ln(dt) + 2p ln(X t ), ( 41) such that ln(dX t ) 2 and ln(dX t ) have a linear relationship.The data list can now be substituted into equation ( 41), from which one can use the linear regression method to get σ and p directly.

Numerical Solution
After the parameter estimation, based on the implicit Milstein method, the monthly precipitation can be simulated by the following formulas: X (1) where X (1) n+1 is the simulated point at the (n + 1)−th step, n = 0...N and N is the step size for one simulation process.

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 equation ( 37), with the exception that the period of this function may be different depending on the climate conditions of the location.Since X t represents the speed function, it must be non-negative to have a meaningful interpretation.Moreover, according to equation ( 35), negative X t leads to no value for X 2p−1 t .However, there is no guarantee that the simulation of equation ( 35) yields positive values, especially when X t is close to 0 at the beginning of the year (Kahl, 2004).If we examine Figure 3, we see that the mean speed of rain is under the time axis at the beginning of the year.
If X t represents a stochastic process with Prob({X t > 0 for all t}) = 1, ( the stochastic integration scheme possesses an eternal life time if 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 (van Emmerich, 2005).
For the mean reversion process given by with α, β, σ, p ∈ R + and p > 1 2 , the Milstein method provides numerical positivity with the following restriction: Origin Simulated Mean Monthly Rain We now use the daily data from Jinan station with n = 1.As Figure 4 shows, the graph of θ(t) for the whole year based on equation ( 48), 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 X t , once a negative X t value is obtained, according to equations ( 42) and ( 43), 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 X 2p−1 n .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  48) that ensures α > 0. Thus we require, If we only add a constant on the right hand of Equation ( 37) to make sure the new θ(t) satisfies Equation ( 48), we still have the same θ (t).Therefore we must find max t=t 0 ,...,t N − θ (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 Figure 5 and the corrector for θ(t) is -40.977.
To use the new θ in the algorithm, we need to notice that all the simulated points X t here 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 have to modify the algorithm as follows: 1 , 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. 42) and (43) to

change equations (
and 3. After calculating all the simulated points, we need to change the whole list back down to the original level: Note that if p < 1 2 , then positivity is not yet guaranteed for a large size simulation.

Data
Based on Section 2, we need the data for temperature, precipitation, day hour length, and crop coefficient curve.In this article, the climate data are from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/).They have data for more than 700 climate stations since the year 1951.
In order to minimize the error, we use those station's observing data where the percentage of missing points is no more than 11.Missing values are imputed by substituting random numbers from the normal distribution.We choose data from Jinan station first, which has a history of drought.It is located 37

Option Price with Different Data Lengths
In this part, based on data from the Jinan station, we set a drought put option contract for three different periods with different strike level: from K=0.5.Because for some locations we cannot obtain enough observations, Tables 1-3 show the option prices based on different data lengths for each period for all three models.
After calculating the RDI list, we have the values of aridity index for the three periods are 0.51751, 0.52899, and 0.31727.This value would, of course, change if we instead chose to use only some subset of the data.Therefore we only calculate the option price with a constant strike level.
The results for the three models are generally similar for each column.When we estimate the option prices with data from 56 years, the results from the stochastic simulation model are not as close as the other two.It is reasonable because the inaccuracy of the simulation of the speed of monthly rain.
In the last column in these three tables, the option prices change significantly regardless of the model used, which means the time limitation would always exist in our models, even though at the beginning of Section 3, we hope the index value simulation model and stochastic simulation model can predict more possible values of RDI for Equation (9).Comparing the results from data over 20 years and those over 56 years for each model, we could see that the difference from historical burn analysis is always bigger than those from the other two.Comparing the prices in each row in Table 3, we see that the prices based on a stochastic process is more stable than the other two even though the original option prices based on data over 56 years is not accurate enough.
Finally, if sufficient data are not used, the estimated option prices will generally be inaccurate.However, even though the lack of data would always affect the accuracy of the estimation of option prices, index value simulation model and stochastic simulation model can generate reasonable results compared to historical burn analysis model.Even though our stochastic simulation model cannot describe the behavior of precipitation well, it is still the most stable result.

Conclusion and Future Work
In this article, we have used three models to estimate the value of the drought contract for the period from January to December, April to August and April to June in Jinan with different strike levels.The three methods are very similar to one another if we have enough historical data, suggesting that each method is reasonable.
If we wish to put this contract into practise, however, there may be other things that need to be considered.First, we should consider the relationship between the drought index and the influence of drought on farmers' incomes.Second, in developing countries with poor irrigation systems, it is difficult to obtain sufficient historical data sets, making it impossible to get an accurate description of drought in those locations and, therefore, impossible to use the historical analysis model.The index value simulation and daily index process can provide more information about future possible values of the drought index.The results from the stochastic simulation model are not as good as index value simulation because of the inaccurate precipitation model.
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, given that the market is not complete (Cao, 2004).
A number of improvements could also be made to the weather models we have used to price drought contracts.For example, we have used the mean reversion process to simulate monthly rain.Although this method yields reasonable results when we simulate behaviour in the long term, we cannot say that the model is ideal given that the standard deviation from our simulations fluctuates considerably from the historical figures.In addition, we have not considered the locational limitations of the climate models.Climate stations that collect precipitation data are typically located in large cities, far from the farmland where the drought contract would be used and where the precipitation may be significantly different.Future work could incorporate spatial basis risk into our model.
Finally, 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 will be developed.K=0.7 1951-20061977-20061987-2006 Historical Burn with j = 1, ..., N µ , where { j } Nµ−1 j=1 follow the standard normal distribution.If we let Tj be T j − (T m j − T m j−1 ), then equation (21) becomes Tj

Figure 1 :
Figure 1: Mean Speed of Monthly Rain over 56 years

Figures 2
Figures 2 and 3 compare the simulated mean curve with and without the simulation of distance.It is clear that the that the inclusion of simulated distance provides results that are much closer to the historical record than the original one.

Figure
Figure 2: Simulated Daily Temperature Comparison

Figure 3 :
Figure 3: Simulated Mean Speed of Rain Comparison

Figure 4 :
Figure 4: Simulated Speed of Monthly Rain in Jinan Compared to equation (35), equation (46) will guarantee positivity with

Figure 5 :
Figure 5: Adjusted Simulated Speed of Monthly Rain in Jinan • N in the east part of China.Data included observations from January 1st, 1951 to December 31st, 2006.

Table 3 :
Option Price Comparison for Months 4-6