The risk and consequences of multiple breadbasket failures: an integrated copula and multilayer agent-based modeling approach

Climate shocks to food systems have been thoroughly researched in terms of food security and supply chain management. However, sparse research exists on the dependent nature of climate shocks on food-producing breadbasket regions and their subsequent cascading impacts. In this paper, we propose that a copula approach, combined with a multilayer network and an agent-based model, can give important insights on how tail-dependent shocks can impact food systems. We show how such shocks can potentially cascade within a region through the behavioral interactions of various layers. Based on our suggested framework, we set up a model for India and show that risks due to drought events multiply if tail dependencies during extremes drought is explicitly taken into account. We further demonstrate that the risk is exacerbated if displacement also takes place. In order to quantify the spatial–temporal evolution of climate risks, we introduce a new measure of multilayer vulnerability that we term Vulnerability Rank or VRank. We find that with higher food production losses, the number of agents that are affected increases nonlinearly due to cascading effects in different network layers. These effects spread to the unaffected regions via large-scale displacement causing sudden changes in production, employment and consumption decisions. Thus, demand shifts also force supply-side adjustments of food networks in the months following the climate shock. We suggest that our framework can provide a more accurate picture of food security-related systemic risks caused by multiple breadbasket failures which, in turn, can better inform risk management and humanitarian aid strategies. This paper was partly funded by the IIASA-IDMC-OFDA Grant No. AID-OFDA-G-17-00285 and by the Austrian Climate Research Program (ACRP) Grant No. KR15AC8K12597. * Asjad Naqvi naqvi@iiasa.ac.at Extended author information available on the last page of the article


Introduction
Climate-related shocks to food systems have recently emerged as a major cause for concern (Mehrabi 2020; Gaupp et al. 2020). Such shocks affect not only agriculturedependent regions but non-agrarian regions as well through demand shocks caused by population displacement and supply shocks caused by value chain adjustments (UNDRR 2019). Highly fertile agricultural regions that produce wheat, rice or other grains are typically referred to as "breadbasket" regions. For example, Punjab region in India and Pakistan, fed by rivers in the Indo-Gangetic plain, is responsible for most of the grain productions used for consumption in South Asia, some of which is also exported. Such regions are highly prone to direct climate shocks (FAO 2018), especially droughts that can vary in degrees of magnitude and spatial coverage (CRED 2018;UNDRR 2019). Breadbasket "failures" can also occur across multiple regions as a result of weather-related inter-dependent risks. In 2010, droughts in Russia and floods in Pakistan that were physically connected through atmospheric blocking (Lau and Kim 2012) directly impacted breadbasket regions in these two countries and affected global food prices (Katsafados et al. 2014). This phenomenon, known as "Multiple Breadbasket Failures" or MBBF, has recently been highlighted as a global security concern by climate experts, agriculture economics and food security experts (Janetos et al. 2017;Gaupp et al. 2020), and with climate change, is expected to get worse in the coming decades (Gaupp et al. 2019).
Climate shocks to breadbasket regions can have several direct consequences, such as crop failure, famine and food insecurity (UNDRR 2019). They can also result in indirect outcomes like internal displacement, economic migration, disruption of food supply chains and volatility in food markets (Haraguchi and Lall 2015). For example, extended drought in the past years in the Horn of Africa caused a sever decrease in crop production, food security and job opportunities which have forced many pastoralists and rural farmers to move to urban areas in search of alternative livelihoods (Bonneau 2013). More recently, Philippines, China and India experienced recently large-scale internal displacement due to climate-related disaster events (IDMC 2018). These results are the outcome of demand-side responses driven by household behavior. For example, households in breadbasket regions affected by climate shocks are likely to lose both the ability to produce food and earn income. This can trigger different coping mechanism responses like running down savings, informal borrowing, or, with a very high likelihood, climate-induced displacement to other regions (Auffret 2003). Displacement also impacts other regions by putting a downward pressure on wages and upward pressure on prices (Naqvi 2017). If markets are not agile enough to respond to sudden demographic changes, they can cause severe fluctuations in supply chain adjustments resulting in price volatility.
In the aftermath of such climate shocks, humanitarian aid organizations and governments need to take policy response measures in a short time span, often with 2 Literature review In this section, we give a brief overview of state-of-the-art modeling approaches that deal with climate shocks. We also discuss the contributions of copulas, multilayer networks and agent-based models in this context. Okuyama (2007) highlights three aspects of climate shocks that models need to incorporate: time, space and nonlinear feedback effects. In the last few decades, several waves of modeling efforts have taken place that led to considerable methodological developments in the field of natural disaster analysis in that regard. This includes the early input-output (I-O) models which conducted extensive research on lifeline losses and disruption of transportation networks but focused exclusively on high-income countries like the USA and Japan (Dacy and Kunreuther 1969;Wilson 1982;Rose and Miernyk 1989;Rose et al. 1997). These models were later extended toward more developing countries and emerging economies (Okuyama 2004(Okuyama , 2011. Better availability of data also allowed I-O models to include regional impacts (Okuyama and Santos 2014;Cavallo et al. 2014). However, these models have been criticized for their strong assumptions of linear and fixed relationships across parameters and purely supply-side-driven adjustments.

Economic modeling approaches
Therefore, in addition to I-O models, computable general equilibrium (CGE) and dynamic stochastic general equilibrium (DSGE) models were developed and are now the most commonly used framework for analyzing post-shock outcomes (Cole 1995(Cole , 1998(Cole , 2004. Both CGE and DSGE models focus on long-run equilibrium using rational optimizing agents. Models in this field include national-level (Ueda et al. 2001;Rose and Guha 2004;Rose and Liao 2005) and regional-level studies (Tsuchiya et al. 2007;Hallegatte and Dumas 2009). Both CGE and DSGE models assume perfectly functioning markets, perfect foresight and smooth transitions to a new equilibrium after shocks. While they can now handle regional analysis both due to availability of better computational power and advances in numerical solutions (Ishiwata and Yokomatsu 2018), they still suffer from certain drawbacks. For example, strong assumptions of fixed and exogenously defined parameters remain unchanged during and after the shock and the exclusive focus on long-run equilibrium with little or no discussion on the transition processes. Additionally, the mathematical structure of these models makes it extremely challenging to include a large number of sectors and regions, or other nonlinear behavioral aspects, for example endogenous parameter changes (Rose 2004;Okuyama 2007).
Some recent models have taken non-mainstream approaches. For example, Hallegatte and Ghil (2008) and Hallegatte et al. (2007) introduce business cycles in a mixed demand and supply-driven framework. Albala-Bertrand 1993 shows that while impacts might be negligible at the macro level, they can completely disrupt regional economies. Albala-Bertrand 1993 also highlights the role of labor in estimating overall losses, an aspect that is missing in mainstream models which are 1 3 The risk and consequences of multiple breadbasket failures mostly driven by capital stock accumulation through investments (Skoufias 2003;Hallegatte and Przyluski 2010).
While the above methods are relevant from a long-term policy planning perspective where economic trends are assumed to stabilize and behave "normally," shortrun outcomes, which are relevant for immediate policy response under an uncertain environment need alternative modeling tools (Skoufias 2003;Naqvi and Rehm 2014b). The above literature points out to several gaps, especially two worth noting. First, little has been done to include probabilistic climate risks in macroeconomic models with only some recent literature touching upon this topic (Hochrainer-Stigler et al. 2011;Ishiwata and Yokomatsu 2018;Poledna et al. 2018). Second, the evolution of short-run direct and indirect losses is not clearly understood with competing and conflicting views in the empirical literature (Kahn 2005;Loayza et al. 2012;Horwich 2000;Strobl 2012;Skidmore and Toya 2002;Hallegatte et al. 2007;Crespo et al. 2008;Hallegatte and Dumas 2009).
We discuss in the following subsections on how copulas and agent-based models can be used to model the short-run outcomes following climate shocks to breadbasket regions.

Copulas
While the above-mentioned models focus on economic relationships and the impact of climate shocks, copulas are especially useful to address issues of tail dependency usually not looked at in such models. For example, a copula can explicitly address the probability of a large loss in one region given another region faces a large loss (Sklar 1959;Joe 1997;Nelsen 2007). Initially developed for the financial sector (Aas 2004;Dißmann et al. 2013), copulas are now increasingly applied to the field of natural hazard and disaster risk modeling. For example, the phenomena of increasing tail dependency are noticeable in hydrological processes such as water discharge levels across basins. In the work of Jongman and colleagues (Jongman et al. 2014), it was shown that large-scale atmospheric processes can result in strongly correlated extreme discharges across river basins leading to extreme flooding across large regions over Europe. Also in the case of drought events, the risk of simultaneous large-scale events that cause drought-induced crop losses in several regions at once are found to exhibit tail dependency (Prudhomme and Genevier 2011;Ratnam et al. 2016;Vicente-Serrano and López-Moreno 2006;Gaupp et al. 2019).
Copulas are able to model nonlinear codependencies of risks as in the case of joint climate extremes. They can also be used to model the dependence structure between multiple variables such as climate indicators and crop yields. In agricultural sciences, copulas are typically used to depict the joint effects of temperature and precipitation on crop yields (Cong and Brady 2012). A common class of multivariate copulas used in this field is the so-called Vine copulas (Aas et al. 2009;Bedford and Cooke 2002;Czado et al. 2013) which consist of a cascade of bivariate conditional and unconditional pair-copulas. For our case study, we use the most advanced approach currently available, namely regular vines, or RVines, which can model complex dependencies even in large dimensions (Dißmann et al. 2013). As we will see further below, copulas are a flexible statistical tool to model the dependencies between risks. Due to Sklar's Theorem (Sklar 1959), the underlying univariate marginal distributions, for example of individual risks, can be modeled separately from the dependence structure (Aas et al. 2009). In this way, droughts in different regions, and with varying degrees of duration, severity and distributions, can still be modeled jointly using copulas (Shiau 2006). For our case, we use logistically de-trended wheat yield data in the Indian breadbasket states to demonstrate the importance of incorporating tail dependencies in drought risk estimations. Subsequent cascading impacts are modeled through an agent-based multilayer network approach as discussed in the next section.

Multilayer networks and agent-based models (ABMs)
While inter-dependent risk can be modeled through copula approaches, economic modeling approaches, as discussed above, do not capture the intrinsic dynamics of shocks within systems, and therefore, agent-based modeling approaches are coming nowadays to the forefront (Elsner et al. 2014;Naqvi and Rehm 2014a). Agentbased models, or ABMs, are a bottom-up methodology that allows agents to make decisions based on past and current outcomes as well as on the behavior of others (Bowles 2006). ABMs have the ability to handle a large set of different agent sets allowing for heterogeneity within each agent set. Within such a framework, the interaction of agents determines meso-and macro-outcomes, which can further feedback into the micro-decision-making of individual agents resulting in nonlinear path-dependent outcomes (Arthur 2006;Elsner et al. 2014;Tesfatsion 2006). Thus, ABMs are well equipped to handle large parameter spaces, nonlinear thresholds, boundary conditions and out-of-equilibrium states thus going beyond state-of-theart modeling tools typically used in disaster analysis (Farmer and Foley 2009;Axtell 2005).
However, it is not enough, as currently done, to look at the economic system as a single-layer entity with heterogeneous agents. We suggest applying multilayer network approaches instead (Kivelä et al. 2014;Battiston et al. 2016). Multilayer networks are represented by unique but interconnected layers. Assuming regions can be represented as nodes in a network structure, these nodes can interact with each other in different ways. For example, goods flow can be represented by trade networks and population flows by migration networks. The nodes interact with and across layers based on behavioral rules usually derived from the economic theory or empirical literature based on the region under study (Alfarano and Milaković 2009;Naqvi and Rehm 2014a;Naqvi 2017). This network structure allows us to track changes resulting from climate and economic shocks that might occur in one part of the network, but can cascade throughout the system through behavioral interactions. The specific network structure is typically determined by the number of locations or nodes, the number of layers, endowment of stocks and the strength of connections across the layers.
The multilayer approach was first applied to financial networks after the 2008 global crisis to understand how shocks cascade and can cause "systemic risk" (Bardoscia et al. 2015;Battiston et al. 2012;Poledna et al. 2015). Unlike financial networks, where all transactions are in monetary units across different types of financial institutions or assets (Hurd and Gleeson 2011;Stolbova et al. 2018), food production involves networks comprising of different monetary and non-monetary units. For example, land is used to produce a physical output of crops, which is later converted into monetary units through market transactions. Similarly, agriculture labor, which is typically a large proportion of the total population in developing countries, engages in farm work in exchange for usually low wages. A climate shock like flood or drought to such a system implies that food output is lost through less land being available. A more catastrophic shock like earthquakes can even cause high casualties reducing labor power in a region. Thus, losses in different layers cannot be treated in the same way as losses within financial networks. Additionally, unlike financial networks which are governed with relatively well-defined homogeneous rules and regulations, the agriculture-food web and its corresponding agrarian and farming population need to be set up with layer-specific behavioral rules. In order to do so, the dynamics within each dimension need to be modeled through the behavioral interaction of different agent sets belonging to different layers across a network of spatially explicit locations.
In the next section, we describe how the three methods of copulas, multilayer networks and agent-based models can be combined together to tackle the challenges addressed above, that is, how to model simultaneous risks and cascading impacts across regions.

Modeling climate shocks through copulas
Copula approaches gained their importance due to Sklar's Theorem (Sklar 1959). It states that every multivariate distribution F with margins F 1 , … , F d can be written as: with C as a d-dimensional copula, a multivariate distribution on the unit hypercube where c is the copula density, f, the multi-dimensional probability density function of F, and f k , the marginal probability density function of F k , k = 1, … , d . d in our example refers to the number of breadbasket states in India. The advantage of this approach is that the copula, and thereby the dependence structure, can be chosen independently from the marginal distributions. There are different copula families capturing different dependence structures (Jaworski et al. 2010) but generally speaking, any multivariate distribution can be used as the basis for a copula (Aas 2007). To structure our multivariate copula model, we used regular vines (RVines) (Aas et al. 2009;Bedford and Cooke 2002;Kurowicka and Cooke 2006;Dißmann et al. 2013), a flexible graphical model that consists of a cascade (or tree structure) of conditional and unconditional bivariate copulas that decompose the multivariate probability density. A d-dimensional RVine model consists of d − 1 trees. Each tree consists of nodes and edges joining the nodes. The first tree identifies d − 1 pairs of variables with a directly modeled distribution. The second tree identifies d − 2 pairs with a distribution conditional on a single variable, modeled by a pair-copula. Continuing in that way, the last tree consists of a single pair-copula with a distribution conditional on all remaining variables. For a graphical representation of a RVine tree, see (Dißmann et al. 2013) or (Brechmann and Schepsmeier 2013). The RVine tree structure in this study is selected using the maximum spanning tree approach with Kendall's as edge weights (Dißmann et al. 2013): with ̂m ,n as pairwise empirical Kendall's (that is, a rank correlation coefficient, see Embrechts et al. 1997). As will be discussed in the results section, this approach was adopted for India on the state level, which showed large tail dependencies for some regions and therefore is an excellent example for showing the advantages of a copula approach when it comes to determine systemic risks.

The multilayer networks and agent-based models
Going to the ABM modeling component, we first want to note that every ABM needs to be contextualized, that is, it must be set up for a specific purpose of analysis (Gilbert 2008;Klabunde and Willekens 2016;Naqvi 2017), which in our case are drought-related production shocks in India. A shock in the food system through a strong reduction in crop yields will cascade through the affected population as well as through the unaffected one. For example, in the case of India, a state's economy can be assumed to produce output across different sectors that can be divided into agricultural or other goods (physical goods and services). Each state holds a stock of workers which are employed across the different sectors that interact with each other to form a simple circular flow economy. Cross-state interactions result in the exchange of goods (including food) and services. Underlying this production and trade interactions are labor markets, which determine how much labor is employed across different sectors. Like trade which is driven by price and profit signals, income differentials result in migration across locations. This, in turn, determines population and income distributions. Migration decisions are also influenced by social signals, for example social-or community-based networks. The household The risk and consequences of multiple breadbasket failures and the production layer combined represents a multilayer network structure as shown in Fig. 1.
Each dot in Fig. 1 represents a location divided across two layers. A production layer employs workers to produce output. This output is in demand not only in one location but in other locations as well. The second layer represents the demand side, where households contribute via two channels: by providing labor for production in exchange for income and demand for consumption purchased through income and savings. Labor can also move around based on market signals or when faced with high-distress scenarios as in the case of shocks.

A multilayer network
In order to operationalize this framework for India, we set up a multilayer network model, where the 36 Indian states are represented as nodes and two layers shown in Fig. 1. The first layer represents low-income households, and a goods layer represents the production of food and other goods (see Fig. 1). Each node is calibrated using actual data to determine the level of production of food and other goods in each state (see Table 1 in "Appendix A").
The nodes which produce food are subjected to a supply-side breadbasket failure shock, which means simultaneous production losses in multiple Indian states. In our case, two specific shock scenarios are tested, (a) a shock assuming independence between food production across states and (b) a copula-based shock derived from a probability distributions that includes production inter-dependencies across the states. In the copula scenario, the number of agents affected simultaneously increases with the intensity of the production shock. This impacts other layers as well through employment, income and consumption decisions. For both shocks, we simulate a no-migration and a migration scenario to test the impact of the shock Fig. 1 A multi-dimensional network layer across the region. The no-migration scenario represents a single-layer simulation where only the production side (or the supply-side) is simulated.
Summarizing, our approach represents two improvements the over existing literature. First, shocks are not assumed to be independent, uncorrelated events, but follow a copula-based structure. Second, food production shocks impact other layers, notably the household layer, resulting in a transmission of the shock to non-affected regions. The ABM model itself is discussed in more detail next.
It should be noted that for the sake of simplicity, the model assumes that all layers adjust completely to demand and supply interactions. Mismatch in markets is corrected through prices and all output, and income is fully exhausted such that there are no leftover stocks. This mimics a general equilibrium framework where prices drive the model outcomes. In a more realistic setting, path dependencies and institutional barriers will also very likely play a role as well. Thus, a full adjustment model is a best-case scenario model and distributions resulting from the simulations shown here therefore represent a lower bound for the level of disparities created from the climate shocks.

The agent-based model
In order to operationalize the model, we need to endow the nodes and links with behavioral rules. These are introduced in the model in a very basic form using standard economic assumptions described below.
Each node i is endowed with land and capital stock that allows for a maximum production of food (F) and other consumable goods (G), respectively. The production level of each node is derived from Table 1. Assuming j is the index of goods produces at each location i, which in our model is j = F, G , the value of the total output produced at a node i equals This two-good model framework can be easily extended to include a larger set of goods. To produce these goods, all the labor supply L i available at each node is employed. Assuming each unit is produced at a cost , the total wage bill is simply y i . The wage rate ( w i ) of workers at a node i equals the total wage bill over the total labor supply, or The income earned by workers is fully spent on the goods produced at the node where a fraction is spent on food F and 1 − on other goods G. Or more formally, the consumption of each good at a location i is defined as: The risk and consequences of multiple breadbasket failures From this, the price of the two goods can be derived as: For the sake of simplicity, the average price at location i, p i , is taken as a simple arithmetic average of the two prices. This can be made more complex to represent goods baskets with price weights to generate average price indices at different locations.
In nominal terms, the value of the total output produced at a node i can be derived as If the nodes are allowed to interact, then a gravity model-like specification determines how much is exchanged across the nodes (Anderson 1979). This is estimated in the model as a logistic function Π iq between node i and its q connected neighbors following a joint probability distribution of the type: In short, the probability of moving to a connected neighbor q or Π iq is positively affected by relative economic gains, Π h iq , at destination q. For migration, the economic gains are determined by real-income differences ( h = w q ∕w i ) and for trade by potential profit gains defined by relative prices ( h = p q ∕p i ). The choice of destination q is also negatively affected by distance Π r iq to neighbor q. The generic logistic function is of the type: where a and b are shape parameters of the logistic function. Additional migration-decision factors can also be added, for example community at destination and other social linkages (see (Naqvi 2017) for the full model specification).
Two thresholds are introduced in the model to make the post-shock decisionmaking behavior more realistic. First, households adjust their propensity to consume food relative to a minimum consumption value c . If their income goes down or food prices go up, then will adjust endogenously to ensure that they stay at least at the c level. Formally, Since households cannot spend more than what they earn, is bounded such that ̄≤ ≤ 1 . Like households, firms only sell in markets which give them positive profits. Unlike households, they can diversify their portfolio to sell in several markets. This also allows them to dynamically adjust their supply to changes in demand in their vicinity. Since firms have production costs (wages in this model), the difference between the costs and the market prices determines the profit margins. Since the unit cost of production equals and the market price equals p j i , then the condition p j i ≥ defines the sorting for selling in markets. If some markets fall below this threshold, or p j i < , then they are excluded from the list. These two thresholds evolve dynamically since all variables are time indexed.

A multilayer vulnerability index
Following the advances in multi-dimensional network measures (Magnani et al. 2013;Kivelä et al. 2014;Gemmetto et al. 2016), a similar measure is adapted here to estimate multilayer vulnerability. Vulnerability is simply defined as the populations falling below the minimum consumption threshold ( c ). In spirit of DebtRank, a key centrality measure of multilayer systemic risk in financial networks (Bardoscia et al. 2015;Battiston et al. 2012;Thurner and Poledna 2013) that itself is shaped after Google's PageRank (Brin and Page 1998;Halu et al. 2013), we develop a new vulnerability measure to represent multilayer food insecurity, a key issue in disasterprone regions (FAO 2015). We will call this measure Vulnerability Rank or VRank, which is estimated as follows: In other words, the vulnerability of a node i is not only defined as the purchasing power of food of populations at that node, but also all its neighboring q nodes. The parameter is the dampening factor, typically set at = 0.8 , which determines the impact of the neighbors on node i. This formula captures two key aspects. First, it encapsulates the dependence of one node on another in terms of food supply. For example, a node which does not produce any food would be considered vulnerable in isolation, but if it is connected by several food-producing regions that can service it in times of shocks, it should have a much lower vulnerability than a food-producing node that is completely isolated. Second, it captures the displacement-based linkages as well. A node that produces enough food for its current population might see a sudden influx of populations from shocked neighbors, making it potentially very vulnerable in the short-run following the climate shock.
1 3 The risk and consequences of multiple breadbasket failures

Results
In this section, we now integrate our suggested modeling approach as discussed above to the case of India. The country serves very well for our purpose. Around 55% of the Indian workforce lists agriculture and associated activities as their primary source of employment. Agriculture also contributed 17.4% to India's real gross values added (GVA) in 2016-17 (Department of Agriculture India 2018). In rural areas, agriculture still generates more than half of the total income (Kadiyala et al. 2014). After rice, wheat is the most important cereal in India (40%) with cereals accounting for 94% of total food grains (Department of Agriculture India 2018). Only around 50% of sown area are irrigated which means that agriculture is still largely dependent on monsoon rainfall which occurs, depending on the region, between June and September. Additionally, in dry years low reservoir levels and fuel shortages can hamper irrigation efforts (Krishna Kumar et al. 2004). Winter crops such as wheat benefit from the monsoon rainfall through residual soil moisture. If the monsoon fails, India experiences risks of crop failure in wide parts of the country.

Risk assessment of a MBBF in India
To give an understanding of the underlying principle of dependent crop losses, we start with an example of two dependent Indian states, namely Uttar Pradesh and Rajasthan, both part of the Indian wheat breadbasket. Figure 2 shows the univariate marginal distributions of logistically de-trended (Gaupp et al. 2017) wheat yields in Uttar Pradesh and Rajasthan. Both de-trended yields follow a normal distribution tested with the Shapiro-Wilk test.
In the middle, their scatter plot is shown as well as the contour plot of the Frank copula which was selected as the bestfitting copula type. It joins the marginals of de-trended wheat yields in Uttar Pradesh and Rajasthan and has the following form: with as the copula parameter determining the strength of tail dependence. To stress the importance of dependence, we take a look at the 1974 food crisis in India. A combination of adverse weather conditions and a rise in oil prices with consequent shortages of petroleum-based fertilizer led to unexpectedly low yields in India (Weinraub 1974a, b;Joerin and Joerin 2013). Based on our data, the likelihood of yields as low as the 1974 yields occurring simultaneously in Uttar Pradesh and Rajasthan equals C [F UP (−0.29), F R (−0.16)] = 0.027 . If we would not take into account the dependence structure between the two states, the likelihood of both states experiencing a 1974-type event would be F UP (−0.29), F R (−0.16) = 0.007 . To estimate the joint probability of all eight breadbasket states, we use an RVine tree structure as described in the methods section. Figure 3 shows the first tree of the The risk and consequences of multiple breadbasket failures de-trended wheat yields in the Indian breadbasket. Edges are labeled with the copula families connecting the nodes via bivariate copulas. The full tree structures which include all conditional copulas are shown in Fig. 7 in "Appendix B". In 2003, following the 2002 monsoon failure, India experienced a major wheat yield shock (with shock defined as deviation from the long-term logistic trend). Based on our analysis in eight breadbasket states (see Table 2), the wheat yield shock in 2003 in the Indian breadbasket was a 1-in-17 years event. The importance of not only considering the yield distributions in the different states in an agricultural risk analysis but also the dependence structure between them becomes clear when we compare the return period of the 2003 yield shock with the return period if we would ignore the dependencies. If we assume independence between the states, the likelihood of a 2003 event would decrease to a 1-in-370 years event, a 20-fold increase in the magnitude. Therefore, by ignoring linkages between the production in different states, the risks of crop failure are underestimated dramatically. In addition, further spillover effects of those spatially inter-dependent yield shocks have to be modeled to fully grasp the systemic risks of our current food system.
The above framework is applied to India to estimate the impact of a breadbasket failure on internal displacement. Figure 4b shows the states of India of which eight (shown as green dots) are major food crop producing states.This represents 40-45% of all employment in India, most of which is at extremely low end of the income distribution. The 2003 crop failure resulted in a rise in food prices. Since India is not a net importer of food items, the loss in crops resulted in an internal market adjustment to the climate shock.

Simulating the effects of a MBBF in India
We conduct four experiments to show the usefulness of our approach, using a combination of no-displacement and displacement scenarios and independent and copula-based probability distributions. The displacement scenario allows populations to move from one state to another if they see better real-income opportunities. Due to a lack of state-to-state migration or displacement data, which can allow us to calculate migration propensities, all states are given equal weights in terms of destination choices, where the only negative factor is average distance where farther away locations are considered less attractive than nearby ones. In a more real-world scenario, propensity to move would also be determined by community and language ties, work opportunities and other social and cultural factors (Black et al. 2013;Klabunde and Willekens 2016;Cattaneo and Peri 2016). Thus, our displacement scenarios represent a highly stylized baseline outcome, which, in reality, is likely to be much worse. The 36 states are set up in the simulation model using data from Table 1 for calibration. The table provides information on low-income populations, total output produced in each sector, and the share of agriculture in this output. Low-income 1 3 workers are divided into these sectors accordingly. The simulations run till states achieve stable trends in prices and population distributions (Fig. 4b).
The breadbasket states, shown in green in Fig. 4, are subjected to a production shock to simulate the drought event. As described above, the shocks can either be modeled as independent or follow an inter-dependent copula structure. The probabilities of changes in yields using both independent probability distributions and a copula-based probability distribution for a 100 year production shock event are estimated and provided in Table 1. Figure 5 shows the temporal trends for the four scenarios for four different indicators: real income, average consumption, marginal propensity to consume and the VRank index. Figure 5a shows that the real income, or the price of labor relative to the price of food, falls across all scenarios. The independent, no-displacement scenario shows the smallest change, while the copula with displacement scenario shows the highest impact. This can be explained by the fact that wheat yields are strongly positively correlated between the states which means that a shock in one state leads to a shock in surrounding states with a high probability. A copula-based shock shows on average a larger decline in relative price changes. Since income falls, and food becomes more expensive, due to an overall decline in output, the average consumption levels, shown in Fig. 5b, also fall.

Fig. 5
Temporal trends from the four simulation runs. Note: All graphs are shown as percentage changes from baseline no-shock scenario As consumption levels decline, households allocate more of their income toward food in order to stay above the minimum consumption threshold. Figure 5a, b indicates that including and excluding dependencies between states have a much higher impact than with or without displacement. Figure 5c shows the development of this indicator, which rises steeply if one assumes tail dependence with displacement. This scenario also takes the longest to stabilize. Figure 5d shows the VRank index, which on average increases for all scenarios with the copula scenarios showing the highest increase. This last graph highlights that underestimating both the joint probability risk and the multilayer displacement effects can result in underestimating the full extent of the shocks. Figure 6 shows the results as spatial choropleth maps at the end of the six-month simulation period, when the indicators stabilize. The graphs show changes relative to the baseline no-shock scenario. Displacement scenarios, shown on the right column in Fig. 6, exhibit a relatively higher vulnerability as the The risk and consequences of multiple breadbasket failures shock spreads from the affected states to the non-affected states through migration. In the non-displacement scenario, the copula-based distribution of the shock shows a higher increase in vulnerability.
The results provide important insights in terms of dynamics as well as policy-wise. For example, while an aggregated macroeconomic performance measure such as the GDP exhibits similar patterns across two displacement scenarios they can completely differ in the distributional effects. Furthermore, such different outcomes also imply that completely different policy responses are needed to tackle such risk. As recently suggested by Pflug and Pichler (2018), the difference in distributional impacts across the dependent and independent case can be used for determining the share of systemic risk (due to the dependency between states) to overall risk. If we conducted a complete probabilistic assessment of all possible agriculture-related climate risk scenarios, which in principle is possible as both the copula and marginal distributions are available, it would be possible to apply classic risk management strategies as well as to analyze humanitarian aid responses. However, this would be extremely difficult to operationalize in short time span due to the high computational requirements (for example, thousand of risk distributions need to be generated and analyzed via the ABM model). Alternatively, so-called risk-layer approaches which focus on risk reduction for more frequent risks and risk financing for infrequent ones  can provide guidance on best strategies forward. As we were looking at independence and 100-year event cases, from a risk perspective, these events would still belong to the "relatively frequent" set. Therefore, risk reduction measures such as crop storage or subsidies during time of crisis could be, in principle, applied and studied with our approach. However, there are limitations due to data scarcity, especially as displacement patterns on such large scales are difficult to be empirically tested and validated. Nevertheless, current efforts are underway to get a more complete picture (Kc et al. 2018).
While these results are derived from a highly stylized model, where adjustments in markets are perfect, a fully calibrated model would include more frictions where the results would show much worse outcomes and distributions. Thus, our model results can be taken as an upper bound for disaster impacts in a best-case scenario.
Model calibration and application to countries are always a challenge for agent-based models. A model on a subregional level, for example provinces or states, has relatively abundant data, allowing for calibration since such regions typically represent standard administrative units (like provinces or districts) at which data are collected. Thus, a subregional level allows the model to be flexible enough to be adopted to different countries. The model can be scaled further down, for example to the village or even household level, but data limitations and the model scale would make the results both highly intractable and potentially subject to extensive sensitivity testing (Naqvi 2017).

Conclusions and directions for future research
Climate shocks result in humanitarian crises that usually affect locations beyond the areas of direct impact. Understanding and responding to these events, given limited resources, are critical for policymakers and aid organizations that are typically making decisions under uncertainty and in a short span of time. State-of-the-art modeling tools usually deal with long-run outcomes and are not adequate to explain short-run nonlinear dynamics that emerge from complex interactions following a climate shock scenario.
To close this research gap, in this paper we proposed that risks associated with climate shocks and their subsequent outcomes can be modeled by combining copulas with multilayer networks and agent-based models. Copulas provide the dependent structures of climate risks across multiple regions, multilayer networks describe how these regions are interconnected, and agent-based models allow us to embed behavioral rules that define socioeconomic interactions. While most of the advances in the above three methods evolved mainly in the fields of banking and finance, they have strong applications in the domain of humanitarian operations.
We discussed how such a model can be constructed and applied to drought shocks in agriculture-dependent breadbasket regions in India where weather-related interlinked dependencies across various Indian states can result in multiple breadbasket failures (MBBFs). Through the modeling exercise, we explored how these shocks can spread across other regions. Using the copula approach, we showed that with higher food production losses, the number of agents that were affected increased nonlinearly which led to cascading spillover effects in other locations via network layers through employment, income and consumption decisions. The simulation results showed large-scale displacement which led to large-scale demand and supply-side adjustments in the months following the climate shock. We summarized the food insecurity arising from interactions across the multilayer economic linkages through a new risk measure, Vulnerability Rank or VRank, which showed that risk increases significantly if inter-dependent structure of climate shocks and the resulting displacement are fully accounted for. Such a modeling framework can be easily applied and calibrated to low-income countries using regional-level data that is typically available from statistical offices. While the current model focused on the short-run (that is, the six-month transition phase that can occur after shocks), it can also be extended to include long-term effects of disaster financing (Mechler et al. 2008) and demographic shifts resulting from long-term displacement or permanent migration.
The suggested framework can be used for testing and employing traditional, as well as, novel risk management strategies, for example food storage and insurance schemes, respectively. It can also be used to quickly identify highly vulnerable locations which require humanitarian aid. Thus, our suggested framework, if effectively utilized, for example in the case of MBBFs, can help minimize or prevent secondary-level post-shock outcomes like displacement and food shortages.

3
The risk and consequences of multiple breadbasket failures A coupla-based multilayer agent-based model can also be deployed to simulate different policy scenarios, for example price controls, allowing policymakers to be better informed about the positive and negative economic and human impacts, changes in risks over time, as well as how to address such issues in the future. These applications are not just limited to climate shocks to agricultural regions but can also be extended to various other aspects of humanitarian crises, for example water and energy issues, early warning systems, potential market failures, spread of diseases, demographic shifts and spending on reconstruction and rehabilitation. We suggest that in all of these domains, where operations research is a key management tool, such modeling innovations can be utilized to manage and minimize various forms of systemic and interconnected risks.