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

3.1. 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.

3.2. 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 $Spendin{g}_{i,w,Y}$ is the total weekly economic consumption calculated in a 7-day time window in municipality $i$, in the week $w$ of the year $Y$.

Finally, we define our main dependent variable $Y$ 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).

Fig. 2 (bottom panel) shows the temporal distribution of weekly Y-o-Y, computed for each municipality. The upper panel in Fig. 2 provides the geographical distribution of our dependent $Y_{i,t}$ 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.

3.3. 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., 2021, Polyakova et al., 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 et al., 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 et al., 2021, Di Porto et al., 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 et al., 2022, Espinoza and 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 $k$ ($Employee{s}_{i,k}$) in municipality $i$ weighted by a coefficient approximating the percentage of employees whose job may be performed remotely in sector $k$ ($Teleworkshar{e}_k$) 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 $k$ in municipality $i$ and the total number of employees in municipality $i$. 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 et al., 2020, Baker et al., 2020, Fabbri et al., 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 et al., 2021, Weaver et al., 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.

3.4. 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 et al., 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 $G(V,E)$ where vertices (or nodes) $V$ represent Italian municipalities and edges $E$ 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 $i$ is the $i$th element of the vector $x$ defined by the equation $Ax=\lambda x$ 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.

3.5. 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 Table 2, Table 3.

3.6. 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).

3.7. Evaluation metrics

We use the following evaluation metrics to assess the performance of machine learning classifiers (TP $=$ true positives, FP $=$ false positives, FN $=$ false negatives):

• 1.

  Precision$=\frac{TP}{TP+FP}$, the ability of a classifier not to label as positive a negative sample.

• 2.

  Recall$=\frac{TP}{TP+FN}$, the ability of a classifier to retrieve all positive samples.

• 3.

  F1-score$=2\frac{\text{Precision}\cdot \text{Recall}}{\text{Precision}+\text{Recall}}$, the harmonic average of Precision and Recall.

• 4.

  Area Under the Receiver Operating Characteristic curve (AUROC): the Receiver Operating Characteristic (ROC) curve, which plots the TP rate versus the FP rate, shows the ability of a classifier to discriminate positive samples from negative ones as its threshold is varied. Larger AUROC values correspond to better classifiers, with the perfect classifier yielding AUROC $=$ 1.

• 5.

  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.

4.1. 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 $K$ $\in$ [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 $K=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., $K=10\Rightarrow 936$, $K=15\Rightarrow 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 two findings, as shown in the bottom panel of Fig. 5: first, for all choices of $K$, 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 $K$ (e.g. bottom 10% versus top 10%) yields generally higher purity scores compared to large values of $K$ (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 $K$. We show AUROC and F1-score for all classifiers and for increasing values of $K$ respectively in Table 4, Table 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 $K$, as expected per the preliminary clustering analysis. Besides, the Random Forest classifier exhibits the best performance (with the highest values for $K=15$) in both evaluation metrics and for all choices of $K$, 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 $K$, with a positive gain in the range $[+1.28\text{\%},+4.93\text{\%}]$.

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 ($k$$=$ 10 folds and 80/20 train/test size ratio) on the data with $K=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 income-related 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 $K=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 $K=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 ($P<0.001$, $\alpha =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 (sorted by their Shapley value in the case when $K=15$) and for different values of $K$. As shown in Fig. 8, the classifier yields the best performance for $K=15$ and the entire set of available features, in accordance with Table 4, Table 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.

4.2. 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 AUROC, 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 equivalent to each other, with the Random Forest yielding the highest AUPRC value and the MLP exhibiting the worst performance for all metrics.

We compute SHAP values as in the previous classification task, finding a very similar ranking, which we omit for brevity.

4.3. 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$\sim$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 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 $K$, 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 $K$, which, however, is less pronounced compared to Table 4, Table 5. Interestingly, we do not observe a maximum at $K=15$ but rather a minimum in the North setting.

Finally, we compute SHAP values for $K=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.

5.1. 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.

5.2. Conclusions

We presented a machine learning framework to predict the economic resilience to crises of Italian municipalities based on local pre-existing 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.