2.1. Overview of current crop yield forecast methods

There are two types of crop forecasts: survey-based forecasts and regression-based forecasts. Survey-based forecasts tend to be more accurate, especially when the harvest date is approaching, usually available shortly before or around harvest time, but they are also more expensive and labor intensive. Regression-based forecasts are more cost effective and can be available largely ahead of harvest; however, their accuracy may be compromised.

Survey-based forecasts that are used by the U.S. Department of Agriculture, National Agricultural Statistics Service (USDA-NASS) are made by conducting annually an agricultural yield survey (AYS) and an objective yield survey (OYS), the details of which can be found in USDA-NASS (2012). In the AYS, farmers are asked to self-report their anticipated yields, which may become the actual yields if harvest has begun. In the OYS, NASS sends technical personnel to the field to take objective measurements and counts of the plants. Both AYS and OYS are conducted monthly from May to November, but soybean yield data are collected and soybean yield forecasts are published from August to November. The final forecast is released in January of the next year. The typical cycle of soybean production in the major producing states in the United States is as follows: planting is in May and June, flowering is in July (which is its moisture/temperature-sensitive stage), filling is in August, maturation is in September, and harvesting is between October and November.

The second type of crop forecast is the regression-based forecast. This type of forecast is used mostly by private agencies and occasionally as a supplementary forecast by public agencies. For example, the World Agricultural Outlook Board (WAOB) releases the World Agricultural Supply and Demand Estimates regression-based forecasts, which use trend analysis and crop weather regression models. Unlike the forecasts released by NASS at the end of the year, the WAOB releases early forecasts throughout the growing season, from May to August (Irwin, Sanders, and Good, Reference Irwin, Sanders and Good2014). The comparison between NASS and WAOB yield forecasts and an evaluation of WAOB forecast accuracy can be found in Irwin, Good, and Sanders (Reference Irwin, Good and Sanders2015). The crop weather model (also known as the modified Thompson model) utilizes a year trend variable, monthly weather variables, and an indicator if the crop is planted late. The crop condition model utilizes a year trend variable, the proportion of the crop planted after a certain date (e.g., May 30 for soybeans; Irwin, Good, and Tannura, Reference Irwin, Good and Tannura2009), and the proportion of the crop rated as good or excellent by USDA (Crop Progress Report). The model we propose in this study is based on the crop weather model but also adds NDVI variables. According to the literature, the modified Thompson model produces a good fit but performs poorly when events (such as insects and diseases) that cannot be captured by a weather variable negatively affect crop yields. We hypothesize that using NDVI can also monitor for insects and diseases because NDVI is a direct indicator of the greenness/health of the vegetation, with the additional benefit that NDVI data are immediately available at a low cost compared with the methods that rate crop conditions. Because regression-based forecasts typically rely on aggregate-level information, such as climatological variables at the county or regional level, a limitation of the regression-based forecasts is their inability to incorporate farm-level characteristics such as managerial skills or soil characteristics. However, regression-based forecasts can become useful when farm-level data are lacking, which is prevalent in many cases, especially in yield forecasting in developing countries.

2.2. Crop yield forecasting using remote sensing

There have been numerous studies documenting the correlation between NDVI and crop yield forecasts, at the national (Maselli and Rembold, Reference Maselli and Rembold2002), regional, county (Bolton and Friedl, Reference Bolton and Friedl2013), and field level (Ferencz et al., Reference Ferencz, Bognar, Lichtenberger, Hamar, Tarcsai, Timar and Szekely2004). Tucker (Reference Tucker1979) determined that a time-integrated NDVI is largely correlated with crop yields when the vegetation is at the maximum level of greenness. Some studies focus on intra-annual variability showing how the correlation between the vegetation index and crop yields varies by the crop cycle and planting date (Basnyat et al., Reference Basnyat, McConkey, Meinert, Gatkze and Noble2004). These studies suggest choosing NDVI data over a specific period for each type of crop in order to produce better forecasts. The weekly availability of NDVI data makes this crop-specific specification achievable. Lv (Reference Lv2014) suggests using earlier May NDVI and the change in NDVI over the crop planting and harvesting season for the most accurate yield forecasting. Johnson (Reference Johnson2014) finds that crop yields are highly correlated with NDVI and daytime land surface temperature. The author conducts a regression of crop yields on NDVI for every week of the growing season and finds that the week in which the correlation is at its peak is at the beginning of August.

In addition to NDVI derived from the National Aeronautics and Space Administration’s (NASA) Earth Observing System (EOS) Moderate Resolution Imaging Spectoradiometer, called eMODIS, other indexes and images have been used. For example, Doraiswamy and Cook (Reference Doraiswamy and Cook1995) is one of the earliest studies that used Advanced Very High Resolution Radiometer (AVHRR) imagery. AVHRR data are coarser, whereas eMODIS data are finer; AVHRR data are available for an extended period, whereas eMODIS data are only available after 2000. Later, Ferencz et al. (Reference Ferencz, Bognar, Lichtenberger, Hamar, Tarcsai, Timar and Szekely2004) also used AVHRR and a vegetation index called general yield unified reference index. Bolton and Friedl (Reference Bolton and Friedl2013) suggest to incorporate crop phenology and use a combination of the EVI2 (two-end enhanced vegetation index), NDVI, and normalized differenced water index (NDWI) for crop yield forecasting. They distinguish between semiarid and non-semiarid areas. They find that vegetation indexes are the best type of indexes for predicting in non-semiarid areas, whereas the NDWI is the best index for prediction in semiarid areas, because the water index is sensitive to irrigation in these semiarid areas.

Instead of using traditional statistical models, Bose et al. (Reference Bose, Kasabov, Bruzzone and Hartono2016) utilize spiking neural networks from machine learning to analyze a remote sensing spatiotemporal relationship. Their work focuses on finding the optimum number of variables (or “features” in machine learning) to be included in regression analysis using machine learning techniques. They find that this type of prediction can be made 6 weeks before harvest with an average accuracy of 95.64%. They find that the year 2002 had the largest forecast error because of the 2002 drought. Adrian (Reference Adrian2012) applies the Bayesian hierarchical model. This model is suitable for modeling data with clusters. It produces unique estimates for each state while requiring the estimates from each state to also follow a prior distribution. Johnson et al. (Reference Johnson, Hsieh, Cannon, Davidson and Bédard2016) focus on comparing forecast performance using linear versus nonlinear machine learning techniques and find that nonlinear models are not necessarily advantageous compared with linear models. Li et al. (Reference Li, Liang, Wang and Qin2007) find that neural network techniques improve corn predictions compared with multivariate analysis. Kaul, Hill, and Walthall (Reference Kaul, Hill and Walthall2005) find that a nonlinear model only outperforms the linear model for barley. Mkhabela et al. (Reference Mkhabela, Bullock, Raj, Wang and Yang2011) categorize the Census Agricultural Regions (CARs) into three distinct agroclimatic zones; however, even within CARs, there might be multiple soil types. Bolton and Friedl (Reference Bolton and Friedl2013) emphasize the importance of delineating the boundary between farmland and nonfarmland, such as grassland and forests, because nonfarmland may contaminate the NDVI–crop yield relationship. Delineation can be done by using a land cover map such as Landsat Thematic Mapper (TM) data (Bolton and Friedl, Reference Bolton and Friedl2013). Another method of delineation is to identify single pixels as agricultural or nonagricultural vegetation using statistical correction analysis (Maselli and Rembold, Reference Maselli and Rembold2002). Among those studies, there are soybean forecasts in the United States using remote sensing (Lobell and Asner, Reference Lobell and Asner2003; Prasad et al., Reference Prasad, Chai, Singh and Kafatos2006). Chang et al. (Reference Chang, Hansen, Pittman, Carroll and DiMiceli2007) focus on using NDVI to map corn and soybean farmland.

Fieuzal, Sicre, and Baup (Reference Fieuzal, Sicre and Baup2017) make corn yield forecasts using both a real-time approach and a diagnostic approach. The real-time approach updates the estimates dynamically after the newest image is acquired, whereas the diagnostic approach utilizes all the image data throughout the season. The authors find the two best estimates perform comparably. Burke and Lobell (Reference Burke and Lobell2017) regress the agreement between satellite-based yields and field-reported yields as a function of farm size and find the vegetation index can most accurately predict crop yield when the field size is large.

All of the abovementioned studies employ a global model to produce the regression results that fit all observations, with the major difference among the studies being the specific model they use. To the best of our knowledge, this study is the first one to employ models that produce site-specific regression results, allowing heterogeneous responses of soybean yields across counties. This is also the first study to our knowledge that applies the FFT model to examine the yield-NDVI relationship.

3.1. Data

We use data for 797 counties from 10 major soybean-producing states in the United States from 2000 to 2016. According to NASS, the soybean production from these 10 states accounted for 78.5% (in 2016) and 79.8% (2000–2016 average) of the total soybean production in the United States (see Table 1 for soybean production and yield by state). Mkhabela et al. (Reference Mkhabela, Bullock, Raj, Wang and Yang2011) state that if a crop is not the dominant crop in the region, NDVI would give a poor prediction of crop yield because it cannot distinguish between different crops. The soybean yield data are obtained from the USDA-NASS QuickStats (https://quickstats.nass.usda.gov [accessed December 1, 2017]). This database provides official published aggregate statistics on U.S. soybean yields and the value of soybean production. Soybean yield is measured in bushels per acre. The NDVI data we use are from eMODIS onboard NASA’s EOS Terra satellite. Landsat TM and eMODIS are two mainstream imagery sources. Though Landsat TM has a better spatial resolution (30 m) than eMODIS (250 m), the latter provides a better temporal resolution (daily) than the former (16-day cycle). For monitoring purposes, we chose the eMODIS data. The eMODIS instrument onboard the Terra satellite achieves global coverage on a daily basis and provides 7-day composited data sets for its suite of products. Each data set provides NDVI information in GeoTIFF format that contains the reflective indices captured by Terra satellite at the resolution of 250 m from 2000 onward. Ag-Analytics converts the 250-m-resolution raw images to county-level NDVI. Ag-Analytics is an open-source, open-access database that provides data on agricultural finance, environmental finance, insurance, and risks (Woodard, Reference Woodard2016). We calculate county-level monthly NDVI values by taking a monthly average of the weekly NDVI values provided by Ag-Analytics. Climatological data are obtained from PRISM (parameter-elevation regressions on independent slopes model) Climate Data from Oregon State University and Ag-Analytics. We include two weather variables: maximum temperature over a month and average monthly precipitation. County boundary shapefiles are obtained from the U.S. Census Bureau. We obtain a sample of 12,027 county-year observations for the FFT analysis.

Table 1. Soybean production and yield in 10 major producing states

^a Soybean production is measured in 1,000 bushels. Soybean yield is measured in bushels/acre.

3.2. Flexible Fourier transform model

When estimating crop yield response to input variables, traditional models use regional and temporal dummies to capture spatial and intertemporal heterogeneity. Adding dummy variables can only capture the difference in the value of the dependent variable across locations and time; it does not take into account how the relationship varies according to site-specific and time-specific characteristics. Another type of model uses a quadratic functional form to estimate the relationship between crop yield and weather variables, assuming that crop yield is nonlinearly related to the weather variable. However, these models may suffer from model misspecification, especially if there is a threshold effect, driven by environmental risks such as drought and flooding (Cooper, Nam Tran, and Wallander, Reference Cooper, Nam Tran and Wallander2017).

Gallant (Reference Gallant1984) first proposed flexible Fourier functional transform to generate unbiased production function approximation and proved its mathematical validity. Cooper, Nam Tran, and Wallander (Reference Cooper, Nam Tran and Wallander2017) applied an FFT function to estimate the relationship between crop yield and temperature. We follow the approach and modeling in Cooper, Nam Tran, and Wallander (Reference Cooper, Nam Tran and Wallander2017) for the flexible Fourier function, which can be presented as follows:

$${\rm{Soybean\,yield}} = {\rm{}}{\beta _0} + \mathop \sum \limits_{m = April}^{August} ({\beta _{1m}}MaxTem{p_m} + {\beta _{2m}}MaxTempSquar{e_m}) + \mathop \sum \limits_{m = April}^{August} ({\beta _{3m}}Precipitatio{n_m} + \,{\beta _{4m}}PrecipitationSquar{e_m}) + \mathop \sum \limits_{m = April}^{September} ({\beta _{5m}}NDV{I_m}) + {\delta _0}TimeTrend + \mathop \sum \limits_{s = 1}^9 {\delta _s}StateDumm{y_s}\, + 2\mathop \sum \limits_{\alpha = 1}^A \mathop \sum \limits_{j = 1}^J \left\{ {{v_{j\alpha }}{\rm{cos}}\left[ {jk_\alpha ^{'}s\left( {NDVI} \right)\left] \ { - \ {w_{j\alpha }}{\rm{sin}}} \right[jk_\alpha ^{'}s\left( {NDVI} \right)} \right]} \right\} + error$$ (1)

In this model, the dependent variable is soybean yield in a county for a given year. β ₀ is the constant term. MaxTemp_m , Precipitation_m , and NDVI_m are the maximum temperature, the average precipitation, and the average NDVI in month m, respectively. We include the weather variables from April to August, following the standard specification in the literature (Cooper, Nam Tran, and Wallander, Reference Cooper, Nam Tran and Wallander2017). We include NDVI variables through September, following the remote sensing literature (Li et al., Reference Li, Liang, Wang and Qin2007). The advantage of the FFT function is that it not only allows for model flexibility but also incorporates multivariate estimation, which is difficult to achieve through other nonparametric models such as kernel regression.

PrecipitationSquare_m and MaxTempSquare_m are the squared terms of MaxTemp_m and Precipitation_m . TimeTrend equals the year minus 1999. StateDummy_s is the state dummy variable. NDVI is a vector with each element being NDVI_m . s(NDVI) is the scaled version of NDVI such that each element of s(NDVI) is in the range of [0, 2π]. In our case, only NDVI variables are transformed.

The ${\beta _0} + \sum\nolimits_{m = April}^{August} {({\beta _{1m}}MaxTem{p_m} + {\beta _{2m}}MaxTempSquar{e_m})} + \sum\nolimits_{m = April}^{August} {({\beta _{3m}}Precipitatio{n_m} + {\beta _{4m}}PrecipitationSquar{e_m})} + \sum\nolimits_{m = April}^{September} {({\beta _{5m}}NDV{I_m})}$ terms represent the quadratic regression part. β _1m, β _2m, β _3m, β _4m, and β _5m are parameters to be estimated. The $2\sum\nolimits_{\alpha = 1}^A {} \sum\nolimits_{j = 1}^J {\{ {{v_{j\alpha}}\cos [ {jk_\alpha ^{'}s(NDVI)} ] - {w_{j\alpha }}\sin [ {jk_\alpha ^{'}s(NDVI)} ]} \}} $ term models the functional flexibility using FFT. Similar to the Taylor expansion, which uses a series of polynomial terms to approximate the true function, the Fourier function uses a series of trigonometric terms to approximate the true function. The Fourier functional form is believed to be the only known functional form that satisfies the Sobolev condition, meaning that the difference between the approximated function and the true function approaches zero as the sample size becomes arbitrarily large. For a proof that the Fourier function satisfies the Sobolev condition, refer to Gallant (Reference Gallant1994). In the model, k_α (α = 1, 2, …, A) is the elementary multi-index vector, whose dimension equals the dimension of x_FFT, whereas A is the total number of elementary multi-indexes. The vector k_α can be obtained in the following way: first, exhaust the list of k_α, such that k_α has only integer elements and the sum of the absolute value of each element in k_α is no greater than K, where K is predetermined; second, delete any k_α whose first nonzero element is negative; and third, delete any k_α whose elements have a common integer divisor. Monahan (Reference Monahan1981) introduced a Fortran code to produce the set of elementary multi-index vectors. Also in the model, J is the order of the Fourier transformation, whereas v_jα and w_jα are parameters to be estimated. We use the following parametrization: K = 2, J = 2, which are chosen such that the rule of thumb—the number of variables after transformation is roughly the square root of the number of observations (Fenton and Gallant, Reference Fenton and Gallant1996)—is satisfied. Because there are 12,027 observations in the data we use, we include a total of 120 variables after the adding the transformed NDVI variables.

The model degenerates to the traditional OLS model when v_jα = 0 and w_jα = 0. In the following discussion, the OLS model refers to equation (1), with v_jα = 0 and w_jα = 0 imposed. By testing the statistical significance of variable v_jα and w_jα, we can decide whether the traditional quadratic model should be rejected in favor of the more flexible FFT model.

A review of the relevant literature reveals that the FFT model has been used/tested by scholars in different studies, fields, and situations. Chang et al. (Reference Chang, Kim, Miller, Park and Park2016) used the FFT to model the nonlinear effect of temperature on electricity demand. Becker, Enders, and Lee (Reference Becker, Enders and Lee2006) proposed a unit root test with a Fourier functional transform. Enders and Li (Reference Enders and Li2015) approximated structural breaks in U.S. GDP trends using Fourier forms. Jones and Enders (Reference Jones, Enders, Ma and Wohar2014) provided a summary on using Fourier forms to model structural breaks.

3.3. Prediction and forecast

We compare the prediction performance of the FFT model versus the OLS model. We conduct out-of-sample predictions and evaluate prediction performance by comparing prediction errors measured by the root-mean-square error (RMSE) and the mean absolute error (MAE), between FFT and OLS, for three schemes: time-series prediction, cross-sectional prediction, and panel prediction. RMSE and MAE are defined as follows:

$${\rm{RMSE}} = \sqrt {{1 \over {\rm{N}}}\sum\limits_{i = 1}^{\rm{N}} {{{({y_i} - {{\hat y}_i})}^2}} } .$$ (2)

$${\rm{MAE}} = {1 \over {\rm{N}}}\sum\limits_{i = 1}^{\rm{N}} {|{y_i} - {{\hat y}_i}|} .$$ (3)

Both RMSE and MAE are commonly used measures to evaluate prediction performance. They measure the difference between true and fitted values for soybean yield. The unit for both RMSE and MAE is bushels per acre. In a time-series prediction, we first select a year for prediction, then we use observations from all other years to generate the model, and after that we predict the soybean yield for the selected year using the fitted model, weather data, and NDVI data from the selected year. In cross-sectional prediction, similarly, we select a state for prediction, then we use observations from all other states to generate the model, and after that we predict the soybean yield for the selected state using the fitted model, weather data, and NDVI data from the selected state; in panel prediction, similarly, we make the prediction for a selected year and state. Though commonly used, a shortcoming of using RMSE or MAE to measure prediction performance is that we do not know whether the predicted yield overestimates or underestimates the final actual yield.

We make predictions and forecasts using the regression results from the models. In this study, prediction refers to cases where we may use data afterward to predict for a specific time; forecast refers to cases where we only use data up to a certain year to make predictions for that year.