Predicting economic resilience of territories in Italy during the COVID-19 first lockdown

This paper aims to predict the economic resilience to crises of territories based on local pre-existing socioeconomic characteristics. Specifically, we consider the case of Italian municipalities during the first wave of the COVID-19 pandemic, leveraging a large-scale dataset of cardholders performing transactions in Point-of-Sales. Based on a set of machine learning classifiers, we show that network-based measures and variables related to the social, economic, demographic and environmental dimensions are relevant predictors of the economic resilience of Italian municipalities to the crisis. In particular, we find accurate classification performance both in balanced and un-balanced scenarios, as well as in the case we restrict the analysis to specific geographical areas. Our analysis predicts that territories with larger income per capita, soil consumption, concentration of real estate activities and commuting network centrality in terms of closeness and Pagerank constitute the set of most affected areas, experiencing the strongest reduction of economic activities during the COVID-19 pandemic. Overall, we provide an application of an early-warning system able to provide timely evidence to policymakers about the detrimental effects generated by natural disasters and severe crisis episodes, thus contributing to optimize public decision support systems.


Introduction
The necessity to safeguard public health objectives represents the main priority for decision-makers during global emergencies such as the COVID-19 pandemic. Hence, data-driven expert and intelligent systems have been devised to help health organizations predict mortality (Smith & Alvarez, 2021), forecast new cases (Ertem, Araz, & Cruz-Aponte, 2022a) and manage hospital readmissions (Davazdahemami, Zolbanin, & Delen, 2022).
Next to public health, the most relevant concern for governments and organizations is to preserve the economic well-being of individuals while keeping them safe. This requires the collection and integration of multi-faceted data on individuals' behaviors, the formulation of alternative policy scenarios for intervention, the evaluation of their potential outcomes, and their implementation. All these steps identify elements of the development of decision support systems for public governance (Ibrahim & Larsson, 2019) and for evidence-based policy-making (Howlett, 2019).
Furthermore, the emergency context of pandemics increases the uncertainty and complexity of these strategic choices. On the one hand, this is due to the fact that economic signals are not available with * Correspondence to: Department of Electronics, Information and Bioengineering, Politecnico di Milano, Via Giuseppe Ponzio, 34, 20133, Milano, Italy. E-mail addresses: francesco.pierri@polimi.it (F. Pierri), francesco.scotti@polimi.it (F. Scotti), giovanni.bonaccorsi@polimi.it (G. Bonaccorsi), andrea.flori@polimi.it (A. Flori), fabio.pammolli@polimi.it (F. Pammolli). sufficient frequency and territorial granularity. On the other hand, it is difficult to disentangle all the factors contributing to individuals' decisions in the complex and challenging environment resulting from emergencies.
Finally, it is extremely relevant that the employed methodology and algorithms are transparent with respect to their decision criteria and as easily interpretable as possible (Mitchell, Potash, Barocas, D'Amour, & Lum, 2021;Shneiderman, 2020).
In this work, we tackle the problem of predicting the economic resilience of territories inside a specific country during emergency situations.
In terms of regional resilience, we refer to ''the ability of a region to recover successfully from shocks to its economy that either throw it off its growth path or have the potential to throw it off its growth path'', consistently with the definition provided by Hill, Wial, and Wolman (2008). Despite the growing popularity of the concept of regional resilience, due to recent global crises and Grand Challenges affecting our society, main existing contributions focused on theoretical frameworks, discussing the process through which territories can effectively deal with large scale shocks, and potentially transform their https://doi.org/10. 1016/j.eswa.2023.120803 Received 3 March 2023; Received in revised form 19 April 2023; Accepted 8 June 2023 economic structure towards new multiple equilibria (Boschma, 2015;Christopherson, Michie, & Tyler, 2010;Hassink, 2010;Pike, Dawley, & Tomaney, 2010;Simmie & Martin, 2010).
Due to recent frequent and hazardous natural disasters, empirical literature on regional resilience also paid attention to those socioeconomic factors that may reduce the negative impact of natural disasters such as floods (Liu et al., 2020 and earthquakes (Du, Wang, Xie, & Dong, 2023;Lu, Li, Mao and Wang, 2022). Moreover, authors focused on specific sectors, such as Transportation (Potter, Soroka, & Naim, 2022) and Energy (Erker, Stangl, & Stoeglehner, 2017), that offer critical services to regions and that may be particularly vulnerable to the disruption of infrastructure networks. Conversely, only a few studies analyzed the drivers of regional resilience during the healthcare emergency caused by the COVID-19 pandemic.
For instance,  highlight how cultural factors such as dialect diversity may reduce regional mobility resilience, whereas areas characterized by higher cultural tightness experienced stronger social connection during the health crisis. Hu, Li, and Dong (2022) show how a larger share of secondary and tertiary activities, a lower diversification of the economic structure and a higher degree of economic openness characterize a lower regional resilience in terms of contagions. Available studies, thus, mainly consider regional resilience in terms of resistance to new infections, whereas they neglect the economic performance of territories. Furthermore, Costa, Sallusti, Vicarelli, and Zurlo (2022) show that the size, sector and digitalization level of Italian companies are key drivers of their resilience. However, such analysis is performed at the firm level, while policymakers may benefit more from an effective decision support system developed at the territorial level, a granularity that is more consistent with the implementation level of restrictions and support policy packages. 1 To this end, we propose an almost in real-time indicator of the economic performance of Italian municipalities using the variation in consumption of citizens during the COVID-19 pandemic. We then collect different socioeconomic features of these areas and predict their performance during the healthcare crisis with a wide set of machine learning techniques. In particular, we cover the economic, social, environmental, infrastructure and demographic dimensions coherently with the empirical approaches adopted in a set of recent studies assessing urban and regional resilience (Cardoni, Noori, Greco, & Cimellaro, 2021;Shi et al., 2021;Yi, Wang, Li, & Dong, 2023;Zhang & Wang, 2023). Finally, we explore the prediction results with the use of Shapley values, in line with Explainable AI principles (Kim, Park, & Suh, 2020;Rai, 2020;Shin, 2021). In this way, we provide a practical application of how policymakers may benefit from the usage of decision analytical tools to increase the efficiency and effectiveness of their functionalities and decision processes (Phillips-Wren, Daly, & Burstein, 2021).
We show that network-based measures combined with economic, social, demographic and environmental variables are relevant predictors of the economic performance of Italian municipalities during the COVID-19 pandemic. In particular, territories with higher income per capita, soil consumption, concentration of real estate activities and network centrality in terms of closeness and page rank display a lower level of economic resilience to the healthcare emergency. We conduct classification experiments in three different settings, first focusing on predicting the most and least affected territories overall, then experimenting with an unbalanced scenario, and finally considering a use case based on the North/South divide.
1 For instance, before the start of the lockdown (9th of March 2020) specific restrictions targeted some Italian municipalities displaying high levels of contagion. Moreover, since the 6th of November 2020 the Italian government differentiated the stringency level of restrictions at the regional level based on the local pandemic intensity, adopting a three-tier risk framework. For further details, see also the following link: https://www.gazzettaufficiale.it/ eli/gu/2020/11/04/275/so/41/sg/pdf. Overall, our approach allows a decision-maker to identify which specific territorial features are reliable predictors of the performance of territories during emergencies such as the COVID-19 pandemic. Furthermore, our findings may be informative for a wide set of public stakeholders interested in early warning systems providing evidence on the resilience and vulnerability of economic activities at the local level during relevant crisis episodes.

Related work
The intersection of social, economic and environmental crises has enhanced the perceived sense of vulnerability of our society paving the way to new research paths on the concept of regional resilience (Hudson, 2010). In the field of regional resilience, when referring to exogenous unexpected and disruptive shocks, the concept of resilience is associated with the resistance to an external disturbance measured as the speed of return to the equilibrium point and the ability to resume either their pre-shock levels or their pre-shock trajectories (Pendall, Foster, & Cowell, 2010).
As shown in Fig. 1, our paper contributes to the literature on regional resilience by providing a decision support system able to predict the vulnerability of territories to the COVID-19 pandemic based on their pre-existing socioeconomic characteristics, thus supporting policymakers in the design of more effective packages of support in periods of emergency.
Since climate change led to more frequent extreme events, the impact of natural disasters on regional resilience has been the object of recent scrutiny. Extant literature shows a more limited resilience to seismic events for territories with lower GDP per capita and higher levels of interconnections and interdependencies among local infrastructures (Du et al., 2023;Lu, Li et al., 2022). Similarly, higher levels of investments in flood prevention projects, a higher level of wealth, lower precipitation levels and larger availability of hospital beds are connected with better resilience performances to floods at the regional level (Liu et al., 2020. On the other hand, more limited empirical evidence studies the economic resilience of territories to recent shocks such as the COVID-19 pandemic Lu, Zhang et al., 2022). Consistently with our research framework, we aim to fill this gap, by providing an analysis of those factors contributing to the economic regional resilience of Italian municipalities during the COVID-19 pandemic (see Fig. 1).
The concept of regional resilience has been investigated under the economic perspective, highlighting that larger immigration flows, economy diversification, governance innovation, technological variety, knowledge network, quality of institutions, availability of skilled human capital and supportive financial systems are key drivers to reduce the vulnerability of territories to external shocks (Boschma, 2015;Corodescu-Roşca, Hamdouch, & Iaţu, 2023;Hassink, 2010;Huang, 2022;Simmie & Martin, 2010;Tan, Hu, Hassink, & Ni, 2020). In line with such literature, we exploit economic, environmental, demographic and social drivers (Cardoni et al., 2021;Hu et al., 2022;Lin, Peng, Chen, & Zhang, 2022;Liu et al., 2020Liu et al., , 2021Lu, Li et al., 2022;Lu, Zhang et al., 2022;Tan et al., 2020;Ye, Hu, Lu, Dong, & Gu, 2022) to explain the economic resilience of Italian municipalities. Furthermore, differently from extant studies, we use mobility network centrality indicators due to the key role of mobility in the spread of contagion and for the resilience of economic activities.
We mainly differentiate from previous studies, since our measure of resilience is not based on an aggregation of those factors aiming to explain resilience as in Cardoni et al. (2021), Lin et al. (2022), Liu et al. (2020Liu et al. ( , 2021, Lu, Li et al. (2022), Lu, Zhang et al. (2022) and Ye et al. (2022). Furthermore, we do not rely on measurements of resilience based on low-frequency macroeconomic variables as Hu et al. (2022), Martin, Sunley, Gardiner, and Tyler (2016) and Tan et al. (2020). Instead, we use as an outcome of resilience the variation of  We underline and write in bold the dimensions that are included in our research framework and we highlight in blue other contributions addressing in a similar way such points in the regional economic resilience literature. economic consumption at the municipal level obtained through high frequency transaction data.
We use such data since the lack of immediately available and structured information on the economic damages produced by natural disasters and severe crisis episodes hampers near real-time analyses of the resilience of territories hit by shocks and of the speed of their recovery process.
We aim to introduce transactional data in the decision support system literature since empirical evidence has proved how this type of information is able to immediately capture the rapid contraction of household consumption during the pandemic, with a larger reduction in service sectors largely relying on social interactions such as Accommodation, Entertaining, Restaurants and Transportation (Carvalho et al., 2020;Dunn, Hood, Batch, & Driessen, 2021;Hacioglu, Känzig, & Surico, 2020). Furthermore, consumption variation across territories observed through transactional data might be explained by specific characteristics of local cardholders such as the income level and the age profile, with areas concentrating high-income users and younger population cohorts experiencing a lower economic activity during the pandemic (Baker, Farrokhnia, Meyer, Pagel, & Yannelis, 2020;Carvalho et al., 2020;Hacıoğlu-Hoke, Känzig, & Surico, 2021;Martin, Markhvida, Hallegatte, & Walsh, 2020;Sheridan, Andersen, Hansen, & Johannesen, 2020).
In addition, in line with our research purpose, transactional data may provide relevant insights into the resilience of economic systems against the pandemic, since sectors that were less affected by restrictions due to the provision of primary need goods -such as Grocery -and activities where consumers may easily shift to online shopping exhibited a lower economic fragility (Bounie, Camara, & Galbraith, 2020;Chen, Qian, & Wen, 2021;Chetty, Friedman, Hendren, Stepner, & Team, 2020;Dong, Gozgor, Lu, & Yan, 2021). Moreover, transactional data may help explain the impact of pandemic restrictions, based on the analysis of the ex-ante purchasing behavior of cardholders. In fact, individuals displaying larger expenditures in sectors characterized by a stronger restriction level show a higher reduction of economic consumption (Andersen, Hansen, Johannesen, & Sheridan, 2022). For such reasons we expect transactional data to be a relevant source of information to forecast the economic impact of the pandemic across territories.
From a methodological perspective, our paper employs K-Nearest Neighbor (K-NN), Random Forest and Support Vector Machine (SVM) classifiers, coherently with previous studies in the regional resilience literature aiming to explain the factors contributing to reduce the vulnerability of territories to natural disasters (Bassiouni, Chakrabortty, Hussain, & Rahman, 2023;Liu et al., 2020Liu et al., , 2021. Furthermore, we apply Logistic Regression, Multi-Layer Perceptron (MLP) and Network Analysis which have been recently applied also by other contributions in the decision support system literature to investigate COVID-19 emergency department readmissions (Davazdahemami et al., 2022), the impact of the pandemic on the tourism sector (Kumar, Misra, & Chan, 2022), and the variation of the centrality of economic sectors during the health emergency with respect to business as usual scenarios (Mascaretti et al., 2022).
As our main aim is to classify the economic performance of Italian municipalities, we do not rely on simple correlation analysis, factor analysis and linear regression methods, that have been used e.g. by Hu et al. (2022), Tan et al. (2020) and Ye et al. (2022) in the regional resilience literature. Furthermore, in line with Davazdahemami et al. (2022), we use Shapley Values to explain the relevance of various drivers of the economic resilience of territories differently from Cardoni et al. (2021) which exploit survey data.
Finally, we differentiate from the work of Ertem, Araz, and Cruz-Aponte (2022b) since they rely on specific empirical models related to the specific domain of study (regional probabilistic seismic hazard analysis).
Overall, we develop a framework of analysis that allows us to forecast almost in real time the fragility of territories based on their socioeconomic characteristics. Our study may thus support policy makers to identify the main drivers of resilience and investment areas to reduce the vulnerability of territories to specific types of shocks (see Fig. 1). Through the identification of relevant characteristics of territories able to explain regional resilience, measured based on high frequency transaction data, we assess the most relevant factors affecting the economic responses of Italian municipalities during the pandemic. Our study contributes to the development of a decision support system allowing regulators to deal with emergency periods underlining the advantage of data-driven and evidence-based decisions.

Bank transaction data
We build our primary variable of interest (see Section 3.2) from measurements of transactions made by customers of Intesa Sanpaolo, one of the major Italian banks, accounting for a market share larger than 15% in each region. The data tracks the aggregated daily value of transactions made by debit/credit cardholders in Point-of-Sales (POSs) of different Italian municipalities. Each data point represents the number and the value of transactions made by bank clients in all POSs of a specific municipality on a certain day. We use data regarding all merchant categories and all transactions made in-person (i.e., no online purchases). Furthermore, the data comprehend the entire Italian territory and the full period of analysis without missing days, however, we restrict the analysis to municipalities with an overall population greater than 1000 since smaller municipalities often exhibit extreme variability in consumption which is not representative of the entire sample. Overall, our dataset captures over 1 billion bank transactions for a total associated value of approximately 48 billion EUR.

Dependent variable of interest
Even though a considerable effort has been devoted to fostering data-sharing initiatives (Scotti, Pierri, Bonaccorsi, & Flori, 2022), the evaluation of the socio-economic impact of catastrophic events such as the COVID-19 pandemic is complex due to the systemic nature of the threat, affecting multiple systems at the same time, and the availability of real-time data (Hynes, Trump, Love, & Linkov, 2020). Most of the time, in similar situations, the evaluation of the effects of these events is obtained in two ways: either by evaluating them ex-post using official statistics or by using a proxy. In the first case, estimates are reliable but often arrive too late and with insufficient detail for policy-makers to track the evolution of the crisis and act immediately (Chetty et al., 2020). Regarding the second case, instead, valuable proxies of the impact of natural disasters have been obtained with social network data (Kryvasheyeu et al., 2016) and mobility data (Buckee et al., 2020). In these cases, estimates are available in nearly real-time and with high granularity, however, they represent an approximation of the real effects of the disaster and are affected by measurement and representativity biases (Wesolowski, Eagle, Noor, Snow, & Buckee, 2013). In this work, we propose a new source of data to track the effect of catastrophic events in Italy, bank transaction data. This source represents an improvement with respect to both official statistics and usual proxies. In fact, on the one hand, it is available in real-time and with sufficient granularity. On the other hand, it is better than other approximations of the socioeconomic impact of COVID-19 since it directly measures the variation in consumption of bank clients. Variation in consumption does not cover the entire effect of the pandemic but constitutes a relevant part of it since consumption represents a significant share of the overall economic output of a country (Carvalho et al., 2020). Finally, Scotti, Flori, Bonaccorsi, and Pammolli (2023) show that the data tracks closely the evolution of consumption in official sources, hence confirming that the data are representative of the overall Italian population.
In line with a large stream of literature on transaction data (Carvalho et al., 2020;Chen et al., 2021;Sheridan et al., 2020), we first define the weekly Year-on-Year (Y-o-Y ) change in the aggregate consumption, i.e., the variation of the total value of transactions performed by cardholders in POSs located in a given Italian municipality. In formula: , where , , is the total weekly economic consumption calculated in a 7-day time window in municipality , in the week of the year .
Finally, we define our main dependent variable as the average of the Y-o-Y variable in the weeks after the period of restrictions. Let lock be the week in which the lockdown is applied, then for each municipality i we define: We exclude the week when the lockdown was put in place because the decrease in mobility was not abrupt, and the effect of mobility restrictions were more visible in the following week (Bonaccorsi et al., 2020(Bonaccorsi et al., , 2021. Gray areas show municipalities for which data is not available whereas white areas show municipalities present in our data but excluded from our analysis either because their population is below 1000 or because the Y value is an outlier (i.e., smaller than the 1th percentile and larger than the 99th percentile).   Fig. 2 provides the geographical distribution of our dependent , variable. The first COVID-19 case was reported in Italy on February 17th, 2020; social distancing restrictions were quickly introduced in the Northern regions in order to prevent localized clusters of contagion. Finally, on March 9th a full lockdown was imposed on the entire country, restricting movements between and within regions until May 3rd (excluded). We can notice that the mobility restrictions imposed during the full lockdown caused a strong decrease in consumption.

Socioeconomic predictors
In terms of predictors of the resilience of Italian municipalities, we include those factors that may be connected with the impact of the COVID-19 pandemic on the reduction of economic activities.
Due to the heterogeneity in the impact of restrictions across territories, that may have been associated with alternative pre-existing levels of wealth and disparities (Bonaccorsi et al., 2020(Bonaccorsi et al., , 2021Polyakova, Kocks, Udalova, & Finkelstein, 2020-11-10), we include Income per capita and Inequality in the distribution of income in our prediction model. Inequality is computed as the ratio between the mean and the median income per capita of each Italian municipality. Since the distribution of income is usually skewed to the right, the former value will be usually larger than the latter, hence higher inequality will be signaled by a greater mean-to-median ratio. Furthermore, the policy outcome may be also driven by the size of municipalities (Ascani, Faggian, & Montresor, 2020;Hacıoğlu-Hoke et al., 2021). Therefore, we consider the municipality surface (Area) and the number of local inhabitants (Population) as additional predictors.
Restrictions have targeted economic sectors with different levels of stringency depending on the level of human interactions and their expected contribution to the resurgence of the contagion (Ascani, Faggian, Montresor, & Palma, 2021;Di Porto, Naticchioni, & Scrutinio, 2022). Consequently, our model takes into account the proportion of workers employed in essential sectors (Essential Workers (%)) that were allowed to stay open even during the lockdown. 2 Moreover, the digitalization level of business activities, the education level and ICT skills of employees contribute to the possibility to work remotely (Barbieri, Basso, & Scicchitano, 2022; Espinoza & Reznikova, 2020). For such reasons, we consider the estimated proportion of workers that are expected to work remotely (Teleworkers (%)) in our model. The Telework variable represents the sum of the proportion of workers employed in sector ( , ) in municipality weighted by a coefficient approximating the percentage of employees whose job may be performed remotely in sector ( ℎ ) estimated by Espinoza and Reznikova (2020). In formula: Our model contains also a measure of the percentage of residents with a fast internet connection (Broadband) to take into account the extent to which territories may be digitally connected. This is computed as the portion of residents that are reached by an Internet Connection whose speed is larger than 100 Mbps.
Since the technological and labor market structure can significantly influence the policy outcome in terms of contraction of the economic activity (Ascani et al., 2020;Ferraresi et al., 2021;Mascaretti et al., 2022), we include also the proportion of workers employed across different NACE sectors at 1 digit level of detail with respect to the total number of employees in the corresponding municipality. Table 1 reports the list of sectors used in our analysis. Furthermore, we account for the Herfindahl-Hirschman index of employees across different sectors, as a measure of economic diversification (HHI workers). We compute it across different sectors as the squared ratio between the number of employees in sector in municipality and the total number of employees in municipality . In formula: Moreover, we consider the percentage of people below 15 years old and above 65 years old as the local proportion of individuals not involved in business activities due to age constraints (Unoccupied(%)). In this way, we account also for demographic differences across municipalities that may have affected the pandemic intensity due to the heterogeneous impact of the virus across population cohorts with different age (Acemoglu, Chernozhukov, Werning, Whinston, et al., 2020;Baker et al., 2020;Fabbri, Gozzi, & Zanco, 2021;Makris, 2021;Sheridan et al., 2020). We then consider a list of factors linked to environmental characteristics of Italian municipalities, that may have fostered or reduced the diffusion of the virus (Lim, Kweon, Kim, Kim, & Lee, 2021;Weaver, Head, Gould, Carlton, & Remais, 2022) and interacted with other natural threats (Djalante, Shaw, & DeWit, 2020). Such factors may provide relevant insights regarding the anthropization level of local territories. For such reasons, our model encompasses the Soil Consumption (%), the Hydro Hazard, the Landslide Hazard and the Seismic Hazard as additional predictors.

Mobility network centralities
Several contributions have highlighted the role of mobility as a relevant predictor of both the diffusion of contagion (Chang et al., 2020;Chinazzi et al., 2020;Oliver et al., 2020) and the resilience of territories during emergencies (Bonaccorsi et al., 2020;Song, Zhang, Sekimoto, & Shibasaki, 2014). For these reasons, we extend our set of socioeconomic features by computing several node centrality measures that describe the importance of each municipality in the Italian network of commuting mobility.
We leverage official measurements of commuters traveling between municipalities daily, for work or study reasons, produced during the last Italian national census (2011) and made available by the Italian National Institute of Statistics (ISTAT). 3 We represent the network as a weighted directed graph ( , ) where vertices (or nodes) represent Italian municipalities and edges represent the number of individuals that declared to commute between two municipalities on a daily basis. The network contains 8083 nodes and 555,906 edges, for a total of 28,802,640 individuals moving between Italian municipalities (including self-loops), nearly half of the 2011 Italian population. We compute the following node centrality measures: • Betweenness: it is a measure of centrality in a graph based on shortest paths, computed as the sum of the fraction of all-pairs shortest paths that pass through each node. In a weighted graph, the shortest path is the path that connects two nodes with the minimum sum of weights on the edges of the path. This centrality quantifies the probability for a node to act as a bridge along the shortest path between two other nodes (Freeman, 1977).  • Closeness: it is also a measure of node centrality based on shortest paths, computed as the reciprocal of the sum of the length of the shortest paths between the node and all other reachable nodes in the graph. It indicates how close is a node to other nodes in the graph (Okamoto, Chen, & Li, 2008). • Eigenvector: it is a measure of node centrality based on the centrality of neighbors. The eigenvector centrality for node is the th element of the vector defined by the equation = where A is the adjacency matrix of a graph. A high eigenvector score means that a node is connected to many nodes that have high scores themselves (Bonacich, 2007).
• Pagerank: it is a measure of centrality introduced to rank webpages in search engine queries (Page, Brin, Motwani, & Winograd, 1999), which accounts for both the number and quality of links to a node, assuming that more important nodes will likely be more reachable from other nodes.
In particular, in order to compute shortest paths for Closeness and Betweenness centralities we consider the distance between two nodes as the reciprocal of the weight on the edge connecting the two nodes, similar to Bonaccorsi et al. (2020) and Latora and Marchiori (2007). We exclude basic centrality measures such as Degree -the number of edges connected to a given node -and Strength -the sum of the weights on the edges connected to a given node -as they are highly correlated with socioeconomic variables such as population and density, whereas we aim to capture non-trivial characteristics of the mobility network of Italian municipalities.

Dataset summary
The total number of observations (municipalities) in our sample is 4672. Each observation is encoded in a vector X with 34 dimensions (covariates) plus the Y variable described above, which is used to derive the class label as defined next. We show pairwise correlations of these 35 variables in Fig. 4. We also provide their descriptive statistics in Tables 2 and 3.

Machine learning techniques
We rely on several open-source packages to perform our classification tasks. We encode each municipality with a vector of socioeconomic and network-based covariates and we assign them a binary label based on their Y variable in several different settings.
We employ scipy (Virtanen et al., 2020) to carry out a hierarchical clustering of municipalities based on their Z-standardized features, using a cosine distance and average linkage to obtain clusters. Similar results are observed when using euclidean distance or other linkage approaches, e.g., complete or Ward linkage. We use scikit-learn (Buitinck et al., 2013) to train several binary classifiers, namely K-Nearest Neighbors, Logistic Regression, Multilayer Perceptron (MLP), Random Forest and Support Vector Machine (SVM). We apply Z-standardization of features and use default hyperparameters as specified in the package. We also employ automl (Feurer et al., 2015) to find the optimal configuration for the Random Forest classifier.
We compute SHAP (SHapley Additive exPlanations) values to assess the importance of features in the binary classification setting. This is a game theoretic approach that explains the output of machine learning models by using Shapley values from game theory (Lundberg & Lee, 2017a). In particular, we leverage the shap Python library (Lundberg & Lee, 2017b).

Evaluation metrics
We use the following evaluation metrics to assess the performance

Area Under the Precision-Recall Curve (AUPRC): the Precision-Recall curve shows the trade-off between Precision and Recall values as the classification threshold value is varied.
AUPRC is a widely adopted measure of success of prediction in settings when the classes are very imbalanced. A high area under the curve represents both high recall and high precision.
We evaluate the results of hierarchical clustering through the average purity score of clusters obtained at each step of the clustering procedure. Purity is computed by (1) assigning each observation to the class which is most frequent in the cluster, and (2) counting the number of correctly assigned observations within each cluster and dividing by the total number of observations. Higher levels of purity score in a balanced setting indicate more homogeneous clusters.

Results
In this section we present results from three different classification experiments. A detailed workflow of our analysis is shown in Fig. 3.

Classification of territories in the tails of the distribution
We formulate our main binary classification problem by considering the tails of the distribution of Y. Specifically,for ∈ [10,15,20,25,30,35,40,45,50], we define two classes of the most and the least affected territories that correspond respectively to those with Y in the bottom and in the top K% of the distribution. For instance, as shown in Fig. 5, if = 20 we classify territories in the bottom 20% versus those in the top 20% of the distribution of Y. This labeling procedure yields balanced samples of increasing size (e.g., = 10 ⇒ 936, = 15 ⇒ 1402, up to all 4672 observations with Y below/above the median).
As a preliminary step, we perform a hierarchical agglomerative clustering by computing the cosine similarity between Z-standardized values of vectors encoding Italian territories. To understand whether observations from the same class tend to cluster together and be well separated from the opposite class, thus indicating that a classifier could reach high levels of accuracy, we evaluate the purity score. Such an indicator is computed when choosing an increasing number of clusters in the hierarchical clustering, i.e., by lowering the distance threshold to merge clusters in the agglomerative clustering procedure. We report F. Pierri et al. two findings, as shown in the bottom panel of Fig. 5: first, for all choices of , the purity score is directly proportional to the number of clusters, as expected by the definition of the purity score; second, clustering observations labeled using small values of (e.g. bottom 10% versus top 10%) yields generally higher purity scores compared to large values of (e.g. bottom 40% versus top 40%), indicating that it might be harder to distinguish the most and least affected territories as we consider larger samples of the dataset. This is interesting, as it means that municipalities that showed a very different economic performance during the pandemic (i.e., at the two extreme tails of the Y distribution) also exhibit stronger differences in terms of pre-existing socioeconomic characteristics.
Next, we train and test several off-the-shelf binary classifiers, as specified in Section 3.6, for different choices of . We show AUROC and F1-score for all classifiers and for increasing values of respectively in Tables 4 and 5. Notice that the precision, recall and F1-score metrics are all equivalent when the classes are perfectly balanced, hence we show only the last metric. Given the balanced dataset, we first notice that all classifiers perform better than a random baseline, which would exhibit evaluation metrics equal to 0.5 in all cases (we omit the baseline in the tables for clarity). We also observe that, for all models, classification performance degrades as we consider increasing values of , as expected per the preliminary clustering analysis. Besides, the Random Forest classifier exhibits the best performance (with the highest values for = 15) in both evaluation metrics and for all choices of , whereas the single-layer neural network (MLP) exhibits the poorest performance. Linear classifiers such as the Logistic Regression and SVM (linear) also exhibit performance that is only slightly inferior to the Random Forest classifier.
We further explore the performance of the two best classifiers, namely Logistic Regression and Random Forest, by performing a hyperparameter tuning. For the former we use the LogisticRegressionCV function of scikit-learn, which performs an automatic grid search over the parameters. For the latter, we use the autoML Python library and select the best Random Forest classifier with cross-validation based on the AUROC value. In the case of Logistic Regression we do not report differences compared to the model with default hyperparameters. We compare performance of the Random Forest classifier with off-the-shelf parameters versus the optimal configuration obtained through AutoML in Table 6. We notice only a slight improvement for different values of , with a positive gain in the range [+1.28%, +4.93%]. To understand the importance of each feature in the classification task, we compute SHAP values (Lundberg & Lee, 2017a) for the best Random Forest classifier obtained through hyperparameter tuning and evaluated with cross-validation ( = 10 folds and 80/20 train/test size ratio) on the data with = 15, where we observed the best classification performance. We show results in Fig. 6. We can see that the features with the largest SHAP magnitude include incomerelated variables (Income per Capita and Inequality) but also centrality measures in the mobility network (Closeness and PageRank) as well as features describing environmental characteristics (Soil Consumption) and the productive structure of territories (Sectors L and F). SHAP values with a positive (negative) magnitude indicate which levels of features characterize observations that are classified as most affected (least affected) territories; e.g., low values of Income per Capita are associated with territories that perform better in terms of the Y variable. In addition to SHAP values, we also computed the cross-validated feature importance of the Random Forest classifier (for = 15), noticing that the ranking is very similar to the one resulting from SHAP values (especially the most relevant features), thus strengthening the importance of the highlighted variables.
Accordingly, we show distributions of the top 4 relevant features for the two classes obtained with = 15 in Fig. 7. We can see that, in accordance with SHAP values, most affected territories exhibit higher values of the four features. All distributions are statistically different according to Kruskal-Wallis tests ( < 0.001, = 0.05).
We carry out an ablation study to interpret the model performance by evaluating the cross-validation AUROC and F1-Score of a Random Forest classifier trained and tested on an increasing number of features F. Pierri et al. (sorted by their Shapley value in the case when = 15) and for different values of . As shown in Fig. 8, the classifier yields the best performance for = 15 and the entire set of available features, in accordance with Tables 4 and 5. However, we observe a very similar performance when employing as few as 10 features, suggesting that the most relevant covariates -as indicated by the SHAP summary plot (Fig. 6) -are particularly more informative than the others and that a very simple model can accurately distinguish the two classes of territories.

Classification of territories in a given interval of variation in y
In our second classification task, we consider specific (symmetric) intervals in the distribution of Y in order to define the most and least affected territories. Specifically, as shown in Fig. 9, we consider territories with Y in the interval [−0.75, −0.25] (respectively [+0.25, +0.75]) as the most (respectively least) affected territories. This labeling yields an unbalanced dataset with 1725 versus 548 observations in the two classes (almost a 3/1 ratio).
We do not perform the clustering exercise in this case as the classes are highly imbalanced and the purity score is not suitable, as it would be biased towards the majority class. We, therefore, train and test the same set of binary classifiers, as in the previous task, and we show their performance in terms of weighted Precision, Recall and F1-Score in Table 7. Notice that we compute the AUPRC instead of the AU-ROC, given the imbalanced dataset. We can observe that all classifiers perform better than a random baseline, with a gain in the range 15%-45% compared to a baseline classifier. Different classifiers are mostly    We compute SHAP values as in the previous classification task, finding a very similar ranking, which we omit for brevity.

Classification of north versus south territories
In our third classification problem, we address again the prediction of territories in the tails of the distribution of the Y variable for different values of K, but considering two separate tasks (and thus training separate classifiers) for, respectively, municipalities in the Northern and Southern regions of Italy. To classify Italian territories between North and South we use the Nomenclature of Territorial Units for Statistics (NUTS) adopted by the European Union (EU). For each EU member country, all its territories are classified in a hierarchy of three NUTS levels from most aggregated (NUTS 1) to least aggregated (NUTS 3). Specifically, based on NUTS 1, we label as ''North'' those municipalities in the North-West and North-East NUTS regions, and as ''South'' those that belong to the South and Island NUTS regions. These two groups, which comprise respectively 2306 and 1697 observations, are shown in Fig. 10, and exhibit different distributions of the Y variable, as highlighted in Fig. 11, with Southern territories displaying better economic performance, i.e., on average southern territories have experienced a smaller reduction in consumption with respect to northern ones. Moreover, the most and least affected municipalities in the North display less marked differences in terms of variation of economic consumption with respect to South areas.
In Fig. 12, we show the purity score plot for a hierarchical agglomerative clustering computed using cosine similarity between vectors of standardized features, separately for the two groups of territories and for different values of K. We can notice very different results between the two regional areas, with a better trade-off in purity score versus the number of clusters in the South scenario (purity∼0.7 with a dozen clusters) compared to the North (purity < 0.7 with a dozen clusters), suggesting that territories in the tail of the distribution are more distinguishable in the latter setting. Moreover, such a result F. Pierri et al. confirms that municipalities with very different economic performances also exhibit different characteristics.
We show the performance of a Random Forest classifier trained and tested on North and South territories separately, and for different values of , in Fig. 13. We omit other classifiers as they perform worse in all cases. As we can see, the classifier performs much better in the task of classifying territories in tails of the distribution in the South scenario compared to the North scenario. In both cases, we notice a similar decreasing trend in performance for increasing levels of , which, however, is less pronounced compared to Tables 4 and 5. Interestingly, we do not observe a maximum at = 15 but rather a minimum in the North setting.
Finally, we compute SHAP values for = 15 in order to compare results with the first classification setting. Fig. 14 shows the top 10 most relevant features. We can observe a partial overlap (5 out of 10 features) and different rankings, which again indicate very different socioeconomic characteristics of territories in the two geographical macro-regions. Interestingly, we also observe some of the same top relevant features as in the first classification setting (cf. Fig. 6). We show in Fig. 15 the distribution of the top 4 relevant features as in the first classification problem: we can observe that the distributions exhibit a higher discrepancy in the South setting compared to the North, thus explaining why the classification performance was higher in the former setting.

Discussion
This paper introduces a decision analytical framework to predict the economic resilience of Italian municipalities during the COVID-19 pandemic. We highlight that local pre-existing socioeconomic characteristics constitute relevant variables to forecast the resilience of territories to this healthcare emergency. Indeed, the most and least affected municipalities can be classified with adequate performance through off-of-the-shelf classifiers based on a set of features combining network-based measures with factors covering the social, economic, demographic and environmental dimensions. Across different estimated models, the most accurate results tend to be achieved through a random forest approach, although a similar performance is obtained with simpler linear-based models such as logistic regression and SVM.
We show that income per capita, network closeness, soil consumption, Pagerank centrality and employment concentration in specific sectors represent the most relevant drivers to predict the variation of   economic consumption across Italian municipalities. In particular, the least resilient territories are characterized by higher values of such variables. Furthermore, we highlight that our model properly classifies most and least affected territories also in unbalanced scenarios and when we restrict the analysis to specific geographical areas. Interestingly, better classification performances are obtained for municipalities in the South of Italy, where a larger distance in terms of economic performances between the most and least affected territories is associated with stronger differences in terms of pre-existing socioeconomic characteristics between the two groups. Our results confirm previous evidence from Costa et al. (2022) and Ferraresi et al. (2021) highlighting a lower resilience of the Construction sector and the relevance of the local economic structure to explain the vulnerability of economic systems to the COVID-19 pandemic, with sectors more affected by social distancing measures displaying higher economic losses. Our estimates are also in line with Cardoni et al. (2021) showing that more populated areas exhibit lower resilience levels during the COVID-19 pandemic. Furthermore, the link between lower economic resilience and higher income per capita might be explained by the larger number of infections observed in wealthier territories as shown by Ascani et al. (2020). Conversely, such evidence is in contrast with results from Liu et al. (2021) who highlight that areas with larger GDP per capita experience a faster recovery from floods. This difference may be explained by the heterogeneous nature of analyzed shocks, with drivers of resilience that may be different in case of a natural disaster with respect to a pandemic.
Our models corroborate that a larger presence of workers in essential sectors contributes to regional resilience, consistently with Ascani et al. (2021) and Ferraresi et al. (2021). We rather do not find significant lower vulnerability in case of higher economy diversification differently from Tan et al. (2020) and theoretical contribution from Boschma (2015), Christopherson et al. (2010) and Simmie and Martin (2010), thus highlighting the peculiarity of the COVID-19 pandemic with respect to other shocks.
Therefore, due to the specificity of the COVID-19 pandemic, we believe that relevant drivers of resilience to the healthcare emergency cannot be directly extended in case of other types of shocks. Nonetheless, our proposed framework of analysis could be adopted by policy makers to deal with alternative kinds of shocks including natural disasters, pandemics and other man-made extreme events. In particular, regulators should collect different factors influencing the resilience of territories. Using variables covering the demographic, economic, environmental, education, health, infrastructure and social dimensions should provide a flexible group of factors allowing to deal with the concept of resilience in a wide set of shocks. Furthermore, they should select a suitable measure of resilience. We propose a high frequency and spatially granular variable to allow almost real time analysis. However, alternatives may be represented by aggregation of input factors of resilience, depending on the data availability and urgency to obtain timely evidence. Moreover, a wide set of statistical modeling techniques can be employed to obtain robust evidence of the most relevant drivers of resilience and to identify those factors where policy makers could invest to reduce the vulnerability of territories to external shocks.

Limitations
There are some limitations to our approach. First, our dependent variable represents a proxy of the full shock associated with COVID-19, hence there are other socio-economical dimensions that may be better tracked with different indicators. Secondly, despite our dataset covering a significant share of the relevant dimensions of Italian territories, it could be expanded to include a greater set of variables, for instance regarding environmental performance and infrastructure endowment. Furthermore, in this work, we have attempted to balance classification performance with model explainability avoiding methods that are difficult to interpret, despite being more performant. Hence a different methodological approach may be adopted to achieve better classification results. Alternatively, we may design a causal framework that employs regression models to explain the individual performance of each territory as well as quantify how different socio-economic endowments contribute to different levels of economic resilience Finally, despite our choice of a flexible framework, there are aspects of our approach that are specifically tailored to the Italian context and to the COVID-19 emergency. A careful consideration of the variables to include in the model may be necessary to generalize our approach to other territorial contexts.

Conclusions
We presented a machine learning framework to predict the economic resilience to crises of Italian municipalities based on local preexisting socioeconomic characteristics, analyzing the first wave of the COVID-19 pandemic and leveraging a large-scale dataset of bank transactions. Overall, our analysis contributes to the literature on expert and intelligent systems for policy-making by providing an application that can support decision-makers in making optimal data-driven decisions during the early stages of relevant crisis episodes, exploiting available high-frequency real-world information.

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

Data availability
Data about economic transactions are confidential and cannot be publicly disclosed. Municipal data on income, population, employee by sector and commuting patterns are all available from the Istat website (http://dati.istat.it/). Data regarding employment by sector have been obtained from the National Registry of Active Firms (ASIA database) and are available to researchers upon request.