Probabilistic forecasting of remotely sensed cropland vegetation health and its relevance for food security

In a world where climate change, population growth, and global diseases threaten economic access to food, policies and contingency plans can strongly bene ﬁ t from reliable forecasts of agricultural vegetation health. To inform decisions,itisalsocrucialtoquantifytheforecastinguncertaintyandproveitsrelevanceforfoodsecurity.Yet,inprevious studies both these aspects havebeen largely overlooked. This paper developsa methodology to anticipatethe agricultural Vegetation Health Index(VHI) while making the underlyingprediction uncertaintyexplicit.To achievethis aim, a probabilistic machine learning framework modelling weather and climate determinants is introduced and implemented through Quantile Random Forests. In a second step, a statistical link between VHI forecasts and monthly food price variations is established. As a pilot implementation, the framework is applied to nine countries of South-EastAsia(SEA)withconsiderationofnationalmonthlyriceprices.Modelbenchmarksshowsatisfactoryaccuracymet-rics, suggesting that the probabilistic VHI predictions can provide decision-makers with reliable information about future cropland health and its impact on food price variation weeks or even months ahead, albeit with increasing uncertainty as the forecasting horizon grows. These results - ultimately allowing to anticipate the impact of weather shocks on household food expenditure - contribute to advancing the multidisciplinary literature linking vegetation health, probabilistic forecasting models, and food security policy.


Introduction
Since the introduction of agriculture, 7000 to 10,000 years ago, the link between weather and crop yield has been a constant reason for concern due to the potential impact of yield variability on food security (Trnka et al., 2016;Black and Thompson, 1978;Frieler et al., 2017;McKeown et al., 2006). Whilst technological change, such as irrigation, new seeds varieties, fertilisers, greenhouses, and land management techniques have Science of the Total Environment 838 (2022) 156157 dramatically mitigated the direct impact of weather on agriculture (Rasmussen, 1962;Parayil, 1992;Sullivan, 1984), weather shocks remain a major reason for concern for stakeholders in the agricultural value chain (Jägermeyr et al., 2021;Wang et al., 2021). This is particularly true in lower income regions of the world, where technology adoption and agricultural knowledge accumulation have lagged behind substantially compared to industrialised countries (Lybbert and Sumner, 2012). Besides weather variability, several studies have highlighted how anthropogenic climate change is and will increasingly be exerting a substantial, mostly negative impact on crops growth and health (Anderson et al., 2020;Mbow, 2020). Climate change is expected to negatively affect global crop yields, and ultimately food supply and accessibility, especially in the most vulnerable regions with fragile agricultural systems (Mora et al., 2015;Porter, 2014). Under these changing conditions and negative outlook, it has become of utmost importance to monitor and provide reliable predictions of vegetation health to better allocate scarce resources, design early warning systems, and ensure food security.
Remote sensing techniques enable collecting large amounts of granular information over vast areas and at different spatio-temporal resolution scales. By characterising natural features on the ground and monitoring their changes over time, remote sensing applications can ultimately support food-related policy-making at different levels (Sishodia et al., 2020;Vroege et al., 2021). The scientific community has made extensive use of satellite imagery for mapping and monitoring changes in land cover and estimating geophysical and biophysical characteristics of the soil (Shanmugapriya et al., 2019;Weiss et al., 2020), as well as to develop and validate vegetation health indicators (Kogan, 2002;Kogan et al., 2004;Kogan, 2019;Kogan et al., 2005;Xue and Su, 2017). The Vegetation Health Index (VHI) (Kogan, 1997;Kogan, 1987) is one such indicators and it has largely been used to monitor crop vegetation over large areas and predict crop yield (Kogan, 1990;Orlovsky et al., 2010;Rahman et al., 2009;Zuhro et al., 2020). In parallel, learning-based statistical algorithms have opened up new frontiers and enabled the development of tools for analysing satellite imagery, providing better and more nuanced insights thanks to their ability to find patterns underlying the complex nonlinear relations that characterise environmental variables. With regards to forecasting, Probabilistic Machine Learning (ML) (Duan, 2020;Ghahramani, 2015;Gneiting and Katzfuss, 2014;Palmer, 2012) has highlighted the importance of accompanying any statement about the future with a quantification of its uncertainty to give a complete picture of the possible scenarios and better inform decision-making at all levels.
In this context, previous research has applied a variety of statistical methods to anticipate vegetation health. Yet, virtually all studies neither analyse the prediction uncertainty nor seek to investigate the relevance of forecasts for food security. For instance, Nay et al. (2018) developed a gradient-boosted machine learning approach to predict the Enhanced Vegetation Index (EVI) based on MODIS satellite data. A similar strategy is followed by Perera et al. (2020), who expand the approach to a larger set of vegetation indexes and ML models for the case of the MENA region. Lees et al. (2022) test Long Short Term Memory Deep Learning architectures for Vegetation Condition Index (VCI) forecasting in Kenya. The only notable exception that explicitly models vegetation health forecasts uncertainty is the recent work by Salakpi et al. (2021), who apply a Bayesian Auto-regressive Distributed Lags (BARDL) model in the context of arid territories in Kenya to forecast the VCI up to 10 weeks ahead while providing a probabilistic prediction.
Irrespective of this evidence, no previous multi-country, large-scale study has sought to analyse the probability distribution of vegetation health forecasts. Yet, to inform policy it is crucial to quantify uncertainty and demonstrate the relevance of such vegetation health forecasts for food security. To fill those gaps, this paper introduces a probabilistic framework seeking to anticipate remotely sensed agricultural vegetation health based on previous weather conditions while making the underlying uncertainty explicit. This represents an important element of novelty compared to studies based on point estimate prediction of vegetation health. The approach is implemented and validated in nine countries of South-East Asia (namely, Cambodia, Indonesia, Laos, Malaysia, Myanmar, Philippines, Thailand, Timor-Leste and Vietnam) over the 2000-2021 period to demonstrate the suitability of probabilistic ML to predict crop vegetation health at increasing forecasting horizons between the measurement of the weather predictors and the VHI response. The choice of this particular region is motivated by several factors, including but not limited to the vast potential for agricultural production in the region, the high vulnerability to climate change and climate extremes and, the related strong dependency on climatic conditions of their food supply chains. In our study, VHI probabilistic predictions are then statistically linked with historical rice price monthly time series in each country to provide evidence of the relevance of VHI for economic policies for food security and evaluate the food price and expenditure impact of a vegetation health shock. This is of utmost relevance for the development and policy-oriented application of data-driven earlywarning systems. The analysis thus contributes to advancing the multidisciplinary literature linking remotely sensed vegetation health monitoring, probabilistic ML, and food security objectives.
The remainder of the paper is structured as follows: Section 2 provides a theoretical background on vegetation health monitoring, the role of remote sensing for food security, and on probabilistic machine learning. It also describes SEA, the pilot region under analysis. Section 3 then describes the methodology designed and implemented to generate probabilistic VHI forecasts and link them to food prices time series. Section 4 presents the results of the analysis, which are then commented in Section 5. Finally, conclusions and future work are discussed in Section 6.

Vegetation health measurement
The literature on the monitoring of the agricultural sector has taken great advantage of satellite technologies, and several indices have been developed to characterise land cover with a specific focus on vegetation health. As defined by Kogan (Kogan, 2019) [p.52]: Vegetation health is the method designed to derive vegetation condition (favorable, unfavorable, normal, etc.) or health (healthy, unhealthy, stressed,etc.) in response to changing weather (precipitation, temperature, and others) in each ecosystem and climate.
Remote sensing instruments detect light within the electromagnetic spectrum of an object when it is illuminated. Depending on the source of energy used to illuminate the object, there are two ways to collect remote data. Actively, that is, a sensor device emits a signal on the object of interest, and the sensor captures its reflection. Passively, when the natural sunlight reflection of the object is used (Woodhouse, 2017). An import implication for measuring vegetation health using remote sensing instruments is that plants absorb most of the light as part of photosynthesis. Specifically, vegetation absorbs the second-highest amount of solar energy in the visible range of the solar spectrum, which turns to increase the photosynthetic rate and green mass accumulation. Hence healthy plants reflect more near-infrared (NIR) than visible red (RED). When a plant becomes unhealthy, it reflects more RED and less of the NIR. This finding brought to the introduction of the Normalized Difference Vegetation Index (NDVI), measured as the ratio of the difference of the red RED and NIR radiance over their sum as in Eq. (1) (Kriegler et al., 1969).
Based on the NDVI, (Kogan, 1987;Kogan, 1990) proposed three alternative indexes. The Vegetation Condition Index (VCI), the Temperature Condition Index (TCI) and the Vegetation Health Index (VHI). The VCI is expressed as the relation between the value of the NDVI during a chosen composite period of analayis to the long-term minimum NDVI, normalized by the range of NDVI values calculated from the long-term record of the period of interest (Eq. (2)). VCI evaluates the current vegetation health compared to the historical trends and is designed to disentangle the weather-related component of the NDVI from the other environmental components. The TCI is similar to the VCI, but it is based on the brightness temperature (BT) (Eq. (3)) extracted from the thermal band (Ch 4) of the Advanced Very High-Resolution Radiometer (AVHRR) 1 and it captures the vegetation stress induced by abnormal temperature fluctuations and excessive wetness. Finally, the VHI combine VCI and TCI to assess the total moisture and thermal impacts on vegetation health (Eq. (4)). The three indices range from 0 (very unhealthy) to 100 (very healthy), with values from 40 to 0 indicating increased vegetation stress and 60 to 100 suggesting better vegetation health conditions.
The three Vegetation Health indexes developed by Kogan have been used extensively in the literature for monitoring the vegetation activity in response to weather-related drivers such as drought (Zuhro et al., 2020;Baniya et al., 2019;Kamble et al., 2019;Liang et al., 2017;Marufah et al., 2017;Masitoh and Rusydi, 2019;Pei et al., 2018) and to evaluate crop production (Kogan, 1990;Orlovsky et al., 2010;Rahman et al., 2009). The findings in the literature highlight that VHI correlates well with meteorological drought and agricultural drought in monsoonal rainfall areas (Marufah et al., 2017) and, most importantly, is useful to predict the yield of several grain crops such as corn in China (Kogan et al., 2005), wheat in the USA (Salazar et al., 2007), and rice in Bangladesh (Rahman et al., 2009) several months in advance of the harvest with considerable implications for food security. On the other hand, VHI was found to perform poorly in high latitude regions due to the assumption of inverse correlation between NDVI and BT that is not met in northern ecosystems where increasing temperatures support plants growth (Karnieli et al., 2006).

Food security and remote sensing
In the comprehensive definition of the Food and Agriculture Organization (FAO) of Food Security (FAO, 2003), some key features of what is to be in a"state of food security", are far from having being globally reached. A constant, undeniable and powerful force that limits access to"sufficient, safe, and nutritious food", is represented by extreme climate variability. In a recent work Hasegawa et al. (2021), suggest that regions heavily affected by climate extremes such as South Asia, might require to triple the current food production to compensate for the impact of climate change. An increase in production is challenging considering that since 1961 the total agricultural productivity has decrease by 21% globally and between 26 and 34% for regions closer to the tropics (Ortiz-Bobea et al., 2021). On the top of this, the"physical, social, and economic access" to food has been additionally hindered by the global COVID-19 pandemic, stressing the fragility of the world food supply chain as a consequence of the protracted crisis (Aday and Aday, 2020;FAO, 2020).
The complexity of these forces and the difficulties in providing reliable measurements and predictions of food security, have proven to be a stimulus for the development of a wide variety of systems for food security prediction and early warning systems based on different proxies for food insecurities. Using a variety of methodologies such as expert judgment (Funk et al., 2019) and data-driven approaches (Westerveld et al., 2021), the scientifically literature has contributed significantly to the way the status of food security is monitored in any given region of the world. Remote sensing have been widely used by practitioners and researcher as a costeffective quantitative and non-destructive source of data, to monitor and detect the evolution and health of crops over large ares and correlate it with crop yields as primary proxy for food security.
As highlighted by the in-depth review of Karthikeyan et al. (2020), there is a clear and robust relationship between crop yield and remotely sensed vegetation indices with correlations ranging from 60% to 80%. Aside from the use of VHI, other important indices have been explored in the literature such as the well-known (NDVI) (Rouse et al., 1974) and its variations such as the enhanced Vegetation Index (EVI) (Liu and Huete, 1995), the Green-Red Vegetation Index (GRVI) (Tucker, 1979), the Soil Adjusted Vegetation Index (SAVI) (Huete, 1988) and the Red Edge Position (REP) (Filella and Penuelas, 1994). In a recent work Jeong et al. (2022), used a combination of vegetation indices, weather and geographic variables to predict rice yields in North and South Korea within a Deep Learning framework. The target variables was reported rice yields collected at different scales by the local government and international organisations such as FAO as a part of a crop process-based model designed to simulate crop growth. The resulting model reached an R-squared of 0.86. In line with the results in Korea, the study of Son et al. (2020), demonstrated the effectiveness of ML algorithms for regional rice yield predictions based on Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data in Taiwan. Using a less computational intensive approach Nazir et al. (2021), found results ranging from 0.62 to 0.83 based on a Partial Least Square Regression in Pakistan on multiple indices including NDVI, EVI, SAVI and REP). On the other hand, in the review of Wen et al. (2021), the authors points out how some of these indices were found to be inaccurate in characterising crop health under drought and salinity stress and plant functional traits to be better in characterising the impact of these stressors.
It has to be noted that while all these studies indicate the robustness of the relation between crop yield and remote sensing and its direct relation with food security monitoring and prediction, none of them have explored the economic implication and relevance for further food security analysis. This goes in line with the review of Kubitza et al. (2020), where the authors stress the preponderance of studies limited to the biophysical sphere.

Probabilistic machine learning
Uncertainty and its quantification are fundamental aspects of any environmental analysis, and they are now crucial aspects in the face of climate change and its complicated consequences. Epistemic uncertainty (caused by limited data and knowledge) and aleatory uncertainty (randomness or variability in the underlying variables) are virtually present in all projections of a future phenomenon (Helton et al., 2010). In the last decade, researchers in the field of Probabilistic Machine Learning have developed more sophisticated algorithms that provide reliable predictions and account for the inherent uncertainties of Machine Learning algorithms (Duan, 2020;Ghahramani, 2015;Gneiting and Katzfuss, 2014;Kabir et al., 2018;Meinshausen and Ridgeway, 2006). By shifting from a pointwise estimation focused only at the conditional mean of a response variable to a distribution estimation aiming at providing information about the full conditional distribution of the outcome variable of interest, the notion of uncertainty is explicitly integrated in the prediction while providing more comprehensive information to policy-making.
To clarify how decision making and our understanding of a phenomenon of interest could benefit from moving beyond the conditional mean, it is worth introducing some notations. 2 Let Y and X denote a real-valued outcome variable and a set of covariates, respectively. A typical ML algorithm, produce a point estimate as in Eq. (5).
Eq. (5) can be interpreted as way of answering the question "What value Y would take on average, given X?". Answering this question can provide important insights, but -at the same time -such insights will regard exclusively the mean of the outcome of interest Y, while other information inherited from the data will be disregarded. A relatively simple, yet highly informative approach to overcome this limitation is to focus on the quantiles distribution of Y give X instead of the simple mean. Formally this is equivalent to Eq. (6).
Where the conditional distribution function F(y | X = x) is equivalent to the α quantile (Q α ) given X = x. This leads to prediction intervals as a way to encapsulate information about uncertainty into the model. Indeed, a prediction interval for Y given X, with a desired width α, can now be obtained directly from the data. For example, if the 90% prediction interval I(x) of Y for X = x is sought, then: In the classic parametric quantile regression, obtaining the quantities of interest is achieved by minimizing the expected loss E(L α ) (the analogous of minimizing the expected squared error loss for classic regression). If the loss function L α for 0<α<1 is defined by the weighted absolute deviations as in Eq. (8), then the objective is to minimize the expected loss E(L α ) as per Eq. (9).
In a non-parametric setting, Meinshausen and Ridgeway (2006) have shown an alternative approach to retrieve the conditional quantile distribution using the Random Forests (RF) algorithm of Breiman (2001). Quantile Random Forests (QRF) generalise RF by storing all observed responses at the leaf level of each tree, instead of just the mean, allowing to measure empirical quantile estimates. In other words, the classic conditional mean estimated by RF as the weighted mean of the observed outcome of interest, is replaced with its weighted distribution.
The predicted intervals give us an important additional information regarding the possible variation around a new predicted value. The correspondent width of the intervals can be interpreted as a simple yet highly informative measure of accuracy of the model, with wider intervals indicating a higher degree of uncertainty in the predicted value. A more sophisticated measure of performance of a distribution prediction can be found in the family of scoring rules (Gneiting and Katzfuss, 2014;Brier et al., 1950;Gneiting and Raftery, 2007;Krueger et al., 2016). A scoring rule function S(F, Y) is a measure of accuracy of a distribution prediction F given the observed outcome Y. In this paper, the Interval Score is considered, a proper scoring rule designed to score quantile predictions based on the work of Gneiting and Raftery (2007) which is computed as per Eq. (10).
Where Y is the true value of the outcome of interest, l and u denotes respectively the lower and the upper quantiles of the range defined by the value of α.
The usefulness of the QRF algorithm has been investigated in a number of contribution in the environmental literature. Sanderman et al. (2018) used QRF to provide probabilistic predictions of carbon stock in Mangrove finding that the highest level of uncertainty in the predicted carbon stocks at 1 m depth was on average 40.4% of the organic carbon stocks. Taillardat et al. (2016) compared QRF with a parametric probability density function model to forecast surface temperature and wind speed, using hourly data and lagged features from 3 up to 54 h. Their results highlight the suitability and better performances of QRF with sharp and reliable probabilistic predictions. Vaysse and Lagacherie (2017) Compared Regression Kriging for sparse digital soil mapping with QRF. The authors suggest that while the two algorithms reached similar performances on the point-wise prediction of soil properties, QRF is preferable over Regression Kriging since the former provides ore accurate and interpretable probabilistic prediction than latter.

Study area
The SEA region is characterised by a tropical climatic with monsoons and a substantial amount of rainfall alternated with period of dryness. Variation in rainfall an temperature between countries (Table 1) can be observed and have been related with El Niño/Southern Oscillation (ENSO), causing differences in the timing of the rainfall phase between the north-west of the region and the rest of the area.
According to Eckstein et al. (2021) the average temperatures in the SEA have been increasing rapidly since 1960 bringing Myanmar, the Philippines, and Thailand among the most affected countries in world by climate change in the period 2000-2019. Panel A of Fig. 1 shows the 2010-2020 average temperature change with respect to the long-term historical average in each country, highlighting that in most countries the one°C warming threshold has already been surpassed, with a peak of 1.3°C in Myanmar.
Turning to food production, the main commodity are represented by rice, maize, coffee, cocoa and palm oil. Rice is the single most important crop in the region. SEA countries have already experienced extensive losses due to extreme weather and climate change (Zhongming et al., 2021). Panel B of Fig. 1 shows the Coefficient of variation of rice yield coefficient of variation of rice yield (units of production per land cultivated) according to the country-level statistics for the 2010-2020 period provided by the FAOSTAT database (FAO, 2022). This reveals significant heterogeneity in the stability of rice productivity, with the greatest fluctuations observed in Cambodia, Laos, and Malaysia. Among the listed countries, Cambodia, Myanmar and Laos have been labeled as the "poorest of the poor" (Food &amp and of the United Nations, 1998). Overall, the percentage of population living on less than $5.50 a day, ranges from 2.7% in Malaysia to as high as 94% in Timor-Leste (Bank, 2021). The number of undernourished people in 2020 increased by 6% in the region most probably due to the pandemic (Zhongming et al., 2021). Finally, four of the country under analysis have a Global Food Security Index (GFSI) score (The Economist, 2021) falling below 60 with Laos, Cambodia and Myanmar below 50.  distance between the target VHI measurement and the predictor variables measurement: two-week, one-month, two-month and three-month ahead.

Modelling framework
The QRF models (a model for each of the nine pilot countries in SEA, for each prediction horizon) is trained on 15 years of weekly data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and tested on the following five years of sample (2016-2021). Finally, a panel regression model is estimated to appraise the statistical link between country-specifc variations in the VHI predictions and monthly rice prices. The estimated coefficients are used to derive implications for food expenditure and security.

Data sources and processing
3.2.1. VHI, climate, and satellite data The analysis is based on a combination of remotely sensed and atmospheric data detailed in Table 2 from the first week of April 2000 to the last week of July 2021. The primary outcome of the analysis is the VHI over cropland produced from the National Environmental Satellite, Data, and Information (NESDIS) part of the National Oceanic and Atmospheric Administration (NOAA). The data is available from 1981 with 4 km spatial and 7-day composite temporal resolution. To capture the long-term contribution of climate on vegetation health, yearly means of precipitation and temperature are included as covariates. For the short-term weather contribution, monthly rainfall, temperature, humidity, solar radiation means, cloud coverage and the total amount of rain by month are used. Raw covariates data are masked over cropland pixels on the Google Earth Engine platform (Gorelick et al., 2017) using the GFSAD1000 gridded cropland product (Teluguntla, 2015); then, the mean value of each variable is extracted for each day of the study period over SEA.
The masked data are than imported into the R scientific programming environment (Team, 2020), where additional processing steps are carried out to generate a set of features that could enhance the accuracy of the model for each prediction scenario, and capture trends and seasonality in the time series. This includes lagged features of each variable, rolling sum of rain and rolling mean of rain, temperature, humidity and solar radiation over the previous month. Year, month, week of the year, and sine and cosine function of time (for week, month, and year), are included as additional predictors (see Table 3).

Food security analysis data
To examine the relevance of the VHI indicator for food security forecasting purposes,additional statistical analysis based on monthly time series of the wholesale prices of rice in seven countries of East-Asia is carried out.  The price data is drawn from the FAO GIEWS Food Price Monitoring and Analysis (FPMA) tool 3 historical records, retrieving the longest possible span of time for each country corresponding to the VHI period under analysis. In this food security analysis, rice is considered as it is the first staple crop produced and in the region (FAO, 2022), as well as one where historical monthly price data is consistently available across the countries considered. The rice price time series are visualised in Fig. SI-3.
Then, for each of the countries investigated, statistics on the yearly average consumption of rice per capita are retrieved from the OECD-FAO Agricultural Outlook (OECD, 2020) and on the yearly average food expenditure from the USDA Economic Research Service. Combining these statistics with the results of the VHI-rice price statistical analysis, the potential impact of a shock of the VHI (such as a standard deviation decrease from its average long-run value) is simulate to evaluate its impact on household food expenditure.

Probabilistic forecasting of VHI
To explore the suitability of the QRF algorithm (Meinshausen and Ridgeway, 2006) to produce reliable probabilistic predictions of crop health for use in early-warning systems, a modelling framework based on increasing prediction horizons between the environmental predictors and the VHI is defined, as depicted in Fig. 2. More specifically, instead of providing a single estimation over the full SEA region as whole and using a single prediction horizon period, for each country four different prediction leads are considered, corresponding to multiple forecast horizons, namely, two weeks, one month,two months and three months ahead prediction. These lead to a total of thirty-six models trained over about sixteen years and tested over a period of five years. While our main interest is over the predicted quantiles intervals at a probability of 0.1 and 0.9, also provide point-wise predictions of RF are provided.
In this study, the QRF algorithm is implemented using the R packages ranger (Wright and Ziegler, 2015) and caret (Kuhn, 2008). To find the optimal set of hyperparameters, the resample procedure for times series proposed by (Hyndman and Athanasopoulos, 2013) is adopted. To score the result of the probabilistic prediction, the quantiles scoring rule metric as implemented in the package scoringutils (Bosse et al., 2020) is considered, while for the point-wise predictions classic metrics including Rsquared (R 2 ), root mean square percentage error (RMSPE), and mean absolute percentage error (MAPE) are adopted. 4

Linking VHI forecasts with food prices and expenditure
To link time series of wholesale rice prices with the VHI probabilistic predictions, a monthly unbalanced panel (due to price data availability varying across countries, see Fig. SI-3) dataset is assembled. Based on these data: 1. Fixed effects panel data regression is carried out to absorb unobserved country-specific and time-specific heterogeneity from the effect estimation. In particular, both country fixed effects, and month and year fixed effects, as well as their interaction are included. Here: • country fixed effects absorb national market specificities; • month fixed effects capture the seasonality dynamics of prices; • year fixed effect are representative of regional to global variations in prices which affect all countries in the panel; • finally, month by year interactions capture the interplay between the two latter dynamics.
Fixed effects regression thus aims at isolating the true link between VHI predictions and rice prices while controlling for additional underlying trends affecting prices. The models are run on a dataset ensemble of twoweek up to three-month ahead VHI predictions.
2. As predicting variables, the middle quantile of the predicted VHI is included, as well as an array of lagged values (1-3-6-12-24-48-72 months) to control for the short-to long-term impact of vegetation health on prices. 3. Different dependent variable timings are tested for, starting from the current value of prices, going back to lagged values of prices by 1, 2, and 3 months with respect to the VHI prediction. The reason for testing specifications with lagged dependent variables is that lags can capture the anticipation effect of prices, such as market expectations over future food prices 4. Finally, a control variable measuring the prediction interval width at each t is also included.  In mathematical terms, the following is estimated: P cmy ¼ VHI cmy þ VHI cm À 1 y þ VHI cm À 3 y þ VHI cm À 6 y þ VHI cm À 12 y þ VHI cm À 24 y þ VHI cm À 48 y þ VHI cm À 72y þ width cmy þ γ c þ ψ m þ μ y þ ψ m ⋅ μ y þ ε cmy (11) where: • P are wholesale rice prices • c, m and y represent country, month and year, respectively • VHI is the middle quantile of the Vegetation Health Index prediction • width is the VHI prediction interval width • γ are country fixed-effects • ψ are monthly fixed-effects • μ are yearly fixed-effects • ε is the error term Finally, the estimated regression coefficients are used to evaluate the impact of a shock in the (predicted) VHI on rice prices, and thus on rice expenditure given a representative yearly level of per-capita consumption based on historical statistics. To achieve this, a shock S c is simulated: Here, for each country c, n describes the number of SD (standard deviations) from the mean VHI values to the lowest percentile of each distribution, i.e. for the vegetation to become stressed (falling in the range 0 to 40). The VHI shock S c is then multiplied by the regression coefficient from Table 4 (considering the mean value of VHI regression coefficients from specifications 1-4) to derive the impact on rice prices, and divide it by the most recent yearly mean price in each country to get the relative (%) price impact. Finally, based on country-specific average household yearly rice consumption and food expenditure statistics, the potential impact of the shock S c on total household food expenditure due to this potential increase in rice prices is calculated. Table SI-3 reports, for each country, the yearly average household food expenditure, the mean rice price for the most recently available year, and the value of one standard deviation in the VHI. In addition, it is useful to compare this quantity with the month and country specific observed variations in the VHI index between 2000 and 2021 plotted in Figure SI-1.

Probabilistic forecasts of cropland vegetation health
Contrarily to point-wise estimates, a prediction interval offer a way to assess the reliability of a prediction by its width. Tighter prediction intervals indicate less fluctuations of the response variable and consequentially, a more reliable prediction. Fig. 3, displays the range of potential outcome values for VHI against the observed values (yellow line) for the four prediction horizon for each of the nine countries analysed in the five-year validation period (see Methods). Wider fluctuations are observed especially for Malaysia, Myanmar and Thailand even on the two weeks prediction horizon (dark blue ribbon). Overall, the intervals for two weeks and one month prediction horizons are tight around the observed values of VHI while the intervals for two and 3 months are very large with a tendency to under-predict the observed outcome especially for Malaysia and Myanmar. On the other hand, over-prediction is more pronounced for Laos in the period 2019-2020 and at the end of the time series for Timor-Leste.
To quantify the prediction intervals goodness-of-fit, the Interval Score (IS) metric is considered (Gneiting and Raftery, 2007). The IS has attractive properties because it balances coverage and length of the prediction intervals with lower value of the metric indicating better estimated prediction intervals. Fig. 4A reports the mean value of the calculated IS metric, across the 5-year testing period in the analysis, while Fig. SI-2 shows the complete time series for each week in the testing period.
With a prediction horizon of two weeks, for most of the countries under analysis the MIS is below the value of one, suggesting high reliability on average. For each increase in the prediction horizon, all the predictions depart from 0, but Indonesia and the Philippines are the best performing. The predicted intervals for Myanmar, Cambodia, Laos and Thailand are the most sensitive with the one for Cambodia triplicating from the initial prediction horizon of two weeks to the three months horizon.
While the focus and the core novelty of our approach is its probabilistic nature (allowing to produce prediction intervals) the conventional pointwise prediction accuracy metrics are also considered for comparability with the existing literature. Fig. 4B reports the training and testing accuracy metrics -as measured by the R 2 -of the point-wise prediction of the QRF models, 5 by country and prediction horizon. While the results on the training are all very close to 100%, it is evident that increasing the prediction horizon negatively impacts the accuracy of the prediction. Within a prediction horizon of maximum one month the results on the test are still within the range of 60% to 80% which is in line with the results found in the existing literature. An exception is Indonesia, where the accuracy remain constantly above 75% even with a prediction horizon of 3 months. On the other hand, the worst performing prediction is for Malaysia with an accuracy as low as 35% and 18% for two and three months prediction horizons respectively (additional metrics including RMSPE and MAPE are displayed in Table SI-1). Table 4 Panel fixed-effects regression results for the impact of the VHI predictions on national monthly rice prices. Estimated regression coefficients of Eq. (11) with standard errors in parentheses. L1-L2-L3: one, two, and three-month lagged variables.  Table 4 presents the results of the fixed-effects regression models linking VHI predictions with monthly national rice price time series (see Methods). In particular, models with four different dependent variables are tested: rice prices, and the one, two, and three-month lag of rice prices. The lagged specifications are considered to test whether an anticipation effect exists, i.e. if prices adjust in advance to future variations in the VHI index, such as based on future yield expectations. In the regression table, the variable VHI identifies the middle quantile of the probabilistic VHI estimates, whilst a control variable measuring the prediction interval width at each t is also added (see Methods).

Linking cropland health predictions and food security
The results show that the VHI has a strong and significant impact on rice prices in all the considered specifications, and that this effect increases in magnitude as the order of the lag of the dependent variable increases. This means that rice prices are adjusting in advance based on future vegetation health (and, presumably, yield) expectations. Additional results based on the actual VHI values, rather than the predictions, are found in Table SI-2. An additional important finding is that the results based on our VHI ahead predictions find good agreement with the results of the models based on the observed VHI values in terms of their magnitude and statistical significance. This result suggests that VHI predictions are suitable for food price change prediction analysis, such as in early warning systems.
To complement the analysis and quantify the regression coefficients in Table 4, Fig. 5 reports the results of the simulation of a potential shock S c in the VHI (see Methods for details on the definition of the shock magnitude and Table SI-3 for numerical details) and its impact on rice prices, and thus on total household food expenditure. The graphs show that the simulated shock can have a large impact on national rice prices (an average increase of 16.5% of the current mean price), although with great variability, with a minimum of 7% in Laos and a maximum of up to 25% in Myanmar and Vietnam. This price change could have a potential direct impact on total household food expenditure: on average, a 2.5% growth, ranging between a minimum of 1.4% in the Philippines and Thailand, and a maximum of 4.1% in Cambodia.

Discussion and policy implications
Most of the studies in the literature have focused on point-wise estimations, setting aside considerations on the inevitable uncertainty that comes with any prediction. Moreover, compared to most previous studies, the results of this study are based on a large and heterogeneous region of the world in both vegetation, climate and food insecurity terms.
Given the innovative methodological aspects and the wider geographical scope of this work, it is challenging to compare its results with previous studies. To mention the most similar ones, the analysis of Nay et al. (2018) achieved R 2 values raging from 0.76 to 0.86 with a two-week ahead prediction horizon, but it was limited to agricultural land in California and Sri Lanka. Similarly, Perera et al. (2020) achieved an accuracy of 0.93 for the MENA region, but their analysis solely focuses on near-real-time prediction. Lees et al. (2022) reached an average accuracy of 0.83 for one-month ahead prediction and 0.95 for three-month ahead prediction using a Deep Learning architecture but focusing on a circumscribed area in Kenya. The most directly comparable study is the work of Salakpi et al. (2021), which compares different prediction horizons within a Bayesian framework that naturally provide uncertainty in the predicted values. More specifically their 1.5 to 3 months-ahead prediction R 2 ranges between 0.94 and 0.62 but also in this case, the focus remain on a limited region.
In this context, the findings of this paper thus have important implications from a policy-making point of view as they provide a proof-ofconcept for the design and implementation of early-warning systems and contingency plans for crop health and food prices based on satellite data and probabilistic machine learning. In particular, the results are based on openly available remotely sensed data and they are able to capture nonlinear and accumulation dynamics of weather variables. In addition, the proposed framework can be adapted to virtually any geography and crop type. These are attractive features for decision-makers looking for readily updated, inexpensive food-related forecasts (Ehlers et al., 2021). Such system can substantially reduce the need for continuous field data collection to keep track of current and expected changes in cropland vegetation health. In turn, the ability to anticipate national food price variation following weather shocks can have important implications for governmental food trade and stocking decisions, as well as for short-run risk hedging decisions in the commodity stock market. Additional considerations are related with the trade-offs between the short and the long run impact of food security policies. The proposed framework and results are mostly relevant for the short term policy response to support households against food prices shocks (Ansah et al., 2021). On the other hand, this strategy might eventually converge to more structural measures and prove to have long-term moderating effects on climate shocks on household dietary and health conditions.

Conclusions
In this paper, a statistical framework to improve the understanding of the nexus between weather variability, crop health, and food security is designed and implemented. QRF -a ML algorithm designed to estimate quantile values and thus produce probabilistic predictions -is used to train a set of models with twenty years of weekly multi-source and highresolution environmental satellite data in nine countries of SEA.
The work stresses the importance of accounting for uncertainty in prediction when using ML and how this is needed for more actionable insights for policy-making under rapidly changing climate conditions. To concretely show the added value of the probabilistic prediction of future crop health conditions, the paper examines the impact of weather shocks on household food expenditure and finds that the VHI predictions are useful to anticipate food price and expenditure shocks.
Results shows that the algorithm can accurately learn and reproduce the historical patterns of vegetation health based on weather variability, and thus successfully anticipate crop health under different prediction horizons. In a second step, a follow-up econometric analysis to evaluate the link between the vegetation health predictions and food security, relating the VHI dynamics with time series of national monthly rice prices is carried out. Results show that VHI predictions are suitable to anticipate food price and expenditure shocks.
Irrespective of this strong policy-support potential, a number of remarks must be discussed in relation to the results. First, while findings are robust and in line with the performance and results of the wider literature that links satellite imagery with crop health, a limitation of the analysis and its findings stems from the potential error in remotely sensed variables used in the analysis and generally in the data granularity. It also has to be noted that smooth effects between the set of variables under analysis are not well captured by RF hence any further analysis could benefit from the use of a local linear correction as proposed by Friedberg et al. (2020) but at the cost of higher computation and time requirements. Second, in the VHI shock analysis, the estimated household food expenditure increases are solely based on rice prices variations. Hence, they are best interpreted as lower-bound values of the overall impact of a VHI shock on food security. Conversely, these rice price growth impacts on food expenditure do not consider potential substitution dynamics, e.g. temporary switches to other foods less affected by a potential VHI shocks. In other words, the price elasticity of rice demand is assumed to be inelastic (ε = 0).
Overall, this analysis contributes to the growing literature focused on using remotely sensed data and ML algorithms to improve food security analysis. The proposed framework sets the ground for the development of early warning food security systems inclusive of an uncertainty component that could support programs and policies aiming at boosting the resilience and adaptation to climate change in vulnerable countries. In this perspective, future research could extend the here proposed framework to other geographies, consider additional crops, and represent food trade and crop substitution dynamics.

Code and data availability
The code to replicate the analysis and the figures are publicly hosted at https://github.com/athammad/VHI_SEA. The repository includes references to retrieve the input data.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.