Climatological and hydrological disasters have grown by a staggering 3.5% per year since 1980 (Munich Re, NatCatSERVICE 2018, https://natcatservice.munichre.com/). Such disasters and extremes often distort fundamental physiological processes of plants thus affecting crop yields in terms of quantity and quality [1, 2] . Following large-scale agroclimatic events, such as a heatwave or a drought, lower yields may also progressively give rise to sweeping socioeconomic repercussions upon seed prices, the value of crop production, farm income, and food security, to name a few. The implications of this chain of challenges for global food security have given birth to various recent studies that contemplate detrimental effects of climate extremes on crop yields [3] [4] [5] and agricultural and food systems [6, 7] . In a future scenario where climate change is increasingly colliding with major stressors of the global food system, such as population growth, environmental degradation, and trade interdependence, the risk of food-supply instabilities attributable to more frequent and intense climate extremes is expected to increase [6, [8] [9] [10] . For this reason, international organisations that monitor food systems are showing signs of increasing interest to understand the adaptive capacity of agricultural and food markets to extreme shocks (e.g., Chinese Ministry of Agriculture 2016 www.fao.org/china/news/detail-events/en/c/1191739; OECD/FAO 2018www.fao.org/3/I9166EN/I9166en.pdf).
Food security relies on the resilience of agricultural commodity markets to shocks; that is, on the capacity of supply and demand, with or without human intervention by the market participants, to continuously adjust to disruptions while still producing and providing human food and animal feed. In general, agricultural and food markets have many degrees of freedom to achieve resilience against climate variability and change through various channels that complement the dominant goal of yield improvement and production expansion. Seed traders sell new drought-resistant cultivars every year. Farmers buy those cultivars and also maneuver soil and irrigation management practices. Consumers diversify food patterns and reduce food loss and waste. Countries suffering damaging climate extremes temporarily lift trade barriers, such as import tariffs, to avoid skyrocketing domestic prices [6] . Net importers not only develop a variety of stock schemes to guarantee food security and supply stabilisation at the national level but also formalize international emergency reserves (e.g., APTERR) for humanitarian purposes [11] . Such 'market adjustments' occur irrespective of the origin of the disruption. Recent prominent examples of non-climate-driven disruptions that will likely keep impacting globally integrated food markets are the African Swine Fever and Covid-19 pandemics.
Inspired by recent efforts to quantify the resilience of ecosystems [12] [13] [14] [15] [16] , this article illustrates how a robust and simple indicator of resilience can be calculated taking the example of food production systems. The underlying indicator informs us both about the degree of resistance of food production to damage as well as the degree of recovery. Therefore, it helps to assess and rank the capacity of domestic food markets in particular and society in general to cope with supply-chain disturbances attributable to production, trade and price disruptions of any origin. The indicator is flexible in many ways as it can be calculated (i) on variables that measure any dimension (e.g., production value as a proxy of economic performance, self-sufficiency as a proxy of food security), (ii) using freely available historical time-series, and (iii) on historical or even futureprojected data. For this reason, such a synthesised resilience indicator can be highly relevant for national agricultural policy officials, who can use it as a tool to inform decision-making processes about the evolution of the most -or least-resilient food production systems.
The resilience of food production systems can be calculated consistently with the original definition of resilience, firstly introduced in ecology, as the ability of a system to absorb a shock without losing its function [17, 18] . This can be achieved through the annual production resilience indicator [15, 16] :
where µ is the mean and σ the standard deviation of the production time-series. The value of the indicator is inversely proportional to the risk of annual production losses, and it has been tested for annual agricultural time-series [15, 16, 19] , vegetation primary production and water resources [20] .
Real time-series often display non-stationarities and trends that hinder the computation of Rp because µ and σ are not well defined [16, 19] . In this case, time-series need to be normalized in such a way to filter out the low-frequency variability and to isolate the annual production fluctuations whose amplitude and frequency can be related to resilience. This requires computing the smoothed production time-series using LOESS algorithm [21] , for instance, which is used to normalize the production anomalies and standard deviation:
(2) Pi = loess (pi),
where the pi represents the production values of the time-series under evaluation, πi are the normalized time-series with respect to the baseline values, i.e. the smoothed time-series Pi, and σ' is the standard deviation of the normalized time-series. The non-stationary crop production resilience indicator is simply given by the inverse squared standard deviation of the normalized anomalies [16, 19] :
In case the production time-series is stationary, the baseline is constant and R'C is exactly equal to RC.
In this paper, we present a simple and accessible code that computes production resilience from nonstationary time-series and some examples of analysis that can be conducted by stakeholders, policy makers, and any person interested in estimating the resilience of the agricultural production at the local, national and supra-national levels. This code can be easily adapted to any time-series of positively defined values.
PyResPro is composed of the ProSeries class and some functions that are presented in the following frame: ProSeries class consists of the time-series object creation with a label name and production time-series, which is passed through a pandas DataFrame, and includes several methods:
-The length method returns the number of years of the time-series.
-The mean method returns the production average of the time-series.
-The smooth method returns the smoothed time-series through equation 2. The span parameter of LOESS smoothing corresponds to 20 years, which is hardcoded in this initial version of the software. This value has been calibrated in previous studies [10, 22] .
-The norm method normalizes the production time-series (equation 3).
-The p-res method computes the resilience indicator using equations 4 and 5.
-The __add__ methods allows aggregating two individual time-series. The result is a time-series defined by the sum of the production in the years when both individual time-series are defined.
-The copy method returns a shallow copy of the time-series.
-The acor method computes the correlation coefficient between two individual time-series.
-The plot method produces a plot of the original time-series and the smoothed one, both with the same colour.
The functions included in PyResPro are useful to analyse groups of time-series. Given a list of time-series, the function tot_res returns with all necessary information to characterize the total resilience of a diversified production system are the name, the (average) productions and the resilience of the individual time-series, the resilience of the time-series obtained summing progressively the individual time-series, the pairwise correlations between the time-series, the length of the individual and of the aggregated time-series.
The function res_plot produce a plot of all information computed by tot_res (if moreinfo is set to true) or just the average production and the resiliences of the individual and aggregated time-series (if moreinfo is left false).
The examples illustrated in this section are based of the freely available FAOSTAT data named 'Value of Agricultural Production'. Paraphrasing the FAOSTAT documentation (http://fenixservices.fao.org/faostat/static/documents/QV/QV_e.pdf), the value of gross production has been compiled by multiplying gross production in physical terms by output prices at farm gate. Thus, value of production measures production in monetary terms at the farm gate level. Since intermediate uses within the agricultural sector (seed and feed) have not been subtracted from production data, this value of production aggregate refers to the notion of "gross production".
The unit of measure of the downloaded data is in millions of USD, taking years 2004-2006 as a reference. This implies that joint macroeconomic variability (i.e., the variability of all local exchange rates w.r.t. the USD) is subsumed into the European time series of agricultural production values. However, the smooth function filters out the low frequency variability of this effect. Future analysis could consider the Eurozone only, converting the time series in EUR. For, the present paper, aiming aims at presenting some examples of PyResPro package application that can be easily replicated with publicly available data, no country selection and currency conversion is applied.
Two examples of application of the PyResPro package are given:
the resilience estimation for a single time-series or simple combination of them (Section 3.2) and the resilience estimation for a diversified system production composed of different countries or different crops (Section 3.3).
import numpy as np import pandas as pd import statsmodels.api as sm import csv import matplotlib.pyplot as plt import matplotlib as mpl import itertools as it # FAOSTAT data file file = "FAOSTAT_data_2-25-2020_crops-europe.csv" data = pd.read_csv(file,sep=',') The code reported above provides an example of a basic use of the methods of the ProSeries class for wheat in Italy, which is an important European commodity. The output of the print statement is: "Wheat -France : time series length = 56 , P-res = 97". The plot in Figure 1 allows a quick visual inspection of the time-series and of the anomalies with respect to the non-linear smoothing performed through the LOESS algorithm. Wheat production in France increased rapidly till about the 1990s and then started stagnating [23] . This change of tendency roughly corresponds to a climatic shift bringing warmer spring temperatures in France [24, 25] , which shorten the growing season reducing yields [22] . Additionally, the variability increases as well in the recent decades. Worthy of notice are the effects of the 2003 heat wave [26] and the 2016 yield loss for bad weather conditions [22, 27] .
The same analysis can be replicated for any commodity in any country such as Italy, which is the main European durum wheat producer.
# select another commodity and country sel1 = 'Wheat' sel2 = 'Italy' # same as before but define the time-series as "ts2" # … Figure 2 : As Figure 1 , but for Italy. The resilience value for this time-series is Rp=138.
Wheat production in Italy decreased significantly since the 1970s to 2000s. This negative trend can be attributed to reduction of sown wheat areas [19] . The largest loss happened in 1977. Large negative fluctuations are more frequent in the two most recent decades; an example of such event is year 2003, when heat wave significantly hampered growing conditions during sensitive stages in late spring. However, the fluctuations compared to the baseline values are smaller than in France, in fact the annual production resilience indicator computed for Italy is larger than that for France. increasing production diversity using different crops with varying climate sensitivity to spread the risk of unfavourable climate events could be an effective measure to counteract the increase of extreme events and other negative effects of climate change. The effects of diversification on total crop production resilience can be quantified comparing the resilience of individual time series to that of the sum of production time-series in different countries or different crops. Figure 3 shows an example considering the different climates and production values of wheat in France and Italy, which is characterized by larger resilience than the individual countries. The output of the print statement gives: Wheat-France + Wheat-Italy : time series length = 56 , P-res = 159 , correlation = 0.33. If the two time-series were characterized by the same mean and variance, and uncorrelated between each other, the resilience of their sum is expected to double with respect to the individual ones [15, 16] . The two time-series are indeed only partially correlated, with Pearson correlation coefficient equal to 0.33. The resulting compensation effect is cancelling out some of the fluctuations, despite the fact that average production in France is much larger than in Italy, also because the Italian wheat resilience is larger than the French one. The next section shows more complex analyses involving more than two time-series.
The PyResPro package can be used to characterize the overall production resilience of a single crop over different spatial units [19] , or the resilience of a complex production system composed of different crops [16] . These are potentially more complex analyses that can be easily achieved with the methods of the ProSeries class and the tot_res and res_plot functions. Let's consider, for instance, wheat in the whole Europe. The code in the following frame selects the relevant data for wheat in Europe. It sorts the countries from the larger to the smaller producer, limiting the number a maximum of fifteen producers (for example) in this illustrative example. Only time-series with at least 30 production years are considered in order to achieve a sufficient accuracy in the estimation of resilience [15, 16] . The code presented above computes the resilience of the individual countries and of the progressive aggregation of countries, which is useful to understand the effect of spatial aggregation on the overall resilience of wheat production is Europe. The time-series of the top wheat producers is displayed in Figure  4 . France is the most prominent producer followed by Germany, Italy and the United Kingdom.
The code in the following frame, produces the "resilience-diversity plot" (Figure 5 ), which is useful to understand the relative weight of the different countries in building up the overall resilience of wheat production in Europe.
# call the function computing system resilience totres = tot_res(sel1, tss) # call the function producing the resilience-diversity plot fig = res_plot(totres, moreinfo=True, ylabel=ylabel)
fig.suptitle(sel1 + ' commodity production system resilience') Figure 5 : Individual (black dots) and cumulated (red lines) resilience of wheat for the top ten European countries. Cumulated resilience is obtained summing progressively the production time-series and computing the resilience from the sum of the production. For instance, cumulated resilience labelled as 'Germany' is obtained summing the production of France and Germany. The cumulated resilience labelled as 'Italy' is obtained summing the production of Italy to that of the previous two countries, and so on. The bar plot reports the average production of the different countries. The colours of the bars indicate the pairwise correlation coefficients between the production anomalies each country and the sum of the previous countries in the sorted list. Numbers in the x-labels correspond to the number of valid data in the individual time series and in the aggregated ones.
From the plot on Figure 5 it is possible to estimate the effect of the spatial aggregation of national time-series on the overall resilience of the European production system. The contribution of each country depends on the average production values. The effect of the cross-correlations between the different time series is modulating the effect of diversity on resilience [15, 16] . The production resilience of the top wheat producer (France) greatly increases when Germany and Italy -i.e. the second and third producers -production are summed up. Germany is characterized by the largest individual resilience. Italian production is only partially correlated to the sum of France and Germany, which considerably increases the resilience of the aggregated time-series. The United Kingdom resilience is characterized by low resilience, and it is highly correlated to the sum of the first three countries. Thus, it does not bring positive contribution to the overall resilience. The other countries, characterized by lower individual resilience and lower production, contribute negligibly to the overall wheat production resilience.
The effects of crop diversity for the overall resilience of a single country can be conducted simply by inverting the selection labels (see code block below) and then executing that same code reported in the last three blocks presented above. This simple modification yields the results shown in Figure 6 . # same as previous block # … Figure 6 : Same as Figure 5 , but for the main crops in Italy.
In Italy (Figure 6 ), we note that olive and tomatoes are not providing additional resilience with respect to the primary commodity, which is grapes. This happens because these three crops are quite correlated between them, and olives are characterized by lower resilience. The resilience of the Italian commodity production systems gradually increases when the other commodities are considered, even though they are characterized by lower production. The accuracy of the resilience estimation decreases when Kiwi are considered, because they are represented by shorter time-series.
This paper presents the PyResPro package and provides a sensible methodology to evaluate the resilience of the commodity production system at the annual time-scale. The examples proposed in this paper are based on freely available datasets that are provided in a common structure and format, which also allow the replicability of the presented results.
The same methodology can be applied to crop production time-series evaluated in terms of raw quantity [19] or caloric content [16] , but also to natural ecosystem such as vegetation primary production [20] in order to evaluate the reliability of the related ecosystem services and of water resources. The PyResPro software can be applied at different time-scales as well, such as the monthly or daily time-scales, and in different contexts ranging from natural science to economy and society.