Reliable quantification of uncertainty in crop-yield projections under climate change is critical for global food security. Beyond directly influencing food availability, crop production has broader socioeconomic and political implications. Fluctuations in agricultural outputs can trigger global price shocks, leading to social unrest and political instability, as exemplified by the Arab Spring, which was accelerated by climate-induced declines in wheat production and resulting price spikes (Johnstone and Mazo, 2011). Additionally, agricultural forecasts serve as vital inputs to integrated assessment models (IAMs), shaping regional investment policies and adaptation strategies (Nelson et al., 2014). Accurate yield forecasts support effective adaptation strategies, mitigate risks to agricultural production, and guide policymaking in a rapidly changing climate (Rosenzweig et al., 2014;Asseng et al., 2013;Müller et al., 2021). Crop yields respond to numerous interacting factors, including climate extremes, biophysical processes, technological advancements, and socio-economic conditions. Together, these factors introduce significant complexity and uncertainty into yield forecasts (Challinor et al., 2018;Müller et al., 2021). Addressing these uncertainties rigorously remains a fundamental challenge for climateimpact assessments and agricultural policy. Earlier research has examined multiple approaches to reduce uncertainty in yield projections. A common method is to enhance crop-model representations of key physiological processes, such as crop sensitivity to temperature. These improvements have notably reduced simulation errors and uncertainty in crop-yield forecasts (Wang et al., 2017). Zhang et al. (2015) explicitly compared uncertainties from climate and crop models in maize yield predictions. They demonstrated that, over certain future periods, variability in climate projections often sur-passes uncertainty arising from differences among crop models. More recently, studies have combined machine learning with global gridded crop models (GGCMs) and detailed environmental datasets. These hybrid methods successfully capture complex yield responses to extreme weather events, pests, and diseases. Consequently, they reduce uncertainty and reintroduce stressors often overlooked in earlier models (Li et al., 2023b;Heinicke et al., 2022). Despite these advances, most existing approaches remain focused primarily on refining model structures or incorporating additional data. Few studies have comprehensively characterized the full spectrum of uncertainty within a unified framework. Even fewer explicitly consider regional differences or systematically integrate uncertainties from multiple sources simultaneously. Ensemble-based studies (Asseng et al., 2013;Müller et al., 2021) highlight the value of multimodel comparisons, but typically assume uniform error distributions across regions. They also often rely on deterministic or semi-probabilistic methods, potentially underestimating or mischaracterizing the true scale of uncertainty (Franke et al., 2020). This leaves a significant gap: reliable and fully probabilistic quantification of yield uncertainty remains incomplete. To fill this gap, our study introduces a Bayesian hierarchical model embedded in a fully probabilistic framework. Our approach enhances the precision and reliability of uncertainty quantification for crop-yield projections. We build upon the recent multivariate autoregressive crop-yield growth model (Li et al., 2025) and introduce a novel extension that explicitly models error variances separately by country. This improvement relaxes the traditional assumption of homogeneous errors, better capturing real-world regional variability. Unlike previous approaches focused primarily on model refinement and data integration (Wang et al., 2017;Zhang et al., 2015), our framework fundamentally advances the quantification and interpretation of uncertainty itself. In our Bayesian framework, we jointly account for the uncertainties stemming from climate-model projections, emissions pathways, and the parameterisation of the crop model. This approach provides detailed annual probabilistic forecasts and explicitly quantifies tail risks associated with extreme yield outcomes. Specifically, we measure extreme events using two established metrics from risk management: Value-at-Risk (VaR), the yield threshold that our model expects will be met or exceeded with 1α probability in that particular year, and Expected Shortfall (ES), the average yield conditional on falling below that VaR level. Because both numbers are estimated directly from the Bayesian posterior predictive distribution, they carry forward the full uncertainty from climate models, emissions pathways, and crop-model parameters, giving stakeholders a clear, year-by-year gauge of downside agricultural risk. We validate the accuracy of our uncertainty quantification using historical crop and climate data. Through detailed global wheat production analyses, our model consistently outperforms existing methods. It provides robust, interpretable forecasts, offering policymakers clear, risk-informed guidance. This comprehensive framework thus enables more resilient adaptation planning and more informed agricultural decisions in the face of future climate extremes. The remainder of this paper proceeds as follows. Section 2 describes the climate and crop datasets and explains our data processing methods. Section 3 details our modeling framework. First, we review the existing multivariate autoregressive crop-yield growth model (MAR-X). Next, we describe its Bayesian implementation. Finally, we extend this model with a hierarchical structure to allow for country-specific error variances. Section 4 presents the Sequential Monte Carlo (SMC) algorithm used for parameter estimation and model selection. Section 5 provides our empirical results. Specifically, Section 5.1 evaluates historical forecast accuracy and tail risk. Section 5.2 shows probabilistic yield projections under various future climate scenarios. Section 6 concludes the paper by discussing key implications of our results. This study utilizes the same climate and crop datasets as described in Li et al. (2025), where a comprehensive description of the data sources and processing methods is provided. For completeness, we briefly summarize the relevant datasets below. We use historical climate data from the Climate Research Unit (CRU) Time Series version 4.05 (CRU TS 4.05; http://www.cru.uea.ac.uk), which provides monthly observations of surface climate variables, including temperature and precipitation, on a 0.5 • latitude-longitude grid from 1961 to 2018 (Harris et al., 2020). CRU TS is a widely used observational dataset that is constructed from weather station records using interpolation methods. It belongs to a broader category of climate reanalysis and observational reconstruction products, which aim to reconstruct past climate conditions by integrating diverse observational datasets (Kalnay et al., 2018;Dee et al., 2011). To generate future climate projections, we apply the QuickClim method proposed by Kitsios et al. (2023). QuickClim is a machine-learning-based emulator designed to efficiently replicate the outputs of complex global climate models (GCMs). A climate model is a computer-based simulation that represents how different parts of the Earth's climate system-such as the atmosphere, oceans, land surface, and ice-interact over time.
While GCMs provide detailed and scientifically robust simulations of climate responses under various greenhouse gas and land-use scenarios, they are computationally intensive. QuickClim addresses this limitation by rapidly approximating monthly climate variables, including temperature and precipitation, across a wide range of emissions scenarios. To align with standard climate-impact assessments, we focus on three Representative Concentration Pathways (RCPs)-RCP 2.6, RCP 4.5, and RCP 8.5-corresponding to low, intermediate, and high greenhouse-gas emissions. These pathways were adopted by the Coupled Model Intercomparison Project Phase 5 (CMIP5) as common inputs for global climate-model experiments (Van Vuuren et al., 2011). CMIP5 itself is an internationally coordinated ensemble of climate models that simulates future climate conditions under harmonized assumptions about greenhouse-gas concentrations and land use (Taylor et al., 2012). Each RCP traces a plausible evolution of emissions and land use that delivers a specified level of radiative forcing by 2100. Radiative forcing measures the net energy imbalance introduced into the Earth system by elevated greenhouse-gas levels and other agents, and it can be translated into effective atmospheric concentrations of carbon-dioxide equivalents (CO 2e ). Figure 5a plots those CO 2e trajectories for the RCPs alongside the ensemble of 100 QuickClim pathways used in our study. Annual historical crop production data spanning the years 1961-2018 are obtained from the FAOSTAT database (http://faostat.fao.org), covering national-level wheat production. Information about crop planting and harvesting dates, necessary for linking climate variables with crop growing seasons, is taken from the global crop calendar database of Sacks et al. (2010). Using these growing periods, we extract relevant climate variables during the crop growth stages. In cases where countries have multiple cropping seasons, we average seasonal climate conditions within each country before associating them with annual crop production records. Additionally, we apply a standard filtering process, removing outliers identified based on anomalously high volatility in annual crop production. This ensures robust estimates by minimizing potential biases arising from data inconsistencies or extreme events. While this work focuses on wheat in the empirical analysis, the proposed approach can be applied to other crops. The recently developed multivariate autoregressive model with exogenous climate variables for the estimation of crop yield log returns, (MAR-X; Li et al., 2025), takes the form of:

(2)

Here y t = [y 1,t , y 2,t , . . . , y K,t ] ′ denotes the log-returns of wheat yield for K crop-growing regions at year t. The parameters θ 1 , θ 2 , θ 3 , θ 4 are assumed to be scalars, ǫ t is the error term, and ∆T t and ∆T 2 t represent the differences in the annual mean temperature during the crop-growing period and its quadratic form, respectively, for the current year t and the previous year t -1 across K regions. Previous work (Li et al., 2025) evaluated multiple functional forms for the intercept term, θ 0 (t), and found that an exponentially decaying form (Eq 1) provided the best fit to historical data and improved forecasting performance for wheat. While the MAR-X provides valuable insights into estimating wheat yield log returns, it has several limitations that highlight areas for potential improvement. The model assumes a shared error term across all countries, which may not accurately represent the diverse sources of uncertainty in agricultural production. The shared error term overlooks country-specific factors that could influence yield variability, such as localized climate anomalies, policy changes, or infrastructure variations. This limitation potentially reduces the precision and robustness of yield predictions across different countries. Additionally, the model does not consider uncertainties inherent to the agricultural production model itself. This gap suggests a need for a more comprehensive approach that rigorously quantifies and incorporates these uncertainties to enhance the reliability of yield and production projections. The Bayesian framework requires the specification of prior distributions for each model parameter. For the intercept parameters, a and λ, we assume uniform priors bounded within specific ranges. The parameters a and λ in Eq 2 are sampled as: a ∼ Uniform(0, 0.5), λ ∼ Uniform(0, 0.2). (3)

The above prior ensures that both the initial technological advancements and the decay rates remain positive. The upper bounds are selected to include the non-linear least squares point estimates from the non-Bayesian MAR-X model, while maintaining a narrow range to enhance the efficiency of the Bayesian sampling process. The autoregressive coefficients, θ 1 and θ 2 (Eq 1), which relate to the lagged yield log-returns, are assigned normal priors with mean 0 and standard deviation 0.5:

The normal distributions of the above priors reflect the possibility of both positive and negative effects, and a small variance is selected to be consistent with the relatively modest point estimates reported in the existing literature. The temperature-related coefficients, θ 3 and θ 4 , are also modeled with normal priors, assuming mean zero and and standard deviations reflecting more conservative beliefs about their influence:

Finally, the variance of the error term, σ 2 , is assigned a commonly used non-informative inversegamma prior, which is a natural and appropriate choice for positive variance parameters:

The posterior distributions of the model parameters are obtained by combining the specified priors with the likelihood of the observed data. The model specification for each region i at time t is given by

where the conditional mean is defined as

The likelihood is based on the assumption that the error terms are independently and identically distributed as ǫ i,t ∼ N (0, σ 2 ), for i = 1, . . . , K regions and t = 1, . . . , T time periods. The density for each observation is then given by

with Θ = {a, λ, θ 1 , θ 2 , θ 3 , θ 4 , σ 2 } representing the model parameters.
Assuming that the observations are conditionally independent across regions and time, the likelihood for the full dataset is

The Bayesian framework for the MAR-X model (B-MAR-X) provides a flexible and probabilistic approach to estimating wheat yield log-returns while accounting for parameter uncertainty and variability. In order to address the limitations inherent in the single, pooled error specification of the MAR-X, we now introduce a more sophisticated framework. In this section, we extend our analysis by incorporating a hierarchical structure that assigns country-specific error distributions while still borrowing strength across the panel. By nesting the individual variances in a common hyper-prior, the model simultaneously (i) captures genuine cross-country heterogeneity in yield volatility, (ii) mitigates over-fitting in data-sparse settings through partial pooling, and (iii) delivers fully probabilistic forecasts that propagate both parameter and process uncertainty. This hierarchical structure improves information sharing across countries and yields more precise, well-calibrated uncertainty estimates-providing a firmer basis for data-driven food-security risk assessments. Details of the model formulation and prior specifications are provided below. The conditional mean of the crop yield log-returns for K countries at time t follows the MAR-X model specification given in Eq 1, where ǫ t = [ǫ 1,t , ǫ 2,t , . . . , ǫ K,t ] ′ (with ′ denoting the transpose) is a vector of i.i.d. error terms for each country. Unlike the standard MAR-X model, here we assume that the error term for each country follows its own normal distribution:

where σ 2 i represents the variance of the error term for country i. Rather than fixing the variance σ 2 i , we model them as random effects drawn from inverse-gamma distributions. Specifically, the hierarchical structure for the error variances is given by:

In this structure, α σ is a global shape parameter that determines the general behavior of error variability across countries, while β σ is a global scale parameter. The hyperpriors for α σ and β σ are defined as:

This hierarchical setup allows for flexible modeling of country-specific error variances, with global priors capturing the overall distribution of variability. Next, we present the likelihood function for the hierarchical Bayesian MAR-X (HB-MAR-X) model. Let Y = {y i,t } i=1,...,K; t=3,...,T denote the panel of observed crop-yield log-returns, and let the conditional mean µ i,t be defined as in Eq 8. For each time point t, define

Assuming mutually independent Gaussian innovations,

The likelihood contribution for each time point t, conditional on the structural coefficients θ = (a, λ, θ 1 , θ 2 , θ 3 , θ 4 ) and the vector of country-level variances σ 2 = (σ 2 1 , . . . , σ 2 K ), is given by:

Assuming conditional independence across time, the likelihood of the full sample is

This form first groups the independent Gaussian densities by time point and then multiplies across t. It highlights two features of the hierarchical MAR-X specification: (i) conditional independence across countries within each time slice, and (ii) heteroscedasticity introduced by the country-specific variances σ 2 i , which are themselves random quantities in the next layer of the hierarchy. The HB-MAR-X model retains the original econometric framework's simplicity while introducing the flexibility of country-specific random effects. Independent Gaussian shocks remain at the observation level. However, their variances are now modeled using a shared Inverse-Gamma distribution, with hyper-parameters estimated from the data. This structure aims to enhance the adaptability and robustness of the model by accounting for both global tendencies and local variations in crop yield responses. This layered design yields several practical benefits: it disciplines extreme variance estimates, improves out-of-sample predictive accuracy, and furnishes decision-makers with posterior predictive distributions that fully reflect ensemble uncertainties across climate scenarios and technological pathways. In short, the model offers a robust, policy-relevant platform for quantifying both global trends and local idiosyncrasies in future crop-yield dynamics. The Bayesian estimation procedure adopted in this study utilizes Sequential Monte Carlo (SMC), a powerful approach suited for hierarchical models where posterior distributions may be complex or multimodal. SMC efficiently approximates posterior distributions and computes marginal likelihoods, facilitating robust uncertainty quantification and model comparison. The detailed SMC algorithm used for parameter estimation and model selection is described fully in Appendix 1. To provide a visual illustration of how the uncertainties surrounding the crop yield model's unknown parameters are quantified and distinguished from point estimates, Figure A.1 presents the marginal posterior distributions for the B-MAR-X model described in section 3.2, alongside point estimates generated by the existing frequentist approach in Li et al. (2025). While the frequentist method produces point estimates without explicitly incorporating uncertainty in the error variance, the Bayesian approach jointly models both the parameters and the error vari-ance, resulting in a full probabilistic characterization of parameter uncertainty. The posterior distribution of the error variance is also included in the figure to illustrate this aspect. For all parameters of the B-MAR-X model, the non-linear least squares point estimates fall well within the corresponding Bayesian posterior distributions and align with their highdensity regions. This consistency between the two approaches supports the validity of the Bayesian estimates and highlights the effectiveness of the uncertainty quantification achieved through the SMC-based inference. For the HB-MAR-X model, instead of using a single variance parameter as in the B-MAR-X model, the posterior distributions of country-level variances and the associated joint parameters are estimated jointly. The accuracy of the uncertainty estimates is further evaluated in the results section. In addition to parameter estimation, Bayesian methods facilitate rigorous model selection through the computation of marginal likelihoods (model evidence). For the hierarchical Bayesian model introduced in section 3.3, the marginal likelihood involves integrating over both model parameters and hyperparameters. A complete derivation and description of this marginal likelihood, including its estimation via SMC methods, is provided in Appendix 2.
Specifically, for the proposed hierarchical model, the marginal likelihood P (y) of the observed data is defined as:

where Θ = {a, λ, θ 1 , θ 2 , θ 3 , θ 4 , σ 2 i | i = 1, . . . , K} denotes the model parameters and their country level variances, and Φ = {α σ , β σ } denotes the hyperparameters. The likelihood of the observed data y = {y 1 , y 2 , . . . , y T }, given the model parameters Θ, is:

where P (y t | •) follows a Gaussian distribution specified in section 3.3. The prior distribution of the model parameters, given the hyperparameters Φ, is:

The prior distribution of the hyperparameters Φ is:

The marginal likelihood P (y) for this hierarchical model can be efficiently estimated using SMC by incrementally computing the normalizing constants along an annealed sequence of posterior distributions, thereby providing a practical method for quantifying model evidence and performing Bayesian model comparison. To for all model parameters and error terms. These three variants are denoted as HB-MAR-X (a, λ), HB-MAR-X (a, λ, σ 2 ), and HB-MAR-X (all model parameters), respectively, as shown in Table 1. Moreover, as an alternative to employing hierarchical structures with hyperpriors in the HB-MAR-X (a, λ, σ 2 ) specification, we also consider a model in which priors are assigned directly to the country-specific time-varying intercepts and error terms without a hierarchical structure. This model is referred to as MAR-X-IND (a, λ, σ 2 ). The log evidence in Table 1 considers both the full set of 57 wheat-producing countries and a selected subset comprising 40 countries determined via our forward selection procedure. Across all model variants, there is a consistent improvement in model evidence when reducing the number of included countries from 57 to the selected subset of 40. This result highlights the benefit of selecting a subset of countries, demonstrating that removing countries that contribute primarily noise or irrelevant information significantly enhances model parsimony and predictive performance. Specifically, the HB-MAR-X (σ 2 ) model, which applies a hierarchical structure exclusively to country-level error terms, achieves the highest log evidence in both scenarios. This indicates that introducing country-specific variance structures substantially improves the model's capacity to capture heterogeneity in crop-yield data without excessive complexity. Simpler models such as the traditional B-MAR-X and the hierarchical model without variance terms (HB-MAR-X (a, λ)) achieve notably lower model evidence, highlighting the importance of modeling country-specific variances explicitly. Finally, the fully hierarchical model (HB-MAR-X with all parameters hierarchical) and the independent prior model (MAR-X-IND (a, λ, σ 2 )) yield significantly poorer results, particularly when using all 57 countries. This indicates that excessive hierarchical complexity or independent priors without hierarchical constraints may lead to overfitting and reduce overall model performance. To assess the suitability of an agricultural production model, for instance, one might consider its predictive capability. Building upon the Bayesian framework and SMC approach previously introduced, the posterior distribution π t (θ|y 1:t ) facilitates one-step-ahead posterior predictions for data at time t + 1. Within this fully Bayesian setting, the posterior predictive distribution for upcoming observations ŷt+1 is obtained by integrating over the model parameters, thereby naturally incorporating parameter uncertainty without asymptotic approximations (see detailed derivation in Appendix 3). The posterior predictive distributions can be approximated using the previously described SMC algorithm, generating weighted samples from the posterior. Point estimates, such as mean or median, and predictive intervals can then be straightforwardly computed from this weighted set. To assess the quality of the predictive intervals produced by the model, we evaluate the accuracy of the posterior predictive distributions using the squared error of the 95% prediction interval (PI). Specifically, for each country, we compute the squared difference between the empirical coverage rate-the proportion of observed data points falling within the 95% predictive interval-and the nominal coverage rate of 0.95. This provides a direct measure of how wellcalibrated the predictive intervals are, with smaller squared errors indicating better alignment between the forecast distribution and the realized outcomes. relative to B-MAR-X across many countries, with notable improvements in a few cases. In particular, HB-MAR-X exhibits a dramatic reduction in squared error for Australia, where B-MAR-X clearly miscalibrates. While the differences are more modest for most other countries, the HB-MAR-X model maintains consistently low and stable squared errors, suggesting that its hierarchical structure contributes to more reliable generalization and improved uncertainty quantification. We also rigorously evaluate the forecasting performance using leave-future-out cross-validation (LFO-CV) as proposed by Bürkner et al. (2020). Unlike traditional leave-one-out methods, LFO-CV explicitly considers temporal dependence, making it appropriate for sequentially ordered crop yield data. Predictive accuracy of the one-step-ahead forecasts is quantified by the expected log predictive density (elpd; Vehtari et al., 2017), comprehensively evaluating the entire predictive distribution rather than single-point predictions (details of this evaluation are provided in Appendix 3). The proposed HB-MAR-X (σ 2 ) model, which exhibits the highest log evidence as shown in Table 1, achieved a higher elpd of 8.69 compared to 7.02 for the existing B-MAR-X model, indicating superior predictive performance. This improvement suggests that incorporating hierarchical variance components enhances the model's ability to accurately forecast crop yields across the selected 40 countries. In our Bayesian framework for crop yield forecasting, the predictive distributions provide a foundation for assessing extreme agricultural risks. Borrowing concepts from financial risk management, we adapt the Value-at-Risk (VaR) and Expected Shortfall (ES) metrics to quan-tify the potential for adverse outcomes in agricultural production. Specifically, we use these measures both in terms of negative log yield returns (capturing relative yield losses) and in terms of low absolute yield levels, thereby making them intuitive and relevant to practical agricultural risk management. Let Y denote a random variable representing yield-related risk. This may correspond to the negative log yield return, where higher values indicate greater yield loss, or it may represent the actual yield level, where lower values correspond to poorer production outcomes.
Let α ∈ (0, 1) denote a chosen confidence level (commonly, α = 0.0 or α = 0.99). Value-at-Risk at confidence level α, denoted by VaR α (Y ), is defined as the α-quantile of the distribution of Y :

Here, x represents a threshold for the yield-related quantity Y , such that there is at least an α probability that Y does not exceed x. When Y is the negative log yield return, x corresponds to a worst-case loss level not exceeded with high probability. When Y is the absolute production level, x identifies a minimum production level that is met or exceeded with probability α. Expected Shortfall at level α, denoted by ES α (Y ) (also known as Conditional Value-at-Risk or CVaR), measures the average severity of outcomes in the tail beyond the VaR threshold. Assuming sufficient continuity of the distribution of Y , ES is defined as:

where E[•] denotes the expectation operator. This formulation is appropriate when higher values of Y correspond to worse outcomes, such as in the case of negative log yield returns. In this context, the Expected Shortfall represents the average of the worst 1α proportion of outcomes, capturing the expected magnitude of extreme yield losses. When instead Y represents a quantity where lower values are associated with poorer outcomes, such as absolute production levels, the definition is typically adapted to focus on:

In this case, the Expected Shortfall captures the average production among the lowest-performing outcomes that occur with probability at least 1α, and thus quantifies the expected severity of low-production scenarios. We denote this as the (1α) * 100% ES (for example, the 1%

. 1965 1968 1971 1974 1977 1980 1983 1986 1989 1992 1995 1998 2001 2004 2007 2010 2013 2016 In-Sample Year  1965 1968 1971 1974 1977 1980 1983 1986 1989 1992 1995 1998 2001 2004 2007 2010 2013 2016 In-Sample Year Production (tons) Figures 2 and 3 illustrate the in-sample and out-of-sample predictive intervals for wheat yield log returns in Australia and for global wheat production, respectively. Across both figures, the predictive intervals produced by the HB-MAR-X model (dashed blue lines) consistently demonstrate better empirical coverage than those of the B-MAR-X model (dashed black lines). During the in-sample period, the empirical coverage rates of the 95% prediction intervals for Australian wheat yield are approximately 91% for HB-MAR-X and 72% for B-MAR-X, while for global wheat production, the corresponding coverage rates are about 94% for HB-MAR-X and 100% for B-MAR-X. These results indicate that B-MAR-X tends to underestimate uncertainty for Australian yield returns, while it overestimates uncertainty for global production. This interpretation assumes the superiority of the HB-MAR-X model, which is supported by the higher estimated log evidence reported in Table 1. Notably, the most severe negative logreturn years for Australian wheat yields-1967, 1982, 2003, and 2007-coincide with historically documented, severe drought events. These years correspond to most of the observations falling outside the 95% prediction interval under HB-MAR-X, which is both expected and desirable from a statistical standpoint. A well-calibrated 95% interval should exclude approximately 5% of points, and the fact that HB-MAR-X's misses align closely with known climate-induced yield shocks demonstrates that its predictive uncertainty is appropriately quantified. This alignment confirms that HB-MAR-X captures the underlying yield dynamics effectively. Furthermore, the HB-MAR-X model enhances tail-risk assessment through its 1% ES measure (solid blue lines), explicitly quantifying the severity of potential extreme losses. Compared to B-MAR-X, the HB-MAR-X ES consistently suggests deeper potential downside risks. Importantly, the 1% ES estimates from HB-MAR-X better reflect the severity of potential losses during drought-induced crises such as those in 2003 and 2007. This responsiveness allows stakeholders, including agricultural insurers, policymakers, and farmers, to better anticipate and financially prepare for severe yield disruptions due to climatic volatility. Thus, HB-MAR-X's explicit quantification of extreme downside risk significantly improves agricultural risk management strategies, particularly in enhancing resilience against future climate extremes. In the out-of-sample evaluation over the 21-year forecasting horizon, the HB-MAR-X model continues to outperform the B-MAR-X model in terms of predictive reliability. For Australian wheat yield log returns, HB-MAR-X achieves an empirical coverage rate of approximately 90% compared to 71% for B-MAR-X, again suggesting underestimation of uncertainty by B-MAR-X. Similarly, for global wheat production, HB-MAR-X maintains a more calibrated coverage, while B-MAR-X continues to produce overly wide intervals with full (100%) empirical coverage, implying an overestimation of predictive uncertainty. This contrast is further illustrated by the behavior of the 1% Expected Shortfall (ES) estimates. In both the in-sample and out-of-sample periods, the ES band from B-MAR-X, represented by the solid black line, lies consistently below that of HB-MAR-X, indicating a more pessimistic view of potential extreme losses. This tendency suggests that B-MAR-X systematically overstates downside risk, which may lead to unnecessarily conservative decision-making and inflated reserve allocations. In contrast, the ES estimates from HB-MAR-X remain more restrained in most years but respond appropriately to genuine global production shocks, such as the sharp decline in 2003. The HB-MAR-X model's ability to respond to actual production shocks without overstating risk makes it particularly valuable for managing uncertainty in practical settings. For global stakeholders such as commodity traders, food security agencies, and humanitarian organizations, this translates into clearer signals for when to adjust reserves or implement intervention strategies. Rather than prompting excessive caution in relatively stable years, HB-MAR-X offers a more balanced assessment of risk. This helps decision-makers allocate resources more efficiently while remaining prepared for genuine disruptions driven by climate variability. To comprehensively project future crop yields, the proposed Bayesian hierarchical model integrates climate-driven uncertainties by incorporating projections under representative concentration pathways (RCP2.6, RCP4.5, RCP8.5) and an additional ensemble of 100 carbon concentration pathways generated via QuickClim (Kitsios et al., 2023). While the uncertainties associated with different climate models and emission scenarios have been previously investigated (Li et al., 2025), the uncertainty quantification specifically related to the crop modeling framework itself remains unexplored.
Our study fills this critical research gap by explicitly quantifying uncertainties arising from the crop model parameters and their interactions with the uncertainty in the climate projections. The resulting ensemble of predictive distributions allows stakeholders to comprehensively evaluate uncertainties stemming from climate models, emission scenarios, and crop model parameterizations simultaneously. Consequently, this novel approach delivers robust and transparent probabilistic forecasts, substantially advancing beyond traditional deterministic projections (further details are outlined in Appendix 3). Moreover, incorporating tail risk measures within our Bayesian hierarchical framework improves risk-informed decision-making by explicitly characterizing extreme yield outcomes, thereby equipping stakeholders with probabilistic projections that reflect a broader range of potential future conditions. This comprehensive uncertainty analysis offers decision-makers a powerful tool for anticipating and managing extreme agricultural risks amid future climate variability and technological change. To quantify the uncertainty attributable to differences among climate models, we first computed the mean projected yield for each model at each time point based on samples from the posterior predictive distribution. We then calculated the 25th and 75th percentiles of these

This study presents a new Bayesian hierarchical framework, called HB-MAR-X, designed to improve the way we forecast global crop yields, with a focus on wheat. By explicitly modelling country-specific variances and integrating uncertainties from climate models, emission scenarios, and crop model parameters, our approach addresses key limitations of traditional pooled-error models. The HB-MAR-X framework demonstrates improved calibration and probabilistic accuracy for both in-sample and out-of-sample predictions, particularly during periods marked by extreme yield shocks, such as major drought events. This enhances the reliability of tail risk metrics, such as Value-at-Risk (VaR) and Expected Shortfall (ES), thereby supporting more effective and timely risk management and policy responses. A key strength of the proposed framework is its transparent decomposition of predictive uncertainty into contributions from climate models, emissions pathways, and crop model parameters. Our results show that, while climate model uncertainty remains the dominant source of forecast spread, crop model parameter uncertainty makes a substantial and persistent contribution to the total uncertainty throughout the projection period. This is evident in the marked widening of predictive intervals when crop model parameter uncertainty is included, particularly toward the end of the century. Consequently, comprehensive risk assessment and robust decision-making require that both sources of uncertainty are explicitly quantified and communicated. Furthermore, projections across a wide range of emission pathways reveal that higher emissions are associated with lower ensemble mean yields, broader predictive intervals, and a deepening of downside risk, as measured by the 1% ES. This highlights the increasing vulnerability of global food production to adverse extremes under high-emission scenarios and reinforces the need for adaptation and mitigation strategies that explicitly consider tail risks, not just central estimates. In addition to tail-specific metrics, our Bayesian framework supports a full evaluation of expected climate risk under each emissions scenario s, following recent advances in the literature (Villalta et al., 2020). The unconditional expected loss can be expressed as  A commonly used approach is likelihood annealing, where the likelihood function is gradually introduced through a sequence of tempering parameters γ t , evolving from γ 0 = 0 to γ T = 1. Specifically, at iteration t, the intermediate distribution takes the form

thus allowing the posterior distribution to be approximated progressively. At each iteration t, we reweight particles according to the ratio of target densities, possibly resample to discard low-weight particles, and then apply an MCMC move step to rejuvenate diversity in the particle set. This ensures a robust representation of each intermediate distribution. After resampling, an MCMC move step is applied to each particle to mitigate degeneracy and enhance exploration of the target distribution (Chopin, 2002). In this step, an MCMC kernel proposes a new particle, θ * t , for each current particle, θ i t , using a proposal distribution q(θ * t |θ i t ). We employ a multivariate normal random walk kernel:

where Σt is the empirical covariance matrix of the particle set at iteration t, and h is a scaling factor determined adaptively. The proposed move is accepted with probability , ensuring that the updated particle set remains consistent with the intermediate target distribution π t (θ|y). To optimize particle movement, we tune the scale h based on the expected square jumping distance (ESJD; Sherlock and Roberts, 2009;Salomone et al., 2025), defined as

where θ i t -θ * t Σt denotes the Mahalanobis distance between the current and proposed particles. The optimal scale is selected by comparing the median ESJD across particles for a range of candidate scales, and choosing the scale that maximizes this median. Once the optimal scale is determined, the MCMC kernel is applied for adaptively chosen number of iterations, R t , until the ESJD of a specified proportion of particles exceeds a predetermined threshold. This adaptive move step ensures robust local exploration while maintaining the integrity of the intermediate distribution. In our hierarchical crop yield setting, SMC offers several key benefits. It provides robust exploration by navigating the potentially high-dimensional posterior space of crop yield parameters (including country-level variances) more effectively than a single-stage sampler that attempts to sample all parameters simultaneously. Additionally, by gradually updating from a simpler distribution to the posterior, the algorithm achieves adaptive accuracy, reducing the risk of particle impoverishment and local trapping. Finally, SMC naturally facilitates Bayesian model selection and predictive distributions through its evolving particle set, thereby offering a unified framework for posterior inference and out-of-sample yield prediction. Algorithm 1 describes the key steps of our SMC procedure in detail. This method underpins the Bayesian estimation of our proposed hierarchical crop yield model, and it also facilitates robust model selection and prediction, as detailed in subsequent sections. SMC's efficiency and advantages in various econometric contexts have been noted elsewhere (Li et al., 2021(Li et al., , 2023a)), making it a suitable choice for our application. dence, which indicates improved predictive accuracy relative to the complexity introduced, it is retained in the model.
Conversely, if the addition of a country results in a reduction in model evidence, suggesting an adverse impact due to increased complexity or noise, it is excluded from the final specification. This approach strategically departs from existing modeling practices, where typically all available countries are simultaneously incorporated, potentially introducing unnecessary noise and complexity. In contrast to traditional methodologies employing non-linear least squares point estimators, our proposed selection procedure leverages Bayesian model evidence computed through SMC, thereby providing substantial methodological advantages by inherently accounting for parameter uncertainty and effectively penalizing model complexity. MAR-X model, which incorporates a hierarchical structure for country-level error terms only, we also examine several alternative specifications. These include: (i) a model with hierarchical structures applied solely to the country-level time-varying intercepts, specifically the parameters a and λ; (ii) a model incorporating hierarchical structures for both the country-level intercepts and error terms; and (iii) a fully hierarchical specification that includes hierarchical structures

E[Loss|s] = h v e v p(v|h) p(h|s) dv dh = E[e] E[v|h] p(h|s) dh, where p(h|s)