Data-driven determination of plant growth stages for improved weather index insurance design

Abstract

Purpose

This paper aims to improve the delimitation of plant growth stages in the context of weather index insurance design. We propose a data-driven phase division that minimizes estimation errors in the weather-yield relationship and investigate whether it can substitute an expert-based determination of plant growth phases. We combine this procedure with various statistical and machine learning estimation methods and compare their performance.

Design/methodology/approach

Using the example of winter barley, we divide the complete growth cycle into four sub-phases based on phenology reports and expert instructions and evaluate all combinations of start and end points of the various growth stages by their estimation errors of the respective yield models. Some of the most commonly used statistical and machine learning methods are employed to model the weather-yield relationship with each selected method we applied.

Findings

Our results confirm that the fit of crop-yield models can be improved by disaggregation of the vegetation period. Moreover, we find that the data-driven approach leads to similar division points as the expert-based approach. Regarding the statistical model, in terms of yield model prediction accuracy, Support Vector Machine ranks first and Polynomial Regression last; however, the performance across different methods exhibits only minor differences.

Originality/value

This research addresses the challenge of separating plant growth stages when phenology information is unavailable. Moreover, it evaluates the performance of statistical and machine learning methods in the context of crop yield prediction. The suggested phase-division in conjunction with advanced statistical methods offers promising avenues for improving weather index insurance design.

Citation

Zou, J., Odening, M. and Okhrin, O. (2024), "Data-driven determination of plant growth stages for improved weather index insurance design", Agricultural Finance Review, Vol. 84 No. 4/5, pp. 297-319. https://doi.org/10.1108/AFR-01-2024-0015

1. Introduction

Basis risk, namely the difference between the triggered insurance payout and the actual yield losses, poses a significant hurdle to the wide adoption of weather index insurance (e.g. Woodard and Garcia, 2008; Odening and Shen, 2014; Janzen et al., 2021). Thus, much attention has been paid to quantifying and alleviating this challenge (see Norton et al., 2012; Keller and Saitone, 2022). Basis risk can be decomposed into spatial and product risks (see Boyd et al., 2019).

Spatial basis risk arises from variations between on-site farmland weather conditions and weather indices measured at a reference weather station. These differences can be mitigated through the application of spatial interpolation methods (e.g. Cao et al., 2015; Boyd et al., 2019; Leppert et al., 2021) or the diversification of insurance contracts across different regions (Ritter et al., 2014).

Product basis risk, on the other hand, stems from the inability of weather indices to comprehensively account for variability in agricultural yields. Product basis risk reflects the fact that weather conditions are pivotal in crop production, but are not the sole determinant. Other factors, such as pests, controllable inputs (e.g. fertilizer), and managerial practices affect crop yields as well. Albers et al. (2017), for example, estimated that only 43% of the variability in wheat yield can be attributed to weather conditions. Despite this inherent limitation in reducing product basis risk, many attempts have been made to optimize the design of weather index insurance. Two main strands of literature can be distinguished. First, the specification and estimation of the relationship among yield and weather indices. Second, the selection and specification of appropriate weather indices.

The first strand of literature targets reducing product basis risk by estimating the weather-yield relationship with sophisticated statistical modeling. This task is challenging since the functional form of weather-yield relation is not only region- and crop-specific, but, in general, highly nonlinear (e.g. Schlenker and Roberts, 2009). In early works, Polynomial Regression (PR) has been used to model the nonlinear weather-yield relation (e.g. Vedenov and Barnett, 2004). PR is relatively simple in implementation and easy to interpret, however, the resulting model fit is often poor (see Schmidt et al., 2022). In addition, traditional nonparametric regression methods, such as the Generalized Additive Model (GAM) based on B-spline or penalized B-spline (P-spline) smoothers, have also been employed for modeling the nonlinear relationship between weather and yield with a very good quality (cf. Bucheli et al., 2022; Tan and Zhang, 2023; Zou et al., 2023). In addition to conventional statistical methods, more recently, machine learning methods have been applied to weather-yield modeling (Trawiński et al., 2012). In quantifying the effects of weather on crop yield, Random Forest (RF), a tree-based ensemble method (Breiman, 2001), turns out to achieve higher prediction capability than PR (e.g. Feng et al., 2018; Hoffman et al., 2020). Support Vector Machines (SVM) are also regarded as a powerful regression method in crop yield prediction (e.g. Oguntunde et al., 2018; Dang et al., 2021), though they are primarily applied in field crop classification tasks (e.g. Mathur and Foody, 2008; Löw et al., 2013). Moreover, Artificial Neural Networks (ANN) have been used to capture the relationship between high-dimensional weather variables and yield in large data sets (Schmidt et al., 2022; Chen et al., 2023). So far, there is no consensus on the “best” approach for modeling the weather-yield relation and complexity in methods does not always lead to superior performance (Oikonomidis et al., 2022). Some model comparisons have been conducted (e.g. Das et al., 2020; Joshi et al., 2021), but there is no exhaustive comparison of model performance between statistical and machine learning techniques. In addition, machine learning methods are often criticized by their black box properties that imply lack of interpretability, which could probably hinder the promotion of index insurance products among policy holders since contract transparency is often required (Sibiko et al., 2018). To overcome this drawback, visualization with partial dependence plots is considered to increase the explanatory capabilities of machine learning methods (Greenwell, 2017), which further enables the comparison between methods in terms of the interpretation of results.

The second strand of literature about mitigating product basis risk involves choosing generic weather variables and designing weather indices. There is a consensus that precipitation and temperature are pivotal in crop production. Weather indices are usually calculated by aggregating weather variables (e.g. cumulative rainfall or average temperature). Often used weather indices are Growing Degree Days (GDD), Cumulative Rainfall (CR), and the standardized precipitation index (cf. Turvey, 2001; Stoppa and Hess, 2003; Hill et al., 2019; Okpara et al., 2017). Moreover, remote sensing data are used to retrieve vegetation indices as alternatives to meteorological indices in the design of index insurance (cf. Möllmann et al., 2019; Vroege et al., 2021). Important design parameters are the beginning and length of the accumulation period and the aggregation level of weather indices. Turvey et al. (2021) select the start and end months of a crop growing season from all combination candidates based on model performance. If the accumulation period of the weather indices is not aligned with the vegetation period of the insured crop, a so-called temporal basis risk may arise (Conradt et al., 2015a), which can be considered a special case of product basis risk. Moreover, the varying water demand of crops and vulnerability against weather events during different plant growth stages are often ignored (Deng et al., 2007). To address these problems and thus mitigate temporal basis risk, it has been suggested to split the entire vegetation period into distinct plant growth stages and to calculate weather indices for each growth stage separately (e.g. Zou et al., 2023). Usually, the determination of beginning and end dates of growth phases is conducted by experts based on phenology information. However, the required information is not always available, at least not at the desired spatial resolution. This deficit calls for a data-driven division procedure as an alternative.

Against this backdrop, the contribution of our paper is twofold. First, we propose a data-driven method to search for the optimal division points of plant growth stages as an alternative to an expert-based growth stage division. This is accomplished via the minimization of weather-yield relation estimation errors. We investigate the feasibility as well as the costs and benefits of this approach. We find that the benefits are modest in terms of improvements in model fit. However, model development cost and computational costs are outweighed in situations where external information on division timings does not exist. Second, we provide a thorough comparison of classical statistical methods with various machine learning techniques in the context of yield prediction with weather variables. The comparison not only focuses on model fit and estimation errors, but also includes the visualization of the partial dependence of crop yields on individual weather variables in different growth phases. We use relatively simple, commonly used weather indices because our study focuses on a comparison of statistical methods and the optimal decomposition of the vegetation period.

The remainder of this paper is organized as follows: Section 2 motivates the selected statistical and machine learning techniques with a simulation study illustrating the strengths and weaknesses of each method; Section 3 displays the data and study region of our empirical analysis; Section 4 describes model estimation including parameter tuning procedures, shows numerical results, and visualizes and interprets the impact of individual predictors on response variable by means of partial dependence plots; and Section 5 concludes.

2. Statistical methods and machine learning techniques for weather-yield modeling

2.1 Method selection

In this study, we implement five methods for modeling crop weather-yield relation: Polynomial Regression (PR), Generalized Additive Model (GAM), Random Forest (RF), Support Vector Machines (SVM), and Artificial Neural Networks (ANN). The model choice is based on the idea that the methods should represent different model families, specifically statistical methods and machine learning techniques, and should also complement each other. Moreover, the selection considers model performance in previous analyses. PR is a parametric linear model and serves as the complementary benchmark in this study for the evaluation of the machine learning methods. A penalization technique, namely the Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996), is excluded from the regression since feature selection by reducing insignificant redundant variables does not apply to our case given at most eight input weather indices. Note that although PR is used to estimate non-linear relationships with quadratic terms, it is still a linear regression given the linearity in regression coefficients (Hastie et al., 2009a). GAM is a family member of generalized linear models (Hastie and Tibshirani, 1985), which possesses the interpretability of linear models but in a non-parametric setting (Fan and Jiang, 2005) and exhibits good performance (Zou et al., 2023). In our study, the smoothers in GAM are constructed by two-dimensional P-splines, which are considered as a state-of-art statistical method in yield modeling (cf. Tan and Zhang, 2023). The other three models, RF, SVM and ANN, are classical tools in machine learning and have been proven to be able to model general non-linear relationships. RF is a representative of tree-based methods and, turns out to be a robust prediction method compared to many other tree-based models. For example, Extreme Gradient Boosting (XGBoost) proposed by Chen and Guestrin (2016) is prone to overfitting when the sample size of training data is small as pointed out by Folberth et al. (2019). Considering SVM’s excellent capability in training a small-sized data set (Mountrakis et al., 2011) and that yield data are usually scarce (Odening and Shen, 2014), SVM is selected from kernel-based learning methods to solve the regression task in this study. Finally, ANN offers a comparative machine learning technique to RF and SVM, although it is more widely implemented in solving big data problems (e.g. Pasini, 2015). In our application, a pair of indices WI1 and WI2 are calculated simultaneously in each index aggregation period, such as the accumulative precipitation and average temperature in a certain month, during one plant growth stage or across the whole growth cycle. Since this study highlights the phase-division procedure, we do not dive into the selection of weather indices, which can be investigated in further research. We estimate the expected values of yield variables although the estimation of the lower bound of the yield distribution might be useful in an insurance context (e.g. Conradt et al., 2015b). Moreover, in a decision theoretic framework, e.g. expected utility maximization, information about the entire yield distribution is required (e.g. Zhang et al., 2019). Likewise, the decomposition of basis risk into upside and downside basis risk hinges on the entire yield distribution (cf. Schmidt et al., 2022). In what follows, we describe the specification of the five selected methods in greater detail.

2.1.1 Polynomial regression (PR)

PR is a widely used statistical model that represents a linear combination of predictors, their powers, and cross-products, offering the advantages of simplicity and interpretability. For example, Vedenov and Barnett (2004) apply PR to model the nonlinear relationships between weather indices and corn, cotton, and soybean yields. Chen et al. (2017) design rainfall index insurance using PR. In this paper, we employ quadratic regression to estimate nonlinear relations as a benchmark following Schmidt et al. (2022), where PR fails to capture the yield variability as substantially as the machine learning methods. To illustrate this, we assume only one index aggregation period. For example, in case WI1 and WI2 are two different weather indices calculated across the whole growth cycle, the PR model takes the following form:

(1)Ym,l=β+β1WI1m,l+β2WI2m,l+β3WI1m,l2+β4WI2m,l2+β5WI1m,l∙WI2m,l+ϵm,l,

where Ym,l is yield-related response variable, m represents the year, l denotes the farm, β is the intercept, and β1, …, β5, are coefficients estimated by Ordinary Least Squares. Note that when the number of phases is more than one, the cross-products are measured within each aggregation period, not across crop development stages.

2.1.2 Generalized additive model (GAM)

GAM is a regression model encompassing a linear combination of smoothers (Hastie and Tibshirani, 1987). This paper utilizes Penalized B-spline (P-spline) tensor products as smoothers since they are advantageous in modeling nonparametric interactions (cf. Lee et al., 2013). For instance, Tan and Zhang (2023) model the interplay between precipitation and vapor pressure deficits on yield losses using a B-spline tensor product and penalty function. Zou et al. (2023) implement a P-spline tensor product to design index insurance with GDD and CR. Despite the intuitive visualization, the model performance was not satisfactory in their case. Again, we use a B-spline tensor product as the type of smoother, S(WI1m,l,WI2m,l), under the whole growth cycle case. To construct a B-spline tensor product, we first conduct the de-Boor recursion (de Boor, 1978) to calculate the marginal B-spline basis Nj,du of each input weather covariate as follows:

(2)Nj,0(u)={ 1,if uj≤u<uj+10,otherwise and

(3)Nj,d=u−ujuj+d−ujNj,d−1(u)+uj+d+1−uuj+d+1−uj+1Nj+1,d−1(u),

where u represents the element in the weather index vector, the order d=1,2,3 when selecting cubic spline, {uj}, j=0,…,nk-1, is a uniform knot vector generated by the equidistant knots in the range of the weather index vector, and nk is the number of knots (cf. Tan and Zhang, 2023). Secondly, a P-spline tensor product is defined as the Kronecker product of B-spline marginal bases along two dimensions of weather indices with penalization. Being consistent with the variables and denotations in 2.1.1, a GAM that contains one smoother is Ym,l=S(WI1m,l,WI2m,l)+εm,l′.

2.1.3 Random forest (RF)

As introduced by Breiman (2001), the RF algorithm first draws a bootstrap sample B times from the dataset. This involves randomly selecting observations with replacements to create each bootstrap sample. These bootstrap samples are then used to grow individual trees in the RF ensemble T(WIb; Θb), b=1,…,B, where: WIb represent the b th group of input weather indices; Θb is the b th independent and identically distributed random vector constituted of corresponding split variables, cut points at each node, and terminal values of the tree; and B is the total number of trees (Hastie et al., 2009b). Finally, predictions are obtained by averaging the ensemble of unweighted trees, i.e. Ym,l=1B∑b=1BT(WIm,l,b; Θm,l,b). Through this process, the algorithm builds RF and identifies the best-split variable and cut points to divide nodes into two daughter ones, iteratively continuing until reaching the minimum node size. Prasad et al. (2021) found that the performance of RF is reliable and faster in predicting cotton yield. Schierhorn et al. (2021) also implement RF to explore the contribution of weather to winter wheat yields across various crop growth stages. Nevertheless, RF shows poor extrapolation capabilities when the training data are narrowly distributed (Montes et al., 2021).

2.1.4 Support vector machines (SVM)

While initially introduced for classification tasks (Cortes and Vapnik, 1995), SVM has gained popularity in regression applications (Gunn, 1998). By leveraging SVM, data points are projected into a higher-dimensional space using kernel techniques (Meyer, 2023). To better illustrate the application of SVM in weather-yield regression, we take the whole-cycle (non-phase-division) case for denotation simplicity. First, we reprocess the input weather indices WI1m,l and WI2m,l, by mapping them into a feature space (Aizerman, 1964) and the map is denoted as Φ(WI1m,l,WI2m,l). The SVM algorithm depends on the dot product of this map, which is defined as Ym,l=〈Φ(WI1m,l,WI2m,l),Φ(WI1m,lT,WI2m,lT)〉, where the superscript T denotes the transpose of the covariate vector. However, it could be computationally infeasible to derive the explicit forms of both the map function and the dot product, especially in the case of complex nonlinearity (Smola and Schölkopf, 2004). To solve this problem, “kernel tricks” are implemented to replace the dot product with a kernel function (Brereton and Lloyd, 2010), which is represented as K(WI,WIT)=〈Φ(WI),Φ(WIT)〉, where WI denotes the matrix of weather covariates [WI1m,l,WI2m,l]. Therefore, the optimization task can be accomplished with this kernel function and a commonly utilized one is the radial basis function (cf. Dang et al., 2021; Lischeid et al., 2022). For details on the optimization process, refer to Smola and Schölkopf (2004). Some studies conclude that SVM regression can achieve more accurate results in agricultural yield prediction (e.g. Brdar et al., 2011; Shah et al., 2018). Chen et al. (2016) also find that SVM performs dominantly in selecting climate variables that impact yield variability. SVM handles small training data sets effectively, but the estimation results are highly susceptible to hyper-parameter settings (Mountrakis et al., 2011).

2.1.5 Artificial neural networks (ANN)

ANNs are designed to simulate the functioning of neurons in the human brain within a computer system (Krogh, 2008). The network architecture typically consists of the number of layers (such as input, hidden, and output layers), the number of neurons within each layer, the activation function used in each layer, and the training algorithm employed (Benardos and Vosniakos, 2007). The construction component, namely the artificial neuron, is defined as Ym,l=g(∑h=1Hwm,l,hWIm,l,h+bias), where g(∙) is the activation function, h denotes the ordinal number of the input weather index, H is the total number of input variables, wh is the weight coefficient associated with the h th weather index WIh, and bias denotes the biases. Schmidt et al. (2022) employed the rectifier linear unit (ReLU) activation function, i.e. g(v)=max(v,0) (Sharma et al., 2017), to estimate the nonlinear weather-yield relationship. Moreover, the Adam optimizer (Kingma and Ba, 2014) was utilized in the training process. By leveraging ANN, intricate and nontrivial patterns in the data can be captured, allowing for a more comprehensive and accurate modeling of the weather-yield relationship (Schmidt et al., 2022; Chen et al., 2023). However, ANNs are commonly used in big data analysis and often more factors need to be taken into consideration when implementing ANNs in small sample sizes (e.g. Pasini, 2015; Shaikhina and Khovanova, 2017).

2.2 Simulation study

Despite the wide use of the above methods in weather-yield modeling, no model has been recognized as superior. It calls for methodological guidance through a thorough comparison between these methods, which is one of the contributions in this study. To illustrate their characteristics and potential flaws, we conduct a simulation study with artificial data generated by four nonlinear functions (the functional form is specified in the Appendix [cf. Lee et al., 2013]). We sampled uniformly distributed inputs xi,i = 1, …,8 with a sample size of 1,000 and evaluated the performance with the added white noises. The value of the function filled with the sampled observation is considered as a response variable. Note that for the hyper-parameter tuning process of SVM and ANN, identical training and test datasets were employed. The results are in Table A1 in the Appendix. The simulation study results are presented in Figure 1, which showcases the contour plots of true functions in blue and the estimated surfaces produced by the five methods in red fixing x3,x4,x5,x6,x7,x8 = 0.50. GAM demonstrates its robustness in approximating the true surface, while PR, RF, SVM, and ANN exhibit noticeable discrepancies consistent with the descriptions of methods above. First, PR shows the most prominent performance deficiency in predicting nonlinear relations. Second, analyzing the limited sample size, SVM outperforms ANN. Third, as shown by the vacant area around the point (x1 = 0.25, x2 = 0.50) RF is insufficient in extrapolation. These findings highlight the varying performance of the different methods in capturing the underlying relationships between covariates, which motivate the objective of comparing methods in this empirical study. All the models have problems in performing well at the boundaries due to lacking degrees of freedom.

This simulation study was intended to show weak points of the methods. In real-world applications, the function linking weather and yield are unknown. Therefore, we need to rely on other benchmarks or comparisons, keeping in mind the models’ weaknesses.

3. Data and study region

Our study area, the state of Saxony, is located in eastern Germany. According to the Saxon State Office for the Environment, Agriculture and Geology, a total of 12 agro-ecological zones are defined across the state based on soil properties and climate conditions. Figure 2 illustrates the distribution of these agro-ecological zones in polygon form. From north to south, soil is separated into three regions with relatively fertile soils in the central part. For details on the soil conditions, please refer to the map in the Appendix. According to the 2020 agricultural census in Germany, about 6,500 farms in Saxony cultivate around 894,100 hectares of arable land. The most important grains are wheat with around 182,700 hectares and barley with 115,200 hectares. The most important oil crop, winter rapeseed, is grown on 101,500 hectares (StLa Sachsen, 2020). From an email survey conducted in Saxony in 2022 by the Saxon State Ministry for Energy, Climate Protection, Environment and Agriculture, weather risks, e.g. hail, drought, frost, and ice days, are highly agriculturally relevant with drought ranked as the most common agricultural risk that resulted in yield losses. Nonetheless, the surveyed farmers were rather dissatisfied with the market’s risk management instruments. At the same time, they rated insurance as an attractive instrument, especially in view of the increase in weather risks caused by climate change. This real-world demand motivates our research to design a weather index insurance product for crop producers in Saxony.

Our empirical analysis incorporates yield, meteorological, and phenological data. Yield data were obtained from the annual financial statements of the Federal Ministry of Food and Agriculture through the Farm Accountancy Data Network, which is a microeconomic data source that relies on national surveys within Europe. Farm-level yield data provide advantages over county-level aggregations when estimating yield variability, as highlighted in previous studies (Popp et al., 2005; Schmidt et al., 2022). In our case study, we focus on winter barley [1] for two reasons. First, our data set is most complete for this crop, which are measured at the farm level and originate from 217 farms located in 157 municipalities across 12 agro-ecological zones in Saxony. The data spans from 2000 to 2015. For confidentiality reasons, the exact coordinates of farm locations are concealed, but information about the municipality a farm is located in is available. In Figure 2, these municipalities are shown as white points. Second, winter barley had a significant cultivation in Germany between 1989 and 2020 (Ellsäßer, 2022). Moreover, winter barley is a globally important cereal crop, widely used in animal feed, beer, and spirits production (Verstegen et al., 2014). Thus, there is potential demand for weather index insurance that protects winter barley producers in Saxony against weather perils. For a robustness check of our findings, we estimate the weather-yield relation for rapeseed, which is also an important economic crop in Saxony. The rapeseed yield data in our data set cover 202 farms located in 142 municipalities from 2000 to 2015 and the results are presented in the Appendix.

In this study we define yield deviations Ym,l=DYm,l−Y̅l as the response variable, where m represents the year, l denotes the farm, Y̅l is the average yield of farm l, DYm,l are detrended yield observations, and negative deviations indicate farm losses. We detrended the raw yield observations YOm,l using a linear time regression model YOm,l=α+δ(m−2000)+DYm,l (Shi and Jiang, 2016). The detailed descriptive statistics of YOm,l, DYm,l, and Ym,l are displayed in Table 1.

Daily weather data for Germany, including precipitation, maximum and minimum temperatures, are publicly available through the Deutscher Wetterdienst (DWD). Specifically, in the daily Hydrometeorologische Rasterdaten (HYRAS) data sets of the DWD that we use, the precipitation product for Germany is in a 1  × 1 km grid and the spatial resolution of the temperature raster is 5  × 5 km. To unify the spatial resolution of raw weather files, we use resampling techniques to obtain the 1  × 1 km gridded temperature data. To relate the weather data with municipalities, we calculate the average values in the polygons of municipalities and export the means.

In this study, we choose Growing Degree Days (GDD) and Cumulative Rainfall (CR) as weather indices. GDD and CR are calculated according to the following definitions:

(4)GDDm,l,s=∑t=sdsedsmax{0,(TMAm,l,s,t+TMIm,l,s,t)2−TB} and

(5)CRm,l,s=∑t=sdsedsRm,l,s,t,

where s∈{1,2,3,4} denotes the ordinal number of winter barley growth stage, t is the t-th day in a year, sd and ed denote the start and the end dates, respectively, of the accumulation period, whose unit is day of the year (DOY), and TMA, TMI, and R are the daily maximum temperature, minimum temperature, and rainfall, respectively. The baseline temperature TB of winter barley is defined as 0° Celsius (Manderscheid et al., 2009). Start dates of the accumulation periods are not dependent on the year m and we assume no temporal overlapping between two consecutive phases, i.e. sds+1=eds for s∈{1,2,3}. To calculate the weather indices GDD and CR for each farm, we average the gridded daily weather data collected from the DWD HYRAS data sets in each municipality instead of each farm because exact farm locations are unknown. However, the geographical basis risk is negligible due to the relatively small size of the municipality.

DWD also provides phenology reports in Saxony, but the spatial resolution of the data is not satisfactory. Although there are data from 57 phenology observation stations available, the number of available phenological stations varies each year with phenology timing. When we link each municipality with the nearest phenology station to determine the beginning and end dates, a discontinuity of data would be introduced. Thus, for each phase in each year, we define the medians of observations from available stations across Saxony and ignore spatial heterogeneity. Likewise, the temporal resolution of the phenology information is not expedient because it is not available on a weekly basis (Vallentin et al., 2022). To overcome this limitation, we follow the procedures in Schmitt et al. (2022) by deriving dates based on historical time distances between consecutive phenological phases and expert knowledge. Specifically, we first obtain the medians of all observed data as state-level phenology dates as described above. Then, we add a certain number of days to estimate the missing phenological information. For instance, the timing of “flowering” sd2 is defined by experts as the observation data of “heading” plus 11 days, while the occurrence of “fruit formation” sd3 is considered as 21 days after “heading”. Therefore, the temporal distance between “flowering” and “fruit formation” is as long as ten days across 16 years. The determination of phase length is highly dependent on expert experience. Thus, it can introduce additional errors due to individual cognition biases. The lack of appropriate phenology data is exactly the reason that motivates a data-driven growth phase determination.

To test the effectiveness of the data-driven phase division, we specify three distinct models:

1. A model without separation of growth phase (“non-division”, hereafter “ND”).

2. A model with exogenously determined length of growth stages (“expert-division”, hereafter “ED”).

3. A model with data-driven determination of growth phase length (“optimal-division”, hereafter “OD”).

The first two models serve as benchmarks for the third. Models (2) and (3) are referred to as “phase-division” models.

Following the methodology proposed by Schmitt et al. (2022), we divide the entire growth cycle of winter barley into four stages and the intervals of stage 1 (from “shooting” to “flowering”), stage 2 (from “flowering” to “fruit formation”), stage 3 (from “fruit formation” to “ripening”), and stage 4 (from “ripening” to “harvest”) are [sd1,ed1], [sd2,ed2], [sd3,ed3], and [sd4,ed4], respectively. Figure 3 shows DOYs for each growth phase from 2000 to 2015. Throughout the 16 years, the earliest occurrence of the “shooting” stage was on the 93rd DOY in 2014, while the latest occurrence of the “harvest” stage was on the 201st DOY in 2004. Therefore, we define the start and end division points as sd1=93 and ed4=201, respectively. Moreover, we denote the median DOY of internal division dates across 16 years as sd∼2=145, sd∼3=155, and sd∼4=177. sd1 and ed4 mark the beginning and end of the vegetation period in the non-division model. Note that these two timestamps are constant and identical in all models. In the expert-division model, the three internal division points are fixed at sd∼2, sd∼3, and sd∼4. With these specifications, GDD and CR time series are calculated for all ND and ED models. Table 2 displays the descriptive statistics of these weather indices. It indicates that the non-division weather indices vary greatly from expert-division ones regarding mean, standard deviation, minimal value, and maximal value. In addition, significant difference among stages in the statistics of expert-division weather indices is shown (see Table 2).

To conduct an optimal phase division, we parametrize sd2, sd3, and sd4 and enlarge the search space to k possible combinations, which amount to 24,804 combinations of growth phase durations in this study, given the constraints of minimal phase lengths, namely:

(6){sd2≥sd1+ml1−1sd3≥sd2+ml2−1sd4≥sd3+ml3−1sd4≤ed4−ml4+1,

where mls, s=1,2,3,4 denotes the minimal phase length of the s-th stage and ml1, ml2, ml3, and ml4 are 26,10,15, and 10, respectively (see Figure 3). For each method, the optimal combination of sd2, sd3, and sd4 is denoted as (sd¨2,sd¨3, sd¨4), respectively. This combination is selected by maximizing adjusted R2 or minimizing the Root Mean Squared Error (RMSE), i.e. RMSE=∑m=1n∑l=1217Ym,l−Yˆm,l2217n, where n is the length of time series in each farm in the training data set, m denotes the year, l represents the farm, Ym,l are the yield deviations, and Yˆm,l are the fitted values.

4. Model estimation and results

4.1 Parameter tuning

As described in Section 3, three types of models, namely non-division (ND), expert-division (ED), and optimal-division (OD) models, are estimated in this study. Among the five methods implemented, PR, GAM, and RF can be used without parameter tuning. All computations were performed using the R statistical software (version 4.3.1). For GAM, we use the gam() function using the P-spline and REML criterion in the “mgcv” package. RFs were estimated using the “randomForest” package (Breiman, 2001). The number of trees is set to 500 and the number of variables at each split is 2 (cf. Biau and Scornet, 2016). In contrast, the estimations of the two machine learning methods, SVM and ANN, contain an extra process of hyper-parameter selection.

As introduced in Section 3, for each type of model, we calculate accumulative measures GDD and CR and the input data are normalized into the range [0,1] for computational stability. With constant division points, the weather indices of ND and ED models are static and the hyper-parameter tuning procedures required by SVM and ANN are accomplished with a normal cross validation process.

However, for the optimal-division model, we obtain k combinations of both the division points and corresponding weather indices, which implies k times of cross validation and estimation for each method before we determine the optimal-division points in terms of adjusted R2 or RMSE as introduced in Section 3. For the methods that require no parameter tuning, namely PR, GAM, and RF, it is not an overly complex task, but for SVM and ANN, it would be too time consuming, since k=24,804 in our case. To balance the trade-off between costs and hyper-parameter accuracy, we first selected 20 division points combinations randomly from k possibilities. Second, for each selected division points combination, we selected hyper-parameters via cross-validation. We divided the dataset into a training set spanning 11 years and a test set of 5 years, which ensures a substantial test dataset size of 1,085 (the product of 5 years and 217 farms), that avoids misevaluation, e.g. spurious higher accuracy in machine learning implementation (Vabalas et al., 2019). Afterward, results of 20 sets of hyper-parameters were averaged for evaluating k models.

In the cross-validation process, we selected a radial basis kernel function for SVM based on the methodology described in Lischeid et al. (2022). The parameters “cost” and “gamma” were determined through grid search. Regarding ANN, we designed a fully connected ANN with two hidden layers and determined the hyper-parameters using grid search employing the ReLU activation function (Sharma et al., 2017) and the Adam optimizer (Kingma and Ba, 2014). The number of neurons were selected from 4, 8, and 16 in two-dimensional ND models. In the two phase-division models, the options were 8, 16, and 32. We also incorporated ridge regression as a penalty term corresponding to L₂ regularization to prevent overfitting. The results of the parameter tuning process in empirical study are provided in Table A2 in the Appendix.

4.2 Model results

After the estimations with ND, ED, and OD (24,804 rounds) models using each method, we obtained out-of-sample results (displayed in Table 3). Note that in the ND model estimated by ANN, we used the sigmoid activation function instead of ReLU because our parameter tuning tests indicate that the model performance is unstable using ReLU in a two-dimensional case. Moreover, the ND model was not estimated by RF since with two independent variables, the RF algorithm would randomly select at most one predictor from the two predictors as a subset and thus cannot further select the best-split variable recursively.

We compare the model performance and the division points based on the results. First, Table 3 shows that RMSE is around 10 dt/ha in most cases. This indicates a high basis risk from weather index insurance since 10 dt/ha is about 18% of the average detrended yield for winter barley (55.24 dt/ha). This reflects the general limitation of the prediction capability of weather indices for farm level yields. For the ND models, the performance differences of the four methods are negligible. In both phase-division models, the results differ more across methods, but they are still minor. Overall, PR has a relatively weak model fit, which is in line with earlier studies (e.g. Shah et al., 2018). In contrast, SVM ranks first among the five methods using separate growth phases. This confirms the finding of a review study of SVM in agricultural applications (Kok et al., 2021). GAM and RF are rather similar regarding model performance.

Table 3 further shows that phase-division models significantly outperform models that do not disaggregate the vegetation period into growth phases. The relative reduction in RMSE ranges from 14.0% to 35.2%, confirming the advantage of the phenology-based disaggregation in the weather index calculation of the estimation of the weather-yield relation that has been reported in previous studies, albeit with different study areas and crops (Zou et al., 2023; Dalhaus et al., 2018). Comparing the two phase-division models shows that model fit can be further improved by the proposed data-driven phase division, at least in combination with PR, GAM, and SVM. The performance gain form relaxing the internal three points is most significant for SVM. In contrast, for ANN this procedure comes along with a performance loss compared to the expert-based growth division model. Overall, the gains of optimal phase division compared to ED are only modest. However, our results demonstrate the potential of the data-driven phase-division as an alternative to expert-based phase division.

Apart from the model performance, it is interesting to compare the phase-division points derived by the different methods with the expert-based division points, i.e. sd∼2=145, sd∼3=155, and sd∼4=177 .We find that the combinations generated by PR [139, 149, 182], GAM [139, 149, 177], and ANN [151, 160, 177] deviate the least from the ED growth phases. Specifically, the sd¨4 of GAM and ANN is exact to the median sd∼4=177. Another remarkable finding is that under OD, sd¨4 of PR, GAM, RF, SVM, and ANN are 182, 177, 180, 189, and 177, respectively. This suggests that significant consistency exists in the third optimal point sd¨4, whereas it does not in the other two. However, since this study aims to search for the data-driven division points determined by the model selection, the optimal points chosen by each method do not need to be consistent with the expert-division ones.

Table 4 displays the lengths of the growth phases that result from the two approaches. The five statistical methods suggest different lengths, i.e. there is no unique OD. Only for the second phase do all methods suggest a similar length. This can be explained by the minimum length setting of 10 days, as introduced in Section 3. Due to data scarcity of phenology information from “flowering” to “fruit formation”, the period is fixed to be 10 days long across 16 years according to the expert definition. Thus, when defining the search space of (sd¨2,sd¨3,sd¨4), one constraint is that ml2=10 in Inequality (6). The results of the second phase in Table 4 indicate convergence to a short duration for winter barley for blooming and beginning the grain filling. For stage 1 (from “shooting” to “flowering”), the data driven determination results in a shorter length (except ANN) than the ED phase length, while for stage 3 (from “fruit formation” to “ripening”) the opposite is true. The regularity revealed by the data-driven results as opposed to expert knowledge can be explained by global warming induced by climate change. Increasing temperatures accelerate the flowering process (Distefano et al., 2018), slow the fruit formation period since it is more water intensive, and likely harm the soil’s water storage capacity (Werner et al., 2020).

4.3 Model interpretation

Machine learning methods are often considered as “black boxes” with a non-transparent link between input and output variables (e.g. Castelvecchi, 2016). Farmers, however, are interested in understanding the impact of weather events on crop yields, not only to verify indemnity payments from index insurance, but also to support managerial decisions, such as about irrigation and fertilization. For this purpose, we display partial dependence plots that predict the marginal effect of each explanatory variable on the response variable, while controlling for the other predictors through averaging their effects (Friedman, 2001). In the context of crop yield modeling, partial dependence plots have been employed by Jeong et al. (2016) and Beillouin et al. (2020).

Figures 4 and 5 show the partial dependence plot of the two weather indices (GDD and CR) for all growth stages. The y-axis on the left represents the yield deviations corresponding to the weather index value on the x-axis, excluding the effects of other predictors. The y-axis on the right side denotes the density of the weather index.

Figure 4 visualizes the effect of temperature on yield deviations. It is striking that the impact of GDD differs between growth stages. In phase 1 (from “shooting” to “flowering”), higher GDD₁ are associated with greater yields, while the opposite is true at the third stage (from “fruit formation” to “ripening”). These statistical findings can be explained from an agronomic perspective. In the early growth stage, high temperatures facilitate blooming and development after the vernalization process (Amasino and Michaels, 2010). In the second growth phase (from “flowering” to “fruit formation”), partial dependence plots of GDD₂ for PR and GAM are rather flat when the index value exceeds its mean. This is consistent with the results reported in Schlenker and Roberts (2009) where the yield growth curves of US corn, soybeans, and cotton all exhibit relatively flat shapes. In contrast, in our study, the plots of RF, SVM, and ANN show an inverse U-shape, meaning that deviations from an average temperature level in both directions cause yield losses. These losses, however, are rather small within the observed range of the GDD. The negative effect of higher temperature in phase 3 is consistent with the finding from the data-driven phase length in Section 4.2, i.e. decreasing soil water holding ability can accompany heat stress (Wahid et al., 2007; Werner et al., 2020). At the fourth stage (from “ripening” to “harvest”), the relatively flat partial dependence plots of ANN, RF, and SVM imply that impacts of GDD on yield deviations are modest when winter barley matures. The irregular and U-shaped plots produced by GAM and PR appear implausible.

Figure 5 depicts the partial dependence of winter barley yield on cumulative rainfall in the respective growth phases. For the first growth stage, all plots exhibit a negative slope as long as CR₁ is below its mean value. This confirms the well-known fact that insufficient precipitation in the early plant growth phase causes yield losses (e.g. Stephens and Lyons, 1998). The curves of ANN, RF, and SVM are rather flat around the mean values in stages 2, 3, and 4, indicating robustness of winter barley yields against drought as well as excessive rainfall. This “crop yield plateau effect” has also been reported in Turvey et al. (2021). One explanation is that the length of the second growth phase is relatively short (10 days) compared with the first one (52 days). One should also recall that these plots do not account for interactions among growth stages. For example, due to the water storage capacity of soil, rainfall in the first period can be carried over to subsequent periods. Using soil moisture instead of cumulative rainfall might capture these temporal spillover effects. Moreover, we cannot rule out that some farms use irrigation to compensate for rainfall deficits.

Overall, the partial dependence plots show similar patterns across the five statistical methods, indicating the robustness of the results – at least in the neighborhood of the indices’ mean values. However, some abnormal deviations occur for PR and GAM in some growth stages (e.g. GDD₄ and CR₂) if the indices take extreme values. PR behaves wildly because few observations on the boundaries have very little effect on the Ordinary Least Squares, but the shape of the function is determined as quadratic polynomial. Specifically, the curve of GAM shows more variation than the other four methods. In Alimadad and Salibian-Barrera (2011), GAM can be very sensitive to outliers, in spite of their small proportion in the observations. The instability of partial dependence plots can be considered as a reflection of the limitation of these methods. Moreover, some methods result in dependence plots that are almost constant in the tails of the distribution, i.e. RF for GDD₁ < 400 in Figure 4. The RF curve stays constant if there is a lack of data because the fitted values are obtained from average trees that provide constant values at the leaves. These extrapolation failures of tree-based methods, such as RF, together with other data-driven methods, such as SVM and ANN, have been proved in previous studies (cf. Hettiarachchi et al., 2005; Hengl et al., 2018).

5. Conclusions

This research was motivated by the observation that the performance of crop yield models can be improved by a decomposition of the entire vegetation period into distinct growth phase in which weather events may have a different effect on crop yields. This, in turn, facilitates the design of effective weather index insurance. Unfortunately, the required phenological information for the determination of the beginning and end of each growth stage is not readily available. This study proposes a novel data-driven growth phase division that does not hinge on expert-based information.

Our study offers several findings relevant for crop yield modelers, weather insurance developers, and farmers. First, we confirm the expectation that a breakdown of the vegetation period and the disaggregation of weather indices results in a better model fit. Second, the endogenization and optimization of division points leads only to a slight improvement compared with an expert-based determination. This arouses a discussion on the cost-benefit ratio of the optimal division procedure. Regarding the costs of the suggested data-driven determination of growth phases, one has to distinguish between computational costs, which are rather low, and development costs, which may be considerable. These costs have to be compared with the labor costs of an in situ data collection of phenological information. The benefits, on the other hand, depend on the availability of an expert-based phase division. If insurance providers have access to phenological information, there is only little incentive to develop a model with endogenous splitting dates for the growth phases, because the gain in model performance is minor. However, if expert knowledge is not available, which is often the case in developing countries, the data-driven phase division constitutes a feasible alternative. It is also noteworthy that the optimal division points generated by the various statistical methods differ, i.e. there is no unique optimal division. Nonetheless, model performance is rather robust against changes in the timing and length of the phases. Third, regarding the statistical approach, we find that the parametric PR model is outperformed by nonparametric and non-linear ones, i.e. either spline-based GAM or machine learning models (RF, ANN, SVM). In our application, SVM show the best model fit though they have primarily been designed for classifications tasks. Despite the initial constraints of our study, which focused on a small geographical area and a limited time interval, our findings indicate that the application of smooth, flexible methods, namely GAM, RF, SVM, and ANN, yielded comparable results in terms of R² and RMSE. Notably, we acknowledge that we did not explicitly test for statistical significance. Although ANN posed challenges in terms of training and tuning due to its computational intensity, the similarity in partial dependence plots among RF, SVM, and ANN suggests a comparable performance across these models. It would be interesting to examine whether these findings can be generalized to other crops and other climatic regions.

Despite of the aforementioned improvements, one has to acknowledge that the basis risk of a weather insurance product derived from these models is still high. Even in the best case, average prediction errors amount to more than 10% of the average barley yield and undermine the hedging effectiveness of index insurance. In fact, the uptake of index-based weather insurance in the study region is modest and higher subsidies from the government would be necessary to increase the adoption of this kind of index insurance product in reality (Odening et al., 2022). However, there are several opportunities for optimization to further reduce basis risk. First, the inclusion of additional or more sophisticated weather variables (e.g. rainfall deficit, dry spells, soil moisture, excessive rainfall, and heat days) will most likely enhance model performance (e.g. Bucheli et al., 2021). Second, satellite data have received considerable research attention to reduce the geographical basis risk (e.g. Möllmann et al., 2019; Vroege et al., 2021). Third, as pointed out by Bucheli et al. (2021), considerable risk reduction can be achieved if the underlying drought index is tailored individually for each farm. Finally, the parameter tuning process of SVM and ANN could be further refined.

Figures

Figure 1

Contour plots of fˆ(x1,x2,x3=x4=x5=x6=x7=x8=0.5) estimated by Polynomial Regression (PR), Generalized Additive Model (GAM), Random Forest (RF), Support Vector Machines (SVM), and Artificial Neural Networks (ANN), respectively

Figure 2

Agro-ecological zone classification in the state of Saxony

Figure 3

Day of the year (DOY) of phenological phenomenon timing

Figure 4

Partial dependence plots of growing degree days (GDD) in four growth stages

Figure 5

Partial dependence plots of cumulative rainfall (CR) in each growth stage

Table 1

Descriptive statistics of yield data

	Yield (dt/ha)	Detrended yield (dt/ha)	Yield deviations (dt/ha)
Mean	61.76	55.24	0
Standard deviation	17.03	12.74	12.74
Minimuma	0	−4.37	−59.60
Maximum	122.88	122.63	67.39

Note(s): ^aGiven that we use farm yield data and the fact that farms in East Germany encountered severe droughts in several years within the observation period, a zero yield observation is not unlikely. It simply means that the yield was so low that it was not profitable to incur harvesting costs

Source(s): Table created by authors

Table 2

Descriptive statistics of time series of growing degree days (GDD) (in degree Celsius) and cumulative rainfall (CR) (in mm) for non-division (ND) and expert-division (ED) models

	ND-GDD	ED-GDD1	ED-GDD2	ED-GDD3	ED-GDD4	ND-CR	ED-CR1	ED-CR2	ED-CR3	ED-CR4
Mean	1509.06	575.88	155.44	376.70	446.70	234.86	83.86	34.92	53.01	68.74
Standard deviation	105.11	65.44	28.35	39.85	42.73	66.75	28.17	41.31	26.63	35.68
Minimum	1120.83	371.50	78.79	255.64	293.29	71.86	24.43	0.00	0.69	2.21
Maximum	1746.02	749.97	217.01	459.30	557.07	547.86	187.32	250.04	159.80	195.98

Source(s): Table created by authors

Table 3

(sd¨2,sd¨3,sd¨4), RMSE, and adjusted R2 of ND, ED, and OD models

Method	ND	RMSE	Adj. R2	ED	RMSE	Adj. R2	OD	RMSE	Adj. R2
PR	–	12.00	0.11	(145, 155, 177)	10.76	0.28	(139, 149, 182)	10.32	0.34
GAM		11.75	0.15		10.00	0.38	(139, 149, 177)	9.61	0.43
RF		–	–		9.80	0.41	(121, 132, 180)	9.74	0.41
SVM		11.77	0.14		9.01	0.50	(132, 141, 189)	7.62	0.64
ANN		11.75	0.15		9.11	0.49	(151, 160, 177)	9.51	0.44

Note(s): ND, ED, and OD are abbreviations for non-division, expert-division, and optimal-division. PR, GAM, RF, SVM, and ANN are abbreviations for Polynomial Regression, Generalized Additive Model, Random Forest, Support Vector Machines, and Artificial Neural Networks, respectively

Source(s): Table created by authors

Table 4

ED and OD phase lengths (in days)

Division model	Stage 1	Stage 2	Stage 3	Stage 4
ED	52	10	22	24
OD-PR	47	11	34	20
OD-GAM	47	11	29	25
OD-RF	29	12	49	22
OD-SVM	40	10	49	13
OD-ANN	59	10	18	25

Note(s): ED and OD are abbreviations for expert-division and optimal-division. PR, GAM, RF, SVM, and ANN are abbreviations for Polynomial Regression, Generalized Additive Model, Random Forest, Support Vector Machines, and Artificial Neural Networks, respectively

Source(s): Table created by authors

Notes

1. Winter wheat is the most important crop in Saxony and we have considered analyzing this crop, but unfortunately in our data set winter wheat is merged with summer wheat. Estimating a yield model for this aggregate would be inappropriate due to the different vegetation periods of winter and summer wheat.

Appendix

The supplementary material for this article can be found online.