Plant growth stages and weather index insurance design

Abstract

Given the assumption that weather risks affect crop yields, we designed a weather index insurance product for soybean producers in the US state of Illinois. By separating the entire vegetation cycle into four growth stages, we investigate whether the phase-division procedure contributes to weather–yield loss relation estimation and, hence, to basis risk mitigation. Concretely, supposing stage-variant interaction patterns between temperature-based weather index growing degree days and rainfall-based weather index cumulative rainfall, a nonparametric weather–yield loss relation is estimated by a generalized additive model. The model includes penalized B-spline (P-spline) approach based on nonlinear optimal indemnity solutions under the expected utility framework. The P-spline analysis of variance (PS-ANOVA) method is used for efficient estimation through mixed model re-parameterization. The results indicate that the phase-division models significantly outperform the benchmark whole-cycle ones either under quadratic utility or exponential utility, given different levels of risk aversions. Finally, regarding hedging effectiveness, the expected utility ratio between the phase-division contract and the whole-cycle contract, and the percentage changes of mean root square loss and variance of revenues support the proposed phase-division contract.

1. Introduction

Since its introduction in the 1990s, weather index insurance and weather derivatives have been regarded as risk management instruments offering numerous hedging possibilities to farmers and the agricultural sector. Potential advantages attributed to weather index insurance include the remission of loss assessment, transparency of contracts, and mitigation of moral hazard and adverse selection (e.g., Barnett et al., 2008; Kellner & Mußhoff, 2011). However, thus far, weather index insurance products have not met the expectations that they have raised as financial risk management tools in the agricultural sector. As noted by Lin et al. (2015), the uptake of these products is relatively low. Several factors have been identified in the literature that may explain the reluctance of farmers to adopt these products, including nontransparent pricing mechanism (e.g., Xu et al., 2008), high costs due to systemic weather risk (Okhrin et al., 2013), and, most importantly, poor hedging effectiveness (e.g., Mußhoff et al., 2011; Pelka & Mußhoff, 2013).

Hedging effectiveness is closely related to basis risk, that is, the difference between actual yield losses and indemnification (Woodard & Garcia, 2008). Following Dalhaus et al. (2018), basis risk can be separated into three subcomponents: geographical basis risk, design basis risk, and temporal basis risk. Much research has been conducted to improve the hedging effectiveness of weather index insurance and, thus, the willingness to pay for these products by addressing the different components of basis risk. For example, geographic basis risk, which arises from various weather conditions at the farm location and the reference station of the index product, can be mitigated by regional diversification of insurance products (Ritter et al., 2014) or by spatial interpolation techniques (Cao et al., 2015). Design basis risk can be controlled by several parameters and methodological choices when modeling the complex weather–yield relationship. First, appropriate weather variables should be selected. Undoubtedly, rainfall and temperature are essential drivers of plant growth. Still, there is no consensus on the number and exact specification of the weather variables to be included in a yield model. Often used weather variables are growing degree days (GDD), cumulative rainfall (CR), dry spells, and the standardized precipitation index (cf. Turvey, 2001; Stoppa & Hess, 2003; Hill et al., 2019; Okpara et al., 2017). Moreover, although not fully adopted, it is widely accepted that the weather–yield relation is nonlinear (e.g., Schlenker & Roberts, 2009). However, the functional form of this relation is plant-specific. It must be adapted to the regional production environment, which calls for a flexible modeling approach, either parametrically or nonparametrically (e.g., Bokusheva, 2011; Delerce et al., 2016). While early work by Vedenov & Barnett (2004) employed traditional regression methods to model the impact of weather variables on agricultural yields, machine learning techniques, specifically artificial neural networks, are increasingly used to accomplish this task (cf. Schmidt et al., 2022). Most recently, B-spline and P-spline methods have been proposed as a nonparametric alternative to model the weather–yield relation by Tan & Zhang (2020) and Bucheli et al. (2022). An additional insertion point for reducing design basis risk is the specification of the indemnity function. The majority of existing empirical applications, such as Vedenov & Barnett (2004) and Okhrin et al. (2013), assumed a stepwise linear indemnity function, in which the payoff is triggered by a threshold of the weather index, increases linearly with the index, and is often capped at a maximum level. However, Zhang et al. (2019) extended Raviv’s (1979) seminal work and proved that optimally designed weather index insurance consists of a nonlinear indemnity function. Finally, temporal basis risk can be considered as a particular aspect of design basis risk that refers to the inexpedient choice of the insured period (Conradt et al., 2015).

Agronomy and crop science assert that the vulnerability of crops concerning drought and temperature stress varies during the plant growth cycle. It is well known that the demand for water increases with the plant and leaf area until the plant begins to mature (e.g., Jensen, 1968). Nielsen & Nelson (1998) studied the effects of water deficit at various growth stages of black beans. They concluded that the yield reduction is highly sensitive to limited precipitation at the flowering stage. Moreover, exposure to temperature extremes at the reproductive stage greatly affects the production of plants (Hatfield & Prueger, 2015). So far, agronomic knowledge has been mainly used to determine the beginning and end dates of the accumulation period of weather indices, for example, GDD or CR (Dalhaus et al., 2018). Schierhorn et al. (2021) implemented machine learning techniques to explore the contribution of weather and climate to winter wheat yields during different plant development stages. However, only a few articles, such as Shi & Jiang (2016), explored the effect of weather variables and their interaction in different plant growth stages in the context of weather index insurance.

Against this backdrop, this paper aims to investigate whether the basis risk of weather index insurance can be further reduced by a growth-phase-dependent estimation of the weather–yield loss relationship using a flexible statistical modeling approach. To this end, we employ a generalized additive model (GAM) and several B-spline bivariate row tensor product smoothers (de Boor, 1978; Hastie & Tibshirani, 1987), with each smooth function representing a plant growth phase working together to fit the response variable, which is yield loss. So far, applications of B-spline methods in agricultural insurance are rare. Price et al. (2019) used the B-spline method to model the relationship between land quality, insurance coverage rate, and premium rate. Tan & Zhang (2020) and Bucheli et al. (2022) propagated this method to estimate the weather–yield relation. Although B-splines are attractive for nonparametric modeling, their specification comes with various challenges. Specifically, complicated knot selection schemes, including numbers and positions, have been the subject of much research (e.g., Kooperberg & Stone, 1991; Ruppert, 2002). In addition, to control overfitting, the penalties are pivotal, and many variations have been developed to address this issue, such as the thin plate penalty function (Wood, 2003; Eilers et al., 2015). Among the penalty variations, the difference penalty on adjacent coefficients is a particular branch of research called penalized B-spline (P-spline) (cf. Eilers & Marx, 1996). Using the equivalence between P-spline GAM and the mixed model, we apply a P-spline analysis of variance (PS-ANOVA) method. By PS-ANOVA, the knot selection can be avoided, and the smoothing parameter can be estimated in a computationally efficient manner by the decomposition of B-spline row tensor product smoothers, as noted by Lee & Durbán (2011) and Lee et al. (2013). The contribution of our paper is twofold. First, we explore whether breaking down the whole vegetation period into several growth stages can significantly improve the performance of weather index insurance. Second, we investigate the application of nonparametric methods, namely GAM and PS-ANOVA, in a nonlinear weather–yield loss relation estimation.

The remainder of this paper is organized as follows: section 2 introduces the optimal indemnity framework and the statistical approaches, GAM and PS-ANOVA. Section 3 describes the data, study area, and model estimation. Moreover, this section provides results and discusses the performance of a phase-division contract compared to a benchmark contract that includes the entire vegetation period of the crop. Section 4 draws conclusions on the optimal design of weather index insurance and evaluates the proposed statistical procedure.

2. Methodology

2.1. Optimal indemnity framework

In this paper, we take up the theoretical framework developed by Zhang et al. (2019), who extended Raviv’s (1979) seminal work on optimal insurance design. In that framework, an optimal indemnity function is derived from an expected utility maximization problem, in which a risk-averse insurance buyer (a farmer) is exposed to (yield) loss $Y$ and acquires insurance from an insurer at price $P$ that grants an indemnity payment $I$ . Maximizing the insured’s expected utility is carried out under a participation constraint of the insurer. Raviv (1979) proved that the optimal indemnity function has the structure of a coinsurance above a deductible. Moreover, in the case of a risk-neutral insurer, the optimal indemnity function is linear. This setting, however, is not directly applicable to the design of index-based insurance, because it does not account for basis risk. To capture this aspect, Zhang et al. (2019) modeled indemnity payments as a function of the weather index $X$ , that is, $I=I(X)$ . Normalizing the crop market price to one, income loss is equivalent to yield loss $Y$ . Yield loss and the weather index were assumed by Zhang et al. (2019) to have a joint probability density function $f(y, x)$ . Further, taking a risk-neutral insurer, the optimization problem can be stated as follows:

(1) \begin{align} \begin{cases} \max_{\mathit{I \epsilon }\mathcal{I}} J\!\left(I\right)=E\!\left[U\!\left\{w+I\!\left(X\right)-Y-\left(1-\tau \right)P\right\}\right]\\[5pt] \quad\qquad\qquad s.t. P=\gamma E\!\left\{I\!\left(X\right)\right\} \end{cases}\!\!\!\!\!, \end{align}

where $\mathcal{I}\;:\!=\;\{I|I\;:\; \mathbb{R}^{2}\mapsto [0, M]\}$ is the feasible set of measurable indemnity functions with an exogenous upper limit $M$ , $E$ denotes the expectation, $U$ is the concave utility function of the insured, and $w$ denotes initial wealth. A subsidy rate $\tau \in [0,1]$ is included in the utility function since the government often supports agricultural insurance, particularly in its pilot stage (Tadesse et al., 2015). Due to the assumption of a risk-neutral insurer, the participation constraint takes the form $P=\gamma E\{I(X)\}$ with a loading factor $\gamma \geq 1$ . To further simplify the application, we set $\gamma =1$ . In contrast to damage-based insurance, the usual constraint $I\lt Y$ does not apply. In fact, an indemnity payment is possible even if no damage occurs. The insurance premium $P\in [0, \gamma M]$ is assumed to be exogenously specified, for example, by the insurance buyer, and the indemnity payment is optimally adjusted to this premium.

Zhang et al. (2019) proved the existence and uniqueness of an optimal solution for the optimization problem (1). They emphasize the finding that the optimal indemnity is generally a nonlinear function of the weather index, though a linear payoff function is usually assumed in the literature. Typically, the optimal indemnity function has to be determined numerically, but there are closed-form solutions for particular utility functions. Specifically, for the quadratic utility $U_{qua}(r)=ar-br^{2}$ , with revenue $r\leq a/(2b)$ and parameters $a\gt 0$ and $b\gt 0$ , $I(X)$ takes the form:

(2) \begin{equation}I_{qua}^{\textrm{*}}\!\left(X\right)=\left[\left\{E\!\left(Y|X\right)+\eta _{qua}^{\textrm{*}}\right\}\vee 0\right]\wedge M,\end{equation}

where $E(Y|X)$ is the conditional expectation of $Y$ given $X$ .

For the exponential utility function $U\!\left(r\right)=-\frac{1}{\alpha }e^{-\alpha r}$ , with risk aversion parameter $\alpha \gt 0$ , the solution is

(3) \begin{equation}I_{exp}^{\textrm{*}}\!\left(X\right)=\left(\left[\frac{1}{\alpha }\log \left\{E\!\left(e^{\alpha Y}|X\right)\right\}+\eta _{exp}^{\textrm{*}}\right]\vee 0\right)\wedge M.\end{equation}

The constants $\eta ^{*}$ are determined by the constraint $E[I^{*}\!\left(X\right)]=P/\gamma$ . Equations (2) and (3) indicate that the optimal indemnity functions under quadratic or exponential utility are affected neither by the initial wealth $w$ nor subsidy rate $\tau$ . Moreover, it becomes apparent that the estimation of conditional expectations of yield loss, given the weather index, is pivotal to the practical implementation of this kind of weather insurance. In our empirical application, we will focus on the indemnity functions (2) and (3).

2.2. Statistical approach

In the optimal indemnity equation, the conditional expectation $E(Y|X)$ or $E(e^{\alpha Y}|X)$ is essentially a regression model of $Y$ or $e^{\alpha Y}$ , given weather index $X$ . As mentioned before, we employ GAM with B-spline row tensor product smoothers to estimate the conditional expectation and then embed it into the optimal indemnity function. We are not restricted to B-splines as GAM allows for flexibility in selecting smoothers, such as kernels or cubic splines (Hastie & Tibshirani, 1987). In this work, we use splines based on our preliminary analysis. Kernel estimation provides similar results, yet the models are often overparameterized. As introduced by Eilers & Marx (1996), P-spline, defined as the combination of B-spline with difference penalties on the estimated coefficients, has a valuable property that it shows no boundary effects, while many types of kernel smoothers do exhibit these. Aydin (2007) also concluded that smoothing spline regression estimators outperform the kernel ones. To the best of our knowledge, the applications of splines in agricultural insurance are rare, as mentioned in the Introduction.

To investigate the interaction between weather indices $X_{1}$ and $X_{2}$ and their effects on the response $Y$ , we construct a GAM including one smoother $S$ , that is, $Y=S(X_{1}, X_{2})$ , and denote the vectors of these three variables as $y$ , $x_{1}$ , and $x_{2}$ , respectively. Firstly, we obtain the marginal B-spline basis $B_{i,j}^{u}$ of order $j$ for each weather index through de-Boor recursion (de Boor, 1978):

(4) \begin{equation} B_{i,0}(u)=\left\{\begin{array}{l@{\quad}l} 1,& \textrm{if}\,u_{i}\leq u\lt u_{i+1}\\[5pt] 0, & \textrm{otherwise} \end{array}\right. \textrm{and} \end{equation}

(5) \begin{equation}B_{i,j}=\frac{u-u_{i}}{u_{i+j}-u_{i}}B_{i,j-1}\!\left(u\right)+\frac{u_{i+j+1}-u}{u_{i+j+1}-u_{i+1}}B_{i+1,j-1}\!\left(u\right),\end{equation}

where $u$ denotes the element in the weather index vector $x_{1}$ or $x_{2}$ , $\{u_{i}\}$ is a uniform knot vector determined by $k$ which is the number of equally spaced knots over the domain of $x_{1}$ or $x_{2}$ (strictly $k-1$ is the number of internal knots), and $i=0,\ldots, k+j-1,$ which represents the $i$ th basis function (Lee et al., 2013). Throughout this study, we use cubic splines ( $j=1,2,3$ ), which lead to functions with continuous second derivatives. Secondly, the smooth function $S(x_{1}, x_{2})$ is constructed by the B-spline tensor product $B$ and the GAM with one smoother is given by:

(6) \begin{equation}y=S\!\left(x_{1}, x_{2}\right)=B\theta +\varepsilon, \varepsilon \sim \mathcal{N}\!\left(0,\sigma ^{2}I\right),\end{equation}

where $B$ and $\theta$ are the “regression basis” and vector of “regression coefficients,” respectively (cf. Lee & Durbán, 2011) and $\varepsilon$ is a Gaussian error term with variance $\sigma ^{2}I$ . For efficient computation, row tensor product or box product $B=B_{2}\square B_{1}=(B_{2}\bigotimes 1_{c_{1}}^{\prime})*(1_{c_{2}}^{\prime}\bigotimes B_{1})$ instead of the Kronecker product $B_{2}\bigotimes B_{1}$ is used, where $B_{q}$ for $q=1,2$ denotes the marginal B-spline basis constructed from the weather indices $x_{q}$ by (4) and (5) and $c_{q}$ is a vector of ones of length $c_{q}$ which is the number of columns in $B_{q}$ (cf. Eilers et al., 2006; Lee et al., 2013). To control the overfitting and overparameterization, separate difference penalties are imposed on adjacent coefficients in $\theta$ along the two dimensions, which creates P-splines. For methodological details, we refer interested readers to Eilers & Marx (1996), Lee & Durbán (2011), and Marx & Eilers (1998). The number of knots $k$ lying between 20 and 40 is considered to be moderate (Ruppert et al., 2003).

We employ the P-spline ANOVA (PS-ANOVA) method to estimate the penalty terms. The PS-ANOVA method originated from the idea that through re-parameterizing the P-spline GAM defined above as (6) into a mixed model that contains both fixed and random effects (e.g., Brumback & Rice, 1998), the smoothing parameter becomes the ratio between the variance of the residuals and the variance of the random effect in a mixed model. Therefore, the smoothing parameter can be estimated using ANOVA decomposition (Currie & Durbán, 2002; Lee & Durbán, 2011). This transformed mixed model is

(7) \begin{equation}y=F\beta +Z{\delta} +\varepsilon,\; \varepsilon \sim \mathcal{N}(0,\sigma ^{2}I), \delta \sim \mathcal{N}\!\left(0, G\right), \end{equation}

where “ $\;:\;$ ” distinguishes the blocks in a matrix and the original box product $B=[F\;:\; Z]$ becomes a block matrix consisting of fixed effects matrix $F\equiv [1_{d}\;:\; x_{1}\;:\; x_{2}\;:\; x_{2}\square x_{1}]$ and random effects matrix $Z\equiv [Z_{1}\;:\; Z_{2}\;:\; Z_{2}\square x_{1}\;:\; x_{2}\square Z_{1}\;:\; Z_{2}\square Z_{1}]$ . $Z_{q}=B_{q}U_{q}$ for $q=1,2$ and $U_{q}$ are eigenvectors corresponding to the positive eigenvalues of the singular value decomposition of $D_{q}^{\prime}D_{q}$ . The original coefficient $\theta =(\beta,{\delta} )^{T}$ becomes a vector including the fixed effects coefficient $\beta$ and the random effects coefficient $\delta$ that is assumed to be Gaussian with covariance matrix $G$ . Note that the response variable in the mixed model is assumed to be Gaussian in this paper. Moreover, an isotropic penalization is conducted for simplification, which means that the same amount of smoothing is used for all covariates (Rodríguez-Álvarez et al., 2015). There are two main advantages of the PS-ANOVA method. First, one penalty for each covariate is attractive for tuning multiple penalty parameters, especially when the GAM includes several smooth functions (Lee et al., 2013). Second, the variance components of a mixed model can be estimated based on restricted maximum likelihood (REML), which is proved empirically to be preferable to other selection criteria, such as generalized cross-validation (GCV) or Akaike’s information criterion (AIC) (Schall, 1991; Reiss & Todd Ogden, 2009; Wood, 2011). ¹

3. Empirical Application: Designing Weather Index Insurance for Soybean Production

3.1. Data and study area

In our empirical study, we use the methodology described in section 2 to design a weather index insurance product for soybean producers in the US state of Illinois. Over one-third of soybeans traded on the global market are produced in the USA, and the value of US soybean exports reached a record of $27.4 billion in 2021 (USDA, 2022). Illinois contributed more than any other US state to soybean output in 2021 (The Illinois Soybean Association, 2022). Thus, significant disruption to soybean production in Illinois will not only affect local farmers but can also influence the international market substantially (Boyer et al., 2013). Soybean production is sensitive to weather conditions and specific weather events like hot-dry summer conditions, which can severely affect US soybean yields (Hamed et al., 2021). Insurance policies are available that cover yield losses for US soybean producers, such as the Supplemental Coverage Option (SCO) and Enhanced Coverage Option (ECO) (Schnitkey et al., 2022). Nevertheless, it appears reasonable to design a weather index insurance product to complement existing contracts given the current situation in which Rainfall Index Insurance is a single-peril insurance that only covers apiculture, annual forage, pasture, and rangeland (USDA, 2021).

The data used in our empirical analysis are publicly available, including yield, meteorological, and phenological data. County-level soybean yield data from 96 counties in Illinois from 1998 to 2020 are from the National Agricultural Statistics Services (NASS) (Fig. 1, upper panel). We detrend raw yield data by a linear time regression model to obtain detrended yields (Shi & Jiang, 2016) (Fig. 1, middle panel). Further, following Zhang et al. (2019), yield losses ² $Y$ in each county are calculated as the difference between the highest detrended yield during the dataset’s 23 years and detrended yields (Fig. 1, bottom panel). The average raw yield is about 48.46 bushels/acre, and the average detrended yield is approximately 38.46 bushels/acre. Yield losses range from 0 to 35.17, and the average yield loss amounts to 9.53 bushels/acre.

Figure 1 Yield (upper panel), detrended yield (middle panel), and yield loss (bottom panel) for 96 counties in Illinois. Each line represents one county.

Since temperature and precipitation are major determinants of plant growth (e.g., Walter, 1985; Donoghue, 2008), we use GDD and CR as weather-related indices in our analysis. We confine our analysis to a two-dimensional case, since this research aims to investigate the interaction between temperature and precipitation and their co-impacts on yield losses in different growth stages. We average the gridded daily weather data collected from the Daymet data set (Thornton et al., 2020) for the weather indices calculation in each county. $GDD$ and $CR$ are defined as:

(8) \begin{equation}GDD_{m,l}={\sum}_{t=sd_{m}}^{ed_{m}}\max \left\{0,\frac{\left(TMA_{m,t,l}+TMI_{m,t,l}\right)}{2}-TB\right\}\textrm{and} \end{equation}

(9) \begin{equation}CR_{m,l}={\sum}_{t=sd_{m}}^{ed_{m}}R_{m,t,l},\end{equation}

where $m$ and $l$ denote the years and the region (county), respectively, $t$ is the $t$ -th day in a year, $sd$ and $ed$ denote the start and the end date of the accumulation period, respectively, and $TMA$ , $TMI,$ and $R$ are daily maximum temperature, minimum temperature, and rainfall, respectively. The baseline temperature $TB$ of soybean growth is defined as 10° Celsius (Scholtes et al., 2019). GDD and CR indices have been criticized because they do not account for the temporal distribution of weather events within the accumulation period. This weakness, however, is attenuated because we divide the accumulation period into shorter periods. In addition, we implement a Rainfall Deficit Index (RDI) that has been suggested as an alternative to CR (e.g., Odening et al., 2007). A definition of this index as well as model results are presented in the Appendix.

Agronomic research asserts that the weather–yield relationship is time-variant, that is, meteorological factors affect plant development differently during each growth stage (e.g., Delerce et al., 2016). Thus, we divide the whole growth cycle of soybeans into four phases. According to the terms and definitions of the NASS (2018), these four phases are described as “emerged,” “blooming,” “setting pods,” and “dropping leaves.” The phase-division procedure is based on the weekly state-level Crop Progress Report (USDA-Economics, Statistics and Market Information System (ESMIS), 2022). For the first stage, “emerged,” for example, we selected the dates for which germination rates are closest to 1% and 100% as the start and end dates, respectively, for this phase. As a result, these dates may change slightly from year to year. For example, in 2020, the beginning and end dates of the first phase are May 3 and June 28, respectively, which deviate about 1 week from those in 1998. Moreover, temporal overlapping between two subsequent growth stages can occur. For example, “blooming” started on June 28, 1998, before the complete germination of soybeans ended. To explore the performance gain from dividing the entire vegetation period, we also estimate a benchmark model without decomposition whose beginning and end dates of the accumulation period are defined as the start date of the first growth phase and the end date of the last growth phase, respectively. In our study, we obtained the stage division points from the Crop Progress Report, which was only available after all the necessary information had been collected. Unfortunately, this means that the direct application of insurance contracts is not feasible because we do not know the division points for future years that will be insured. The challenge is compounded by our limited understanding of the underlying mechanism behind plant phenology, as discussed by Tang et al. (2016). Given these uncertainties, we opted to rely on ground observations as the basis for our division procedures.

Figs. 2, 3, and 4 depict the $GDD$ and $CR$ values of 96 counties over the observation period for the entire growth cycle at each growth stage. In Figs. 2, 3, and 4, the values are annual and accumulative, and each line is smoothed by observation points across 23 years. Each subplot reflects accumulative values at a particular phase.

Figure 2 Time series of growing degree days (GDD) and cumulative rainfall (CR) in the whole cycle. Each line represents one county.

Figure 3 Time series of growing degree days (GDD) in separate phases. Each line represents one county.

For example, in Fig. 3, upper left plot, each line represents the GDD values in the “Emerged” stage across 23 years of one specific county, where stage 1 of the first year is immediately followed by stage 1 of the second year. Some facts are notable. First, both weather indices show a considerable variation over time, reflecting the prevalence of weather risk. Second, the graphs reveal spatial heterogeneity of weather conditions across Illinois, though all counties follow similar patterns. Third, unsurprisingly, spatial correlation is less pronounced for rainfall than for temperature. Fourth, the patterns for the two weather indices differ across the four growth stages, which motivates our proposed division. Actually, though the observation period is rather short, some changes in weather patterns can be observed which may be interpreted as a manifestation of climate change. There a several aspects to be mentioned in this context. First of all, it is essential to distinguish between changes in levels (e.g., “global warming”) and changes in weather variability. From an insurance perspective, the latter is more important. The question of whether weather becomes more volatile over time has been intensively discussed in the literature (e.g., Wang et al., 2013). Accordingly, models have been developed that capture time-varying volatility of temperature and rainfall processes (e.g., Okhrin et al., 2013). A related question is whether yields become more volatile due to climate change and what amount of historical data should be used to estimate models (e.g., Shen et al., 2018). Admittedly, our modeling approach does not account for changing parameters of the stochastic weather processes. However, this seems acceptable, because our paper focuses on the optimal design insurance contract that depends on the weather–yield relationship and this relationship is not directly affected by climate change.

3.2. Model estimation

We determine the optimal indemnity function under quadratic and exponential utility functions, which requires estimating conditional expectations of yield loss $Y$ and $e^{\alpha Y}$ , respectively, given the weather indices $X$ and a specific combination of $P$ and $M$ according to equations (2) and (3), respectively. The estimation of the conditional expectation is conducted for two settings. First, we estimate conditional yield loss based on the weather–yield loss relation for the entire growth period, for example, $Y_{m,l}=S_{0}(GDD_{m,0,l}, CR_{m,0,l})+\varepsilon _{m,0,l}$ in the quadratic utility case, where the subscript zero denotes the entire growth cycle. This GAM serves as a benchmark in our analysis. Second, we divide the whole growth period into four phases and estimate the weather–yield loss relation for the phase-division GAM consisting of four interaction smoothers, namely:

(10) \begin{align}\textrm{Y}_{m,l} &=S_{1}\!\left(GDD_{m,1,l},CR_{m,1,l}\right)+S_{2}\!\left(GDD_{m,2,l},CR_{m,2,l}\right)+S_{3}\!\left(GDD_{m,3,l}, CR_{m,3,l}\right)\nonumber \\[5pt] &\quad +S_{4}\!\left(GDD_{m,4,l}, CR_{m,4,l}\right)+\varepsilon'_{m,s,l},\end{align}

where the subscripts $1$ , $2$ , $3$ , and $4$ represent the phases of “emerged,” “blooming,” “setting pods,” and “dropping leaves,” respectively, and $s$ is an overall symbol denoting phase division. Moreover, to consider interactions among weather variables within different stages, we attempted to model the interactions between all four phases (one eight-dimensional smoother). Unfortunately, it was computationally highly demanding. Therefore, we allow for interactions within the four variables of the first two phases, and within the four variables of the third and the fourth phase, that is, a two-segment model containing two four-dimensional smoothers as a benchmark model:

(11) \begin{align}\textrm{Y}_{m,l}&=S_{\textrm{I}}\!\left(GDD_{m,1,l}, CR_{m,1,l},GDD_{m,2,l}, CR_{m,2,l}\right)+S_{\textrm{II}}\!\left(GDD_{m,3,l}, CR_{m,3,l},GDD_{m,4,l}, CR_{m,4,l}\right) \nonumber \\[5pt] & \quad +\varepsilon''_{m,s,l}.\end{align}

Only the benchmark method mgcv::gam( ) is used for the estimation of the model (11) due to computational burden. There might be more efficient ways to handle high-dimensional weather variables, such as neural network models (Crane-Droesch, 2018; Chen et al., 2020) or regression tree-based models. However, in our research, we specifically chose the GAM to prioritize interpretability for each separate growth phase. While neural networks may offer potential improvements in model performance, the relatively small dimensions of our study (only eight variables) suggest that nonparametric regression, as employed in the components of GAM, provides qualitatively similar results. Additionally, neural networks often require extensive hyperparameter tuning, which can be time-consuming and challenging. Thus, the investigation of the performance of other methods in phase-division crop yield modeling is beyond the scope of this study.

In the case of exponential utility, the response variable of GAM is $e^{\alpha Y}$ , where $\alpha-$ is the (absolute) risk aversion coefficient. In our analysis, we consider three levels of risk aversion for US soybean farmers: low-, moderate-, and high-risk aversion. To determine appropriate levels for the risk aversion coefficient, we follow Tan & Zhang (2020), who derived $\alpha-$ by dividing a relative risk aversion coefficient of 2, 3, and 4, respectively, by an estimate of the initial wealth of corn producers in Illinois. This results in the following absolute risk aversion coefficients $\alpha- =0.0052$ (low), $0.008$ (moderate), and $0.0103$ (high).

Figure 4 Time series of cumulative rainfall (CR) in separate phases. Each line represents one county.

To circumvent potential overfitting, cross-validation is conducted to select the hyperparameter $k$ , that is, the number of knots. To this end, the entire data set was divided into three. The training data set contains data from 1998 to 2013, and the validation data set includes data from 2014 to 2017, while the test data set contains data from 2018 to 2020. For computational stability, we normalize the values of all variables into the range $[0,1]$ . After estimation, we convert the fitted values back to the original scale to obtain the scale-dependent root mean squared error: $RMSE=\sqrt{\frac{\sum _{m=1}^{n}\sum _{l=1}^{96}\left(\widehat{y_{m,l}}-y_{m,l}\right)^{2}}{96*n}}$ , where $n$ is the time series length of each county in the data set, for example, $n=16$ for the training set. $RMSE$ measures the average deviation between predicted and observed values and serves as our study’s model evaluation criterion (Schmidt et al., 2022). We select the appropriate $k$ values in the range $[20,40]$ by cross-validation. For whole-cycle models, cross-validation resulted in $k=21$ for the one under quadratic utility and the one under exponential utility given $\alpha- =0.0052 ;\; k=30$ for the models under exponential utility given $\alpha- =0.008$ and $\alpha- =0.0103$ . For all phase-division models, cross-validation resulted in $k=39$ . After cross-validation, we estimate the models with the full sample size from 1998 to 2020 using selected $k$ values.

3.3. Results from the weather–yield loss relation

Table 1 shows that in all cases, the models based on phase division significantly outperform the whole-cycle models in terms of $RMSE$ and adjusted $R^{2}$ . ³

Table 1. $RMSE$ and adjusted $R^{2}$ of GDD-CR GAMs estimated by PS-ANOVA & mgcv::gam( ).

Notes: RMSE is the root mean squared error.

Compared with the benchmark method mgcv::gam( ), PS-ANOVA shows a slightly better fit. $RMSE$ of the two-segment model (11) is 3.79 bushels/acre, and the Adj. $R^{2}$ is 0.60, indicating performance gain compared to the original four-phase-division model (10). There is clearly a trade-off between model performance and model parsimony. Despite this, model (10) effectively isolates the impact of each phase on yield losses, as we will see later. In contrast, the interpretability of the two-segment model (11) is limited not only by the arbitrary nature of segment numbers – as we could theoretically conduct a three-segment model containing two three-dimensional and one two-dimensional tensor products – but also by the challenges of visualizing interactions between more than three variables. As a result, the specific effects of each phase may not be clear. Through our analysis and the results presented in the figures and the discussion that we will see in section 3.3, we could clearly separate and quantify the contributions of each phase to yield losses. This interpretability is crucial for understanding the underlying factors affecting crop yield and designing effective index insurance.

A contour plot of the weather–yield loss relation is shown in Fig. 5 for the whole-cycle model. It visualizes the complex interaction of variables $GDD$ and $CR$ as determinants of yield losses in soybean production. Overall, soybean yield loss decreases with increasing $GDD$ because higher temperatures induce a faster crop development provided that a sufficient amount of rainfall is available, in this case, more than 500 mm (Peiris et al., 1996). When the $GDD$ is as high as 2,200° Celsius and the $CR$ amounts to approximately 600 mm, yield loss is at its minimum (about four bushels/acre), which is less than half of the average yield loss (9.53 bushels/acre). However, when the $CR$ is below 500 mm, the contour lines appear symmetric. It means that under a deficit in rainfall, both insufficient and excessive $GDD$ cause yield losses. Specifically, when $GDD$ exceeds the threshold of 1,800° Celsius, the joint occurrence of rainfall deficiency and heat stress causes significant yield loss. Thus, our analysis confirms the negative impact of hot-dry conditions on US soybean, which has been found in previous studies (e.g., Hamed et al., 2021). Furthermore, Fig. 5 depicts that excessive rainfall also induces yield losses, particularly under cold conditions. However, soybean yield losses are less sensitive to excessive rainfall than drought.

Figure 5 Whole-cycle weather–yield loss relation. Notes: CR is cumulative rainfall and is given in millimeters (mm). GDD is growing degree days and is given in degrees Celsius. The subscript 0 after CR and GDD indicates that the whole-cycle model is utilized rather than the models that separate the four growing phases.

Fig. 6 depicts the impact of temperature and rainfall on soybean yield losses in each growth phase. In each subplot, the displayed values measure the contribution to total yield loss, where negative values indicate a reduction of yield losses. Comparing the four subplots reveals that the weather–yield loss relation differs significantly between the four growth phases. The level of the loss contribution, sensitivity to weather events, and the interaction between temperature and precipitation show pronounced differences. In the first phase (“emerged”), small losses occur at the $GDD$ - $CR$ combination of around $(650^{\circ}\textrm{ Celsius},200\textrm{ mm}),$ whereas large losses take place when $CR$ is lower than 100 mm, and $GDD$ is below 550° Celsius. Both $CR$ and $GDD$ influence yield production positively within some range but result in a negative impact after exceeding a certain level. In 2012, for example, the $CR$ of 32 counties added up to less than 50 mm in this development phase. As a result of this countrywide drought, a considerable production disruption was reported, with an estimated soybean yield loss of around 170 million bushels (Boyer et al., 2013; Rippey, 2015). The contour plots of the second phase (“blooming”) show that potential yield losses are comparably large, that is, soybean plants are vulnerable in this growth phase. Minimal yield losses are realized at 650° Celsius and 150 mm. Notably, plants are more sensitive to excessive rainfall than drought in this phase. High temperatures advance blooming (Cooper, 2003), and more precipitation will likely cause a longer blooming period and high outcrossing rate (Qu et al., 2020). Another possible explanation is that under persistent wet weather, the disease sclerotinia stem rot could attack soybeans during reproductive stages and, in turn, cause yield loss (Wrather & Koenning, 2009).

Figure 6 Phase-division weather–yield loss relation. Notes: Each subplot represents the yield loss contribution from the interaction between growing degree days (GDD) and cumulative rainfall (CR) at the corresponding growth phase. The four growing phases are shown in the subscripts on the labels of the axes: 1 indicates emerged, 2 blooming, 3 setting pods, and 4 dropping leaves. CR is in mm and GDD is in degrees Celsius.

The subplot of the third stage (“setting pods”) suggests a pattern that can be summarized as “the warmer and wetter, the better.” When $GDD$ is in the range $[400,650]$ , yield loss remains unchanged if the $CR$ level exceeds 125 mm. This indicates that temperature conditions restrict the soybean growth rate. Once $GDD$ exceeds 650° Celsius, more rainfall means higher yield gain. The plot also shows that $GDD$ affects yield loss only slightly under rainfall deficit conditions. This can be explained by pods’ abortion under drought stress, since sufficient water is critical to fill the soybean pods (Liu et al., 2004). In the fourth phase (“dropping leaves”), the soybean yield is relatively robust against varying rainfall. For instance, when $GDD$ lies in the range $[175,500]$ and $CR$ is above 150 mm, the contour lines are approximately vertical, which implies that rainfall increments have little impact on yield loss. This result reflects the fact that the plant’s demand for water declines when reaching maturity (Jensen, 1968). On the other hand, soybean’s temperature demand seems to be relatively strong until $GDD$ achieves 500° Celsius. The bottom right area of the plot visualizes that yield losses begin to increase under the hot-dry weather. One possible interpretation for this phenomenon is that drought stress increases crop vulnerability to aphids, a transmitter of the Soybean Mosaic Virus, and further suppresses soybean yields (Rice et al., 2007; van Munster et al., 2017).

3.4. Hedging effectiveness

Though estimating the conditional loss expectation for different plant growth phases provides exciting insights into the weather–yield relationship, it is not yet clear if the suggested procedure can improve the design of weather index insurance. To answer this question, we calculate the hedging effectiveness of a phase-division contract and a standard contract based on the entire vegetation period. Hedging effectiveness is defined in terms of the expected utility of the insurance contract. We determine the optimal indemnity functions under the assumption of an exponential utility function and parametrize the level of risk aversion. Under exponential utility, $EU$ is $E\!\left(-\frac{1}{\alpha }e^{-\alpha r}\right)$ and revenue is $r=w+I^{*}(X)-Y-(1-\tau )P$ after the indemnity payment. As we are interested in the relative performance of the phase-division contract compared with the standard contract, we determine the ratio of the expected utility of the two contracts. Thus, the initial wealth $w$ , premium $P,$ and subsidy rate $\tau$ are canceled out, and the expectation of the response variable $e^{\alpha Y}$ is equivalent to the product of the fixed effects matrix and fixed effects coefficient, that is, $E(e^{\alpha Y})=F\beta$ given the zero-mean assumption of the random effect $\delta$ and the model error term $\varepsilon$ in the mixed model (7). We conduct the calculation for each county, that is, we consider each county as a representative farmer in Illinois.

Up to this point, estimating a distribution of weather variables was not required since calculating the optimal indemnity function is based on the conditional loss expectation given the weather conditions. To arrive at the expected utility, we simply average the observed realizations of weather variables in the 23 years in the observation period. With this simplified historical simulation, the $EU$ ratio under exponential utility takes the following form:

(12) \begin{equation}\frac{\widehat{EU_{s,l}}}{\widehat{EU_{0,l}}}=\frac{1}{23}{\sum}_{m=1}^{23}\frac{e^{-\alpha \widehat{I_{m,s,l}^{\textrm{*}}}\left(X_{m,s,l}\right)}F_{m,s,l}\widehat{\beta _{m,s,l}}}{e^{-\alpha \widehat{I_{m,0,l}^{\textrm{*}}}\left(X_{m,0,l}\right)}F_{m,0,l}\widehat{\beta _{m,0,l}}},\end{equation}

where the subscripts $s$ and zero denote “phase-division” and “whole-cycle,” respectively, as defined in section 3.2, and $m$ and $l$ are year and county indices, respectively. Considering that the average yield loss is about 9.53 bushels/acre and the maximum yield loss is approximately 35.17 bushels/acre, we assume that the premium $P$ is 9 $/acre and the maximum indemnity $M$ is $35 per acre, which normalizes the soybean price to $1/bushel without the loss of generality. Under these assumptions, the optimal indemnities $I^{*}(X)$ of phase-division and whole-cycle contracts are obtained according to equation (3). The 96 counties are clustered into nine groups according to the agricultural district attributes defined by the NASS (2022). The nine groups are “Northwest” (10), “Northeast” (20), “West” (30), “Central” (40), “East” (50), “West Southwest” (60), “East Southeast” (70), “Southwest” (80), and “Southeast” (90). Fig. 7 depicts which counties are included in these groups and where they are located.

Figure 7 County-level agricultural district classification in Illinois.

Fig. 8 presents the results of the hedging effectiveness calculation. We find that the $EU$ ratios are above one across all nine regions, which indicates that the phase-division contract outperforms the whole-cycle one. However, the gain is modest. For instance, when $\alpha =0.008$ , the mean of all the 96 $EU$ ratios is about 1.017, that is, on average, the expected utility of the phase-division contract is only 1.7% higher compared to the whole-cycle contract. For more risk-averse policyholders ( $\alpha =0.0103$ ), the relative superiority of the phase-division contract increases: the mean of the $EU$ ratios is about 1.022. Furthermore, we observe that $EU$ ratios vary across agricultural districts. For example, in the East District (50), the $EU$ ratio is considerably lower than in the Southwest district (80) for all three levels of risk aversion. In addition, we implement the mean root square loss (MRSL) (e.g., Vedenov & Barnett, 2004) and variance of revenue (VAR) as alternative performance measures and calculate its percentage change with and without insurance for both the whole-cycle and the phase-division contract.

Figure 8 $EU$ ratio under three levels of risk aversion.

For each county, the MRSL, VAR, and the percentage changes of MRSL and VAR after insurance contract purchasing are defined as:

(13) \begin{equation}MRSL=\sqrt{\frac{1}{23}{\sum}_{m=1}^{23}\left\{\max \left(0, \textit{price}\;\ast\;\overline{\textit{yield}}-r_{m}\right)\right\}^{2}},\end{equation}

(14) \begin{equation}VAR=\sqrt{\frac{1}{23}{\sum}_{m=1}^{23}\left(\textit{price}\;\ast\;\overline{\textit{yield}}-r_{m}\right)^{2}},\end{equation}

(15) \begin{equation}\textit{MRSLPC}=\frac{MRSL_{with}-MRSL_{\textit{without}}}{MRSL_{with}},\end{equation}

(16) \begin{equation}\textit{VARPC}=\frac{VAR_{with}-VAR_{\textit{without}}}{VAR_{with}},\end{equation}

where $\textit{price}$ is set to be $1/bushel to represent the soybean price; average yield across 23 years is $\overline{\textit{yield}} ;\; m$ is the year, and $r_{m}$ is the yearly revenue of each county. In the case of “with” contract, the revenue $r_{m}=\textit{price}\;\ast\;\overline{\textit{dyield}_{m}}$ , where $\overline{\textit{dyield}_{m}}$ is detrended yield in each year. In the case of “without” contract, the revenue $r_{m}=\textit{price}\;\ast\;\overline{\textit{dyield}_{m}}+I_{exp,m}^{*}(X_{m})-P$ , where the premium $P=9\$ /acre$ is equivalent to expected payoff in this paper.

Averages of the MRSLPC and VARPC of GDD-CRI models are provided in Table 2. The results show a significantly better performance of the phase-division contracts. For example, given $\alpha- =0.0103$ , the average decrease in MRSL after purchasing the whole-cycle contract is 18.4%, while the average decline in MRSL for a phase-division contract is as much as 41.9%. Moreover, the average decrease in variance after purchasing the whole-cycle contract is 16.6% given $\alpha- =0.0103$ , and the average decrease of the phase-division one is 37.8%.

Table 2. Average $\textit{MRSLPC}$ and $\textit{VARPC}$ of GDD-CRI models.

4. Conclusions

This paper utilizes a nonlinear indemnity function framework to optimize weather index insurance. We suggest a GAM with P-splines box product smoothers to estimate conditional yield loss functions as a flexible alternative to commonly used regression models or neural networks that recently became popular for modeling weather–yield relations. The statistical framework is applied to the case study of soybean production in the US state of Illinois. Rainfall and temperature and their complex interaction are the main drivers of soybean yield losses. Our analysis contributes to decomposing the whole growth cycle into four stages. Our results indicate that the proposed phase-division contract outperforms the whole-cycle one regarding model fit and hedging effectiveness. This finding is in line with earlier studies emphasizing informational gains from disaggregated weather data (e.g., Schmidt et al., 2022). The division of the growth cycle allows a state-dependent analysis of the collective impact of temperature and rainfall on crop yields. This is not only beneficial for the design of weather index insurance, but it may also support managerial decisions, such as the timing of irrigation measures (Katyal & Pandian, 2019). Moreover, the model performance gain is significant in terms of the percentage change of mean root square loss and VARs despite the relatively minor outperformance in EU ratio. The extent to which our results can be generalized to other regions, crops, and weather variables is suggested for future research.

We must address some limitations of our analysis. First, aggregated county-level soybean yield data are used instead of individual farm-level yield data due to data availability, which inevitably underestimates yield variability and basis risk of weather index insurance (Popp et al., 2005). We conjecture that the relative advantage of the proposed statistical procedure will increase if it were applied to farm-level data with higher variability. Second, weekly state-level data from the Crop Progress Report cause temporal overlapping of plant growth phases. This dilutes the phase division and may result in an estimation error. Finally, the distribution of weather variables and the joint density of yield losses and weather indices have not been estimated, and thus the random nature of weather perils has not been explicitly considered. While this is not necessary for calculating the indemnity payments in our setting, it is essential for evaluating weather-related insurance products. Linking the proposed P-spline method with a parametric or nonparametric estimate of the density of weather indices would be an exciting direction for future research.