key: cord-0529380-ygktav2g
authors: Diem, Christian; Borsos, Andr'as; Reisch, Tobias; Kert'esz, J'anos; Thurner, Stefan
title: Quantifying firm-level economic systemic risk from nation-wide supply networks
date: 2021-04-15
journal: nan
DOI: nan
sha: 7121f842694cb1974bc809fe0b53ca19c63aaef5
doc_id: 529380
cord_uid: ygktav2g

Crises like COVID-19 or the Japanese earthquake in 2011 exposed the fragility of corporate supply networks. The production of goods and services is a highly interdependent process and can be severely impacted by the default of critical suppliers or customers. While knowing the impact of individual companies on national economies is a prerequisite for efficient risk management, the quantitative assessment of the involved economic systemic risks (ESR) is hitherto practically non-existent, mainly because of a lack of fine-grained data in combination with coherent methods. Based on a unique value added tax dataset we derive the detailed production network of an entire country and present a novel approach for computing the ESR of all individual firms. We demonstrate that a tiny fraction (0.035%) of companies has extraordinarily high systemic risk impacting about 23% of the national economic production should any of them default. Firm size alone cannot explain the ESR of individual companies; their position in the production networks does matter substantially. If companies are ranked according to their economic systemic risk index (ESRI), firms with a rank above a characteristic value have very similar ESRI values, while for the rest the rank distribution of ESRI decays slowly as a power-law; 99.8% of all companies have an impact on less than 1% of the economy. We show that the assessment of ESR is impossible with aggregate data as used in traditional Input-Output Economics. We discuss how simple policies of introducing supply chain redundancies can reduce ESR of some extremely risky companies.

(Dated: April 16, 2021)

Crises like COVID-19 or the Japanese earthquake in 2011 exposed the fragility of corporate supply networks. The production of goods and services is a highly interdependent process and can be severely impacted by the default of critical suppliers or customers. While knowing the impact of individual companies on national economies is a prerequisite for efficient risk management, the quantitative assessment of the involved economic systemic risks (ESR) is hitherto practically non-existent, mainly because of a lack of fine-grained data in combination with coherent methods. Based on a unique value added tax dataset we derive the detailed production network of an entire country and present a novel approach for computing the ESR of all individual firms. We demonstrate that a tiny fraction (0.035%) of companies has extraordinarily high systemic risk impacting about 23% of the national economic production should any of them default. Firm size alone cannot explain the ESR of individual companies; their position in the production networks does matter substantially. If companies are ranked according to their economic systemic risk index (ESRI), firms with a rank above a characteristic value have very similar ESRI values, while for the rest the rank distribution of ESRI decays slowly as a power-law; 99.8% of all companies have an impact on less than 1% of the economy. We show that the assessment of ESR is impossible with aggregate data as used in traditional Input-Output Economics. We discuss how simple policies of introducing supply chain redundancies can reduce ESR of some extremely risky companies.

Keywords: systemic stability | production networks | shock propagation | cascading failure | network centrality measures Increasing the efficiency of production processes and corporate supply chains has been a dominating economic paradigm of the past decades. Popular managerial concepts that exemplify this view are reflected in keywords such as supply chain management, lean production (Lamming, 1996) , just-in-time delivery (Kannan and Tan, 2005) , out-and global sourcing (Hummels et al., 2001; Quélin and Duhamel, 2003; Trent and Monczka, 2003) , or supply base reduction (Choi and Krause, 2006; Trent and Monczka, 1998) . Efficiency gains are usually achieved by reducing inventory buffers, shorter lead times, supplier integration, or reducing the number of direct suppliers. Actions like these reduce production costs and increase profits. However, these actions also do have consequences in terms of resilience of the overall economy. It has been argued that increased levels of efficiency go hand in hand with a reduction of resilience (Choi and Krause, 2006; Craighead et al., 2007) .

Supply chains of firms and consequently their produc- * E-mail: diem@csh.ac.at; stefan.thurner@muv.ac.at tion processes are highly interdependent. Supplier-buyer relations between companies lead to so-called production networks. The ongoing transformation of production networks towards higher efficiency has made these networks more vulnerable to shocks (Craighead et al., 2007 ). On regional scales, hurricane Katrina or the Japanese earthquake in 2011 have shown the economic impacts that can arise due to subsequent cascading shock propagation along corporate supply chains (Carvalho et al., 2020; Hallegatte, 2008; Inoue and Todo, 2019) . The COVID-19 pandemic impressively revealed that not only overall economic activity can be affected by interruptions of supply chains, but they also may lead to shortages in basic supplies (Ivanov and Dolgui, 2020) , affecting people directly. This became apparent, for example, in food production (Bloomberg, 2020) , in vaccine supply (The Economist, 2021) , computer chips manufacturing, and car manufacturing (Financial Times, 2021a,b) .

The propagation of shocks through an economy is tightly related to classical input-output economics (Miller and Blair, 2009) , where, however, only industry sectors are studied. For details and relations to this work, see Appendix A. The importance of supply chains at the firm level has been considered by economists, in particular the consequences of individual company failure on the overall economy (Acemoglu et al., 2012; Bak et al., 1993; Gabaix, 2011) . The supply chain and production management literature studied firm level supply networks (Choi et al., 2001) and the spreading of disruptions along supply chains under various names, such as supply chain resilience (Craighead et al., 2007) , snowball effect (Świerczek, 2014) , the ripple effect (Ivanov et al., 2014) , or nexus suppliers (Shao et al., 2018) . Despite all these efforts, reliable and systematic estimates of firm level systemic risk for entire economies are hitherto not available.

When compared to recent progress in the assessment of financial systemic risk in financial networks, the quantification of economic systemic risk (ESR) of individual companies in production networks is still in its infancy. From initial demonstrations of the relevance of network structure in the context of financial systemic risk (Allen and Gale, 2000; Boss et al., 2004a; Elsinger et al., 2006) , by now it is not only possible to assign systemic risk to individual players in financial networks (Battiston et al., 2012; Thurner and Poledna, 2013) but also to individual transactions (Poledna and Thurner, 2016) , and on multiple network layers (Poledna et al., 2015) . These developments allow for novel policy paradigms for systemic risk regulation (Thurner, forthcoming) . To reliably assess systemic risks detailed and correct information on the underlying networks is essential (Battiston et al., 2016) . For more information on financial systemic risk and its relation to the topic of the present paper, see Appendix B.

It is impossible for decision makers to manage economic systemic risks proactively without understanding which companies pose exceptionally high risk to the entire economy in case of their (temporary) failure. It is therefore important to develop methodology to quantify these risks. Here we show how firm level supply transaction data can be used to compute an economic systemic risk index (ESRI) for every company within an entire country.

The most important ingredient for developing a meaningful ESRI is the underlying firm level production network, consisting of the supply relations between companies and their production processes. Country-wide production networks can be reconstructed from firm level value added tax (VAT) transaction data (Borsos and Stancsics, 2020; Dhyne et al., 2015) 1 . The construction of the supply network is schematically depicted in FIG.  1 (a) . For every sales transaction between a supplier i and buyer j, the monetary value of the goods and services sold, V ji , can be inferred, from the tax rate, τ , and the tax amount paid, T ji = τ V ji . We use this as an 1 The value added tax exists in all OECD countries, except for the USA (OECD, 2020).

Reconstructing the production network.

(a) Schematic supplier-buyer transaction between two companies, i and j. Supplier i produces product, pi, and delivers a quantity of Wij to buyer j. The buyer makes a gross payment consisting of the net price, Vji, and the value added tax, Tji. (b) Production network consisting of eight firms. Color represents their industry classification. Delivered goods and services are recorded in the weighted adjacency matrix, where the element Wij(t) records the volume delivered from i to j. The production function f6 for firm 6 is shown as an example. It uses two types of input, p2 = p4 and p7 from firms 2, 4, and 7 to produce the amount W68 of product p6. Firm 6 sells its entire output to firm 8. (c) Section of the Hungarian production network with 4,070 nodes and 4,845 links. Node size corresponds to the total strength (proxy for company size). Many supply "chains" are seen to be connected to form a supply network. estimate for the volume, W ij , of product type p i , delivered from supplier i to buyer j. Note that the proxy (volume = price × quantity) is a basic assumption in economics (Lequiller and Blades, 2014; Miller and Blair, 2009) . For notation, in-links to a node represent supply transactions (buying); out-links are sales transactions. We define the in-strength of node i as the sum of all its in-links, s in i = n j=1 W ji (volume of purchased products), the out-strength is the sum of all out-links, s out i = n j=1 W ij (sales). The (total) strength is defined as s i = s in i + s out i and serves as a proxy for firm-size. To assess the importance of the supply transaction, W ij , between firms i and j, it is essential to know which product type, p i , is exchanged and how it is used in the production process of firm j. We use the term "products" for goods and services. Economists typically resort to the simplifying assumption that every company produces one out of m different products that is determined by the company's industry classification. We use an industry affiliation vector, p, that assigns one of m industries to each firm i, p i ∈ {1, 2, . . . , m}. We use the NACE (Statistical Classification of Economic Activities in the European Community (EUROSTAT, 2021)) classification scheme on the 4-digit level with m = 615 categories. For more details on industry classifications, see Appendix C.

The production process of a company is commonly described with a production function (Carvalho and Tahbaz-Salehi, 2019) , f i , that determines the (maximal) amount of product x i (output) of type p i that firm i can produce with a given amount of intermediate products, Π i = (Π i1 , Π i2 , . . . Π im ) (Π ik is the amount of input of type k of firm i), its employees (labour), l i , and manufacturing equipment (capital), c i . We use production functions that allow us to determine how much firm i can still produce if a supplier j fails to deliver its products of type p j = k to firm i and hence reduces the availability of input Π ik . This is illustrated with a small production network in FIG. 1 (b) . In-links to a node (arrows pointing towards node) represent supply transactions (buying), out-links are sales transactions. By identifying the supply links belonging to the inputs in the production function, it can be seen how production depends on the current state of the network. This is illustrated for firm 6. It uses two types of inputs, product a (supplied by firm 2 and 4) and product b (supplied by firm 7), and produces an amount x 6 (t+1) = f 6 W 26 (t)+W 46 (t), W 76 (t) of product c. The delivery from 6 to 8, W 68 (t+1), depends on the in-links of 6, W 68 (t + 1) = f 6 W 26 (t) + W 46 (t), W 76 (t) . Should one of the in-links diminish or vanish, this impacts the output of 6. This illustrates a so-called downstream shock propagation of a supply shock. Vice versa, if 8 decides to no longer buy from 6, the in-links W 26 (t), W 46 (t), W 76 (t) would no longer be needed. This illustrates the upstream propagation of a demand shock.

The specific choice of the production function, f i , determines the intensity of the downstream shock propagation. Frequently used production functions include the constant elasticity of substitution (CES) (Carvalho et al., 2020; McFadden, 1963; Moran and Bouchaud, 2019) and its special cases, the Cobb-Douglas (Acemoglu et al., 2012) and the Leontief (Inoue and Todo, 2019; Pichler et al., 2020) production functions, see Appendix D for details. Here we take a short term perspective and consider a generalized Leontief production function (GL) that accounts for the fact that even for companies with physical production processes, not all procured inputs are essen-tial in the short term. Note that provides a study on which inputs are essential for 56 industry sectors. For example, in tire production, rubber intermediates are essential inputs, while consulting services are not (short term). The GL treats non-essential inputs in a linear fashion, while essential inputs are treated in the non-linear Leontief way. We denote the set of essential inputs by I es i and non-essential inputs by I ne i , respectively. We define the GL as

where α ik are technologically determined coefficients and β i is the production level that is possible without nonessential inputs k ∈ I ne i ; α i is chosen to interpolate between the full production level (with all inputs) and β i . Both, α and β, are determined by W , I es i and I ne i . Note, that we assume labour and capital to be fixed (in the short term) and omit them in the notation. The Leontief production function is a special case of Eq. (1) if all inputs are essential. The linear production function is the special case when all inputs are non-essential; see Appendix D. An essential aspect that determines downstream shock propagation is how easily failing suppliers can be replaced. For our purposes we assume that companies with a low market share are easier to replace, than firms with large market shares; see Appendix E.

For the empirical analysis we calibrate the GL in Eq. (1) to four scenarios based on the firms' NACE classifications. First, for a hypothetical purely linear production scenario (LIN), all firms have liner production functions, i.e. only non-essential inputs, or I ne i = {01, . . . , 99}). Second, in a purely Leontief scenario (LEO), all firms have Leontief functions, i.e. only essential inputs, I es i = {01, . . . , 99}. Third, in a mixed scenario (MIX) we assume all firms within NACE classes 01-45 (physical production) have only essential inputs I es i = {01, . . . , 99} and all firms within NACE classes 46-99 (non-physical production) have only non-essential inputs I ne i = {01, . . . , 99}. Fourth, in the generalized Leontief scenario (GL) we assume that for firms within NACE 01-45 the set of essential inputs consists of I es i = {01, . . . , 45} (supplied by physical producers) and the set of non-essential inputs consists of I ne i = {46, . . . , 99} (supplied by non-physical producers), while for firms within classes NACE 46-99 we assume they have only non-essential inputs I ne i = {01, . . . , 99}. We assume that firms within NACE 01-45 have a physical production process while firms within NACE 46-99 predominantly trade or provide services. Note that LIN and LEO provide lower and upper bound scenarios for more realistic situations, where firms have distinct types of production functions as in the MIX and GL scenarios.

Given the details of a production network, W , in terms of weighted in-links, out-links, product classification p i , and production functions f , the economic systemic risk index (ESRI) can finally be computed. It captures the consequences of both, the up-and downstream shock propagation following the default of a particular company, i as outlined in Materials and Methods and Appendix F and Appendix G. The quantity ESRI i can be readily interpreted as the fraction of the production network's observed output that is likely to be affected if firm i fails and its output and demand is not replaced by other firms. The value is computed by simulating the propagation of shocks with the recursive algorithm described in Appendix H.

The directed production network W is reconstructed from fully anonymized VAT micro data of the Hungarian Central Bank. We consider all links, W ij , where at least two recurring supply transactions occurred; see Appendix I and (Borsos and Stancsics, 2020) . Figure 1 (c) shows a section of 4,070 companies and 4,845 links of the empirically reconstructed Hungarian production network, W ij . The section is obtained by sampling 1,500 random nodes and considering all their direct suppliers, yielding 6,113 nodes. Only the giant component (4,070 nodes) is shown. The entire network consists of n = 91, 595 companies; it is too large and dense to visualize its structure in a meaningful way. Already this small section reveals the fact that production by no means happens along independent supply chains, but on a tightly interwoven complex network. Merged chains that are quite extended are clearly visible. The network has an apparent core periphery structure. For visualizations of other aspects of the production network, see Appendix J.

We next calculate the ESRI i for the realistic GL and MIX scenarios as well as for the two limiting cases, LIN and LEO for all of the 91, 595 companies in the Hungarian VAT micro dataset in the year 2017 and for the 85, 131 firms in 2016. Rank-ordered distributions of the ESRI are shown in FIG. 2 (a) in linear-log scale for 2017. For 2016, see FIG. 6 (a) in Appendix K. For the realistic scenario GL (blue), we find that 32 companies show extremely high levels of systemic risk, all being at a value of about 0.23, meaning that about 23% of the entire economy is affected should one of these companies fail and its supply and demand is not replaced. 63 and 165 firms have an ESRI larger than 0.05 and 0.01, respectively. For respective numbers for the other scenarios, see Table I For ranks larger than a characteristic rank of 32, the pronounced plateau in the distribution is followed first by a steep decline and then by a slow decay of the ESRI values. The tail (without plateau and steep decline) can be fitted to a power-law with an exponent of roughly α GL = 0.67 for GL, α MIX = 0.63 for MIX, and α LEO = 0.50. For details of the fits, see Appendix L, FIG. 7 (a). The shape of the rank-ordered ESRI distribution (plateau and power-law tail) is similar to what was found in Fig 4. in (Moran and Bouchaud, 2019) . As expected, the reference case LIN (green) where all companies have linear production functions, generates substantially lower systemic risk levels than GL and MIX. LIN neither shows a plateau nor a power-law decay. The situation is very similar in 2016, see FIG. 7 (b) in Appendix L.

To better understand which companies are forming the plateau of extremely risky companies -data protection regulations prevent us from showing company names or actual turnover -in FIG. 2 (b) we show ESRI i as a function of the strength, s i (firm size within the network). Companies on the plateau belong to industry sectors like energy, manufacturing (electrical equipment, chemicals, computer and electronics, vehicles), or repair of machinery. Red color indicates the highly risky plateau companies; symbol size represents strength, s i . Clearly, the ESRI of plateau firms (located in the shaded area) is not changing with size; in the plateau we find large and small companies (note the range of strength of 4 orders of magnitude), suggesting that firm-size is not able to predict extreme ESRI values at all. For the bulk of companies we find a strong statistically significant correlation of log-ESRI and log-strength (R 2 = 0.90 and slope β reg = 1.12 in log-log regression). However, for individual companies strength is not a reliable predictor of ESRI, since the spread of the ESRI extends up to 4 orders of magnitude. A more detailed regression analysis, found in Appendix M, confirms that generally firm level quantities fail to explain ESRI that is a network-based measure. Note the relation to the Hulten theorem (Hulten, 1978) , stating -in simplified terms-that in an efficient economy the effect of a firm level shock on overall output is proportional to its revenue. Our results are in strong disagreement with this statement, adding further evidence for its limited practical validity (Baqaee and Farhi, 2019) . Finally, we checked that all firms on the plateau are within NACE 01-45, meaning that their production functions are purely of generalized Leontief type with a large share of essential inputs.

The reason for the formation of the plateau is twofold. First, about two thirds of its nodes form a strongly connected component based on highly critical supplier relationships (core of plateau). In this core the failure of one firm leads to the failure of the other members in the component, implying the default of anyone has the same consequences on the entire network (same ESRI). The second reason is that the other third of the plateau nodes are suppliers to the strongly connected component, i.e. their failure causes failure of the nodes in the strongly connected component. These peripheral nodes inherit high ESRI values by supplying (in)directly to inherently risky companies. Note that these two features are impossible to explain with only linear production functions Distributions are rank-ordered, meaning the most risky company is to the very left. The blue line shows the result for the realistic GL scenario (production functions according to industry classification and classifying produced goods as essential or not). A plateau exists around an ESRI of 0.23, containing 32 firms for GL. There is a steep decline to ESRIi ∼ 0.05 from rank 33 to 55. 165 firms have an ESRI larger than 0.01. For comparison the MIX scenario (light blue) is shown (production functions according to industry classification only). As the limiting cases we show the scenarios LIN (green), where all firms have linear production functions and LEO (red), where all firms have Leontief production functions. The tail of the ESRI profile decays as an approximate power-law. (b) ESRI plotted against firm strength (firm-size) in log-log scale. Symbol size represents strength si, red symbols belong to the plateau, emphasised by the shaded area. We find large and small companies in the plateau, suggesting that very high ESRI is not determined by size, even though the correlation of ESRI and strength for the bulk of the companies is high. (c) Systemic risk in 2016 vs. 2017. Colors correspond to ESRI in 2017. Blue and red symbols indicate low and high values in 2017, respectively. Note the strong variability in the plateau firms. There is a significant correlation between ESRI in 2017 and 2016, that indicates relatively small temporal fluctuations for the bulk of the companies. (d) Network of the 32 most systemically risky firms (plateau) in 2017. Node size is proportional to the square root of strength. Link colors correspond to the downstream 'criticality', i.e., the percentage of j's production should i stop functioning, Λ d ij . Red thick (blue, thin) links indicate very large (small) losses of production. Small companies predominately supply to large high-risk companies, thereby inheriting systemic risk. the network spanned by the 32 plateau firms. Node size corresponds to the square root of strength, √ s i ; link colors correspond to the direct down-stream impact of the supplier on the buyer node (criticality), Λ d ij , (for the definition, see Appendix H). Red links indicate highly critical supplier-buyer relations; Λ d ij = 1 means that 100% of j's production fails if i fails. Blue links mark non-critical supplier-buyer relations; Λ d ij = 0.01 means that j's production reduces by 1% if i fails. It is visible that almost all links are critical for the production of the buyer. Possible reasons for these strong fluctuations could be that core plateau firms added additional suppliers for the same input, that supply links were discontinued, or that market share dropped leading to higher supplier replaceability. From the data it is hard to decide how these factors contribute.

To demonstrate that it is necessary to use firm level data for any reasonable assessment of shock propagation in production networks, we compare it to a situation where shocks are studied on a more aggregated level. In FIG. 3 (a) we show three different economic shock scenarios that affect the NACE class 2611 (Manufacture of electronic components). All three initial shocks have the same size on the sector level, but affect different companies within the sector. An initial shock is exogenously applied and triggers a spreading of shocks in the production network. The received shock is the shock each sector receives as a result of the propagation of the initial shock. Shock size is measured as the fraction of the sector's overall strength, (s 2611 = i|pi=2611 s i ), affected by the initial shock. The first shock scenario serves as a reference scenario. It applies a homogeneous initial shock of 18% (reduced production) to all of the 69 firms in sector 2611. The received shock of all (568 NACE 4) sectors (shown on the x-axis) marks the reference response. The received shocks are shown on the y-axis (in log-scale) as a fraction of this reference response, and by construction, the received shock for the homogeneous scenario is 1 (black line). In the second scenario a 100% shock (failure) is applied to a single firm, A, the received shocks in the different sectors (given as a fraction of the reference scenario) are shown as the red line. Sectors are sorted w.r.t. received shock sizes in this scenario (A). The third shock scenario applies a 59% shock to firm B, the responses are shown by the blue line. Firm A is a plateau firm, firm B is not.

Clearly, the specific choice of the defaulting company within the sector has a drastic effect on how the different economic sectors are affected. Note that on the sector level the initial shocks are indistinguishable. In other words, homogeneous sector shocks (aggregated view) yield a grossly false picture of the true shock propagation and thus systemic risks. As so often in networked dynamical systems: details do matter. If studying shocks on the sector level was sufficient, the three scenarios should be strongly correlated. However, the figure suggests that the relative deviations from the sector shock (blue and red line) are negatively correlated with -19% (p = 5 · 10 −6 ). In Appendix O we show the actual received shocks (not their deviations from the reference scenario).

One obvious reason for these large differences between scenarios A and B is made clear in FIG. 3 (b) , where we see for all companies in NACE sector 2611 the inputs at the NACE 4-digit level. Firms (rows) are sorted w.r.t. similarity of their inputs (columns). Even though all companies belong to the same industry sector (NACE 2611) the sectors of their inputs vary substantially. The input sectors of Firm A (red shading) and firm B (blue shading) have only a small overlap, indicated by a Jaccard index of 0.25. Depending on which of two companies fails, different input sectors are affected. The Jaccard index for the customer sectors is even smaller (0.13). For details on customer sectors, see Appendix O. Among the 51 firms (that in NACE 2611 have NACE 4 classified inputs) 56% have no pairwise input sector overlap (Jaccard index of zero) and 85% have no pairwise customer sector overlap. This means that when choosing two firms randomly, the probability of having no common input (customer) sector is 56% (85%). Consequently, shocks to different firms within a sector must lead to different economic sectors being affected.

Based on the reconstruction of all relevant buyersupplier relations between companies in an entire country from VAT data, we demonstrate that the real economy can not be viewed as a collection of separate supply chains, but is a tightly connected directed network that has a strongly (weakly) connected component, containing 26% (94%) of all companies. The network allows us to develop and compute an index for the approximate systemic risk of every individual company to the economy, ESRI. We demonstrate explicitly that systemic risk based on aggregated sector level data yields a severely distorted picture.

We find that the 32 top risky companies contribute to 45 % of the entire systemic risk, the top 100 companies contribute to 74 %. Only 165 companies have more than 1% risk. The 32 top risky companies (0.035% of all 91,595 companies) show extremely high systemic risk of about 23%. Of those, only a fraction appears to be inherently risky (e.g., because of size or high market shares). The average systemic risk of a Hungarian company is ESRI i = 0.00018 (mean). The median is 1.7 · 10 −6 .

Approximately a third of the most risky companies constitute the periphery of the network of the plateau companies. These are small and inherit systemic risk from the core plateau companies because they are crit- . The x-axis denotes the 568 shock receiving sectors, the relative deviation of received shocks from the 18% homogeneous initial shock to sector 2611 is shown on the y-axis. The black line marks the reference scenario. The red line shows the relative deviation if firm A (in 2611) receives a 100% shock; the blue line shows the situation for a 59% shock to firm B. Both firm level shocks are equivalent in size to the 18% sector shock (reference). The particular choice of the defaulting companies lead to drastically different results on how other economic sectors are affected. (b) Most ubiquitous inputs (at NACE 4 level, columns) for 69 firms (rows) in NACE class 2611. Clearly, even though all companies belong to the same sector, their input sector vectors are drastically different. Inputs of firm A (red, 50 distinct inputs) and firm B (blue, 49 distinct inputs) differ substantially (Jaccard index 0.25). Note that 18 firms have no input sectors (empty rows).

ical suppliers to the core. They would not show high ESRI values if they supplied to other -non-risky-companies. This means that several of these smaller companies can be made less risky simply by increasing the number of suppliers of the specific good sold to the inherently risky companies. By analyzing changes of systemic risk from 2016 to 2017, we find that ESRI is relatively stable for most companies. However, several smaller companies change their ESRI from marginal to extreme, and vice versa. The reasons for this might be the risk-inheritance of (non-inherently risky) companies that start (or stop) supplying to inherently risky ones, or changes in market shares that affect their replaceability. We confirmed that to a large extent systemic risk is not predictable with strength (firm-size). The position in the supply network matters much more than company-size, similarly to what has been observed in financial systemic risk in earlier studies (Battiston et al., 2012; Markose et al., 2012; Thurner and Poledna, 2013) .

The presented study has limitations. The proposed measure is an economically motivated, straightforward quantitative measure for systemic relevance of companies in a given production network. It captures up-and down-stream impacts independently. This doesn't cause distortions in tree-like supply chains, however, in networks with strongly connected components it does. In reality, the same firm can be affected by up-and downstream shocks; these shocks can potentially interact. For example, the lack of a critical input (down-stream shock) can translate into an upstream shock for suppliers of complementary inputs. This second-order effect is not considered in the algorithm. Its magnitude needs to be estimated in future work and should be taken care of in more subtle algorithms. Nonetheless, ours is the first attempt to calculate company level ESRI for an entire national economy and the resulting index can be seen as a good first-order estimate for the systemic risk in production networks. We assume, somewhat unrealistically, that all customers and all suppliers are treated as being equally important (proportional rationing). Other rationing mechanisms, for example, prioritizing large customers and suppliers, could lead to modified ESRI estimates. For a detailed discussion of the effects of rationing mechanisms at the sector level, see . Further, for the implementation of the algorithm we assumed a data-driven replaceability index for each firm that is based on its market share (fraction of output within its NACE 4 class). This implies that firms with high (low) market shares are difficult (easy) to replace, see Appendix E. Further improvement of the ESRI would be possible by modelling the substitutability of inputs and the replaceability of suppliers in more detail. However, this comes with the cost of more parameters to calibrate, more detailed data needs, and significantly higher computational effort.

For assigning production functions to companies, we used a simple matching based on NACE 2-digit categories. Although we think that this is a good first approximation, we emphasize that this can be a source of error. Given the current data availability, an exact and objective mapping of production functions to companies is not feasible. Large-scale firm level surveys providing information on how strongly firms are affected by the lack of specific inputs could be a useful first step clarifying this point.

We assume a static production network, which in reality, is a temporal network. With this simplification, seasonality effects in production and manufacturing are missed. Also, in the present implementation we don't consider competition between companies that would result in a dynamic restructuring of the network. Moreover, it is technically hard to interpret the time index t in terms of how long shocks really need to spread and converge. For simplicity, we assume uniform spreading, i.e. all periods t are of the same abstract length. Imports, exports, and production networks of other countries are not considered due to lack of data. This is a limitation since Hungary is a small, open economy with significant exposures to shocks in the global economy. In principle the vulnerability to initial shocks from import and export relationships can be considered in our framework by using coarse grained sector level import-export information.

Despite these limitations, our findings have a number of policy implications.

Monitoring. The economic systemic risk index for individual companies, ESRI, allows countries to use VAT tax data to identify critical companies in the economy and their critical supply relations (via their marginal ESRI contributions, see (Poledna and Thurner, 2016) ). Our findings indicate that only a few firms pose a substantial risk to the overall economy. These should be monitored closely if they produce goods of societal relevance. In particular, sharp increases in firms' systemic riskiness can be monitored over time. Several countries are currently implementing "supply chain due diligence" laws, for example, Germany (Financial Times, 2021c) or the USA (Biden, 2021) . The possibilities for systemic risk monitoring as presented in this paper should be kept in mind, when designing new regulation.

Systemic risk mitigation. A straightforward way to increase resilience in the real economy is to introduce supply chain redundancies, where risk is reduced by altering the network structure and thereby reducing the probability for inherently systemically risky firms to suffer failures that result from a lack of inputs. We have shown that a simple network analysis allows us to identify the inherently risky firms. For those, contingency plans should be established to buy time for changing suppliers or to build alternative production capacities. The change of ESRI over time would allow countries to monitor if implemented policies really increase resilience levels of the economy.

Avoiding risk concentration. In (Poledna and Thurner, 2016) it has been shown how incentive schemes can be designed to make financial networks safer without making them less efficient. A similar scheme could be devised for production networks by de-incentivizing critical supply relations. A simple scheme would be to have at least one back-up supplier for every critical product that is shipped to inherently risky firms. Like in regulations of the financial sector, risk concentration can be avoided by demanding that no customer should have more than a certain supply exposure (e.g. 10% of total exposure).

Inventory buffers. A natural possibility is to introduce mandatory inventory buffers for systemically risky companies, to ensure production in situations where a set of critical suppliers default.

Make economic systemic risk visible for firms. Today, companies usually manage their direct suppliers. Most companies do not know their higher-order dependencies in the production network (GEODIS, 2017) . It is conceivable that it could be highly beneficial to many producing companies to obtain a better overview on their actual supply and customer risks. If countries provided a more global view on supply chains to identify critical situations, this could lead to deeper and more proactive systemic supply chain management, that would increase overall economic resilience by using market mechanisms.

Policies and regulatory measures of this kind might run contrary to the hunt for efficiency that dominated the last decades, since supplier relations are time-and cost intensive. More resilience does not pay off in the short term. It would be a fascinating question to see to what extent the economy could be made more resilient by not making it less efficient, i.e. to maximize the resilience-gain per link. For financial networks the potential of designing networks with lower systemic risk, but comparable economic function was recently highlighted in (Diem et al., 2020; . A related question is whether service-based economies are more resilient than production-based ones because services typically require less critical suppliers (due to their linear production functions). Instead of asking whether large economies are more stable than small ones (Moran and Bouchaud, 2019), our methodology, given data from other countries, could contribute to the question, whether production-based economies or economies with high levels of input-redundancy are more resilient.

The economic systemic risk index for firm i is calculated as

s out j n l=1 s out l the out-links by updating

and upstream along the in-links by updating On the industry level studies on the propagation of shocks in production networks exist at least since the famous Leontief Input-Output analysis (Leontief, 1928 (Leontief, , 1991 Miller and Blair, 2009 ). In (Bak et al., 1993) it is shown theoretically that small shocks to individual sectors can have effects on aggregate output. Refs. (Acemoglu et al., 2012; Carvalho, 2014; Gabaix, 2011) confirm this finding. (Hallegatte, 2008) assess higher order shocks from the Hurricane Katrina disaster. More recently, investigated how demand and supply constraints on the sector level affect GDP. In (Colon et al., 2021) it is modeled how shocks spread through the Tanzanian supply and transport network. The sector level analysis produced valuable insights, but it has certain limitations due to the nature of the underlying data. Only recently high-quality large-scale firm level data has become available that allows for new insights.

The methodology presented in this paper addresses three usual shortcomings of sector level analyses. First, the data shows that even within fine grained industry classifications (NACE 4) firms tend to have considerably heterogeneous input sectors and customer sectors. In fact when comparing the pairwise input (customer) sectors of firms we see that for all firms in NACE 2611 (Manufacture of electronic components) 56% (85%) have no overlap of input (customer) sectors even though they belong to the same NACE 4 category, see FIG. 3 in the main text and Appendix O. Second, this intra-sector heterogeneity can lead to inaccurate results when using them for assessing shock propagation in production networks. Especially if the initial crises scenario does not affect all firms within a sector to the same extent (for example in the current COVID-19 crises). Our proposed framework takes this into account and yields different cascades for shocks which would appear to be the same at industry level, but are actually distributed heterogeneously among firms within an industry. FIG. 3 (a) shows in a simple example how two cascades that are the same on the industry level lead to very different impacts on other firms and industry sectors. Third, in contrast to sector level models, each firm in our approach has a specific production function based on its industry classification that is calibrated to the observed individual input vector of the respective firm; see Appendix D for more information on the calibration of production functions. This matters especially since we show in FIG. 3 (b) that firm input vectors, even in the fine grained industry classifications, vary on the sector level. The same is shown for customer sectors, see Appendix O, FIG. 9 (b) .

Since, data has become available also firm level analysis has been performed. Ref. (Atalay et al., 2011) uses a generative model for firm level production networks that matches degree distributions better than scale-free frameworks. Belgian VAT data has been used to study productivity shocks to individual firms and their effects on aggregate output with a computeable equilibrium model (Magerman et al., 2016) . The theoretical study (Moran and Bouchaud, 2019) shows that shocks propagate widely as soon as production function have a "Leontief" component, but not when they are of pure Cobb-Douglas type. Probably the closest study to ours, (Inoue and Todo, 2019) , investigates how firm level shock propagation in response to an initial shock -the great earthquake in Japan in 2011-based on an estimate of the Japanese production network. However, they focus on effects on aggregate output and do not compute firm level systemic risk.

In the area of financial networks systemic risk has been extensively studied for about two decades (Boss et al., 2004a; Freixas et al., 2000) . The importance of being able to measure systemic risk of single firms (banks, insurance, funds) became apparent in the 2008 financial crises and the European government debt crises 2012. Initial systemic risk assessment methodologies have been shown by, for example, (Boss et al., 2004b; Eisenberg and Noe, 2001; Furfine, 2003) they have been refined to macroprudential stresstesting models that can assess the effects of adverse macro-economic scenarios (Elsinger et al., 2006) and to measure the impact single banks have on the entire network (Battiston et al., 2012; Cont et al., 2010) . A wide range of studies revealed various properties of financial networks that improve the understanding of how systemic risk emerges, for example, the role of network topology (Gai and Kapadia, 2010; Nier et al., 2007) , the role of diversification of external assets (Beale et al., 2011) , or the role played by large banks (Arinaminpathy et al., 2012) . These aspects are crucial for the management of systemic risk. The design of regulatory policies to address the adverse economic and societal effects of systemic risk has strongly benefited from these academic contributions. However, there are important differences between financial and production networks and how stress is spreading there, that require adaptions to how systemic risk and the spreading of shocks is defined and modelled. For example, in production networks in-and out-links affect heterogeneous production process in contrast to stock quantities like equity or liquidity buffers in financial networks. Further, production networks tend to be orders of magnitude larger than banking networks. Here we extent the ideas of measuring systemic risk in financial networks (Bardoscia et al., 2015; Battiston et al., 2012; Cont et al., 2010) to a more general framework of companies in a production network at a national scale. To estimate it we use micro level VAT data of Hungary (Borsos and Stancsics, 2020 ).

Around the world there are millions of different products produced, all having small distinctive features. However, to allow for a meaningful statistical and economic analysis these products are usually grouped and categorized into product classes according to some common features. Examples are the Cooperative Patent Classification (CPC) and the EU's classification of products by activity (CPA). They contain 2647 and 3142 classes, respectively. The second kind of classification aims at grouping companies that produce these product groups based on common activities. Industry classifications have various levels of granularity, for example the International Standard Classification of All Economic Activities, ISIC, has 88 categories at the 2-digit and 419 at the 4-digit level; the Statistical Classification of Economic Activities in the European Community NACE has 88 categories at the 2-digit and 615 4 digit level; the North American Industry Classification System, NAIC, features 1057 categories at the 6-digit level. Note that the CPA and NACE classifications can be mapped onto each other. A concise overview can be found in (EUROSTAT, 2021) .

In practice however, these classifications are not linked to the single supply transactions, W ij . As outlined in the introduction of the main-text economists typically resort to the assumption that each company i produces one of m possible different products, identified by the firm's industry classification. In mathematical notation we use the industry affiliation vector, p, to identify the industry affiliation p i ∈ {1, 2, . . . , m} for each firm i in its entries. On the modelling side, this simplifying assumptions amounts to reducing a firm level input-output matrix, with each column determining the necessary inputs for each product firm i produces, to a single input vector. In mathematical terms this simplification reduces the output vector x i1 , x i2 , . . . , x im of company i to a scalar output x i of type p i . Consequentially, the input matrix (Π i -with element Π i kl determining the amount of input, k, used to produce output x il -reduces to an input vector, (Π i1 , Π i2 , . . . Π im ) and the production function which in reality is a map

In short, a multi-layer network is aggregated to a single network layer, where each buyer-supplier relation,

, is simplified to a single transaction, W ij (t), of product, p i . The reduction of a distinct input and customer vector for each product a firm produces to a single vector for each firm leads to a blur of the up-and downstream contagion. This is illustrated with the case where firm i produces two different products, A i and B i , requiring each two types of inputs. A i requires α and γ, B i requires β and γ. If either the suppliers of input α or input β fail to deliver only the production of one of the two products A i or B i is adversely affected. Consequently, a customer that only buys A i is affected only in case input α is becoming scant, while another customer who only buys B i is affected in the case where input β is not available. In the aggregated picture where there is only one abstract product of type p i , both customers are affected by a loss of input α or β. Note in the case of Leontief production functions, if either input α or β is not available, this would affect the entire production of firm i, even though in reality it would only affect the production of either A i or B i . Similarly, for upstream contagion, if only one of the two customers stops buying, just the suppliers of either α or β are affected but not both. The presented framework can be easily extended to account for this as soon as appropriate data becomes available.

The mere existence of a buyer-supplier connection does not say much about how their production processes depend on each other, i.e. how the failure of one affects the other. This exposure depends primarily on the type of goods and services delivered and how these specific intermediate products (inputs) are used in combination with each other. Naturally, some inputs are more critical for specific production processes than others. We need to know the type of goods and services that are delivered and an approximations for the firms' production processes to model these exposures. As outlined in the main text, production processes of firms in the economics literature are commonly modelled with production functions. The standard choices are the constant elasticity of substitution production function (CES) (McFadden, 1963) and its special cases, the Cobb-Douglas and the Leontief production function. Note the linear production function is also a special case of CES.

Intuition. In FIG. 4 we show a schematic visualisation of the three special cases. Figure 4 (a) shows a Leontief production function for firm i, with two types of input, Π i1 (screws) on the x-axis, Π i2 (wood) on the y-axis, and the resulting production level, x i (tables) on the z-axis,

The amount of tables produced depends on the availability of both products according to a fixed production recipe, α i1 screws and α i2 cubic meters of wood. Since, we only know the volumes (= price × quantity), the interpretation 

, is constrained by the least available input. We illustrate this with a carpenter building tables. She requires an exact amount of wood and nails for every table. (b) Linear production function of two inputs, xi = Π i1 α i1 + Π i2 α i2 . For a linear production function the output is proportional to a linear combination of the inputs, here illustrated by a retail company, selling shirts and trousers. If shirts are unavailable as an input, the revenue generated by trousers is not affected. (c) Cobb-Douglas production function, xi = βiΠ α i1 i1 · Π α i2 i2 , with α1 + α2 = 1. The Cobb-Douglas production function can be interpreted as an intermediate case, allowing for partial substitution. Exemplified by the production of crops, the amount of land used for the crops can be reduced if e.g. more fertilizer is invested. However, with the complete lack of either production factors, any production is impossible.

changes such that α i1 is the money spent on screws as a fraction of the table's value and α i2 is the money spent on wood as a fraction of the table's value. In this case, if 50% less of one of the inputs is available, only 50% of the tables are produced. In industrial processes the production recipes encoded in the α ik are commonly used for production planning and management. Figure 4 (b) shows a linear production function of firm i with two different types of inputs, Π i1 (trousers) on the x-axis and Π i2 (jackets) on the y-axis and the resulting production level, x i (sales of trousers and jackets) on the z-axis,

It is clear that "production" in this case is possible with only one of the inputs. The coefficients α i1 and α i2 determine how important each of the inputs is for the output. The ratio, α i1 /α i2 , determines how much of input 1 needs to be added to substitute the one unit loss of input 2 to keep the same level of production. In our example the inputs are equally important. The linear production function can be used to reasonably describe a trades business or services.

In the case of a trade business the coefficients α ik correspond to the markup firm i charges when re-selling input k to customers. Similarly, for the delivery of services, inputs are not used jointly as intermediate inputs to produce a new physical product, but to support the creation process. Illustrative examples are a hairdresser, running out of shampoo can still cut hair, a consulting firm running out of office supplies that still can produce powerpoint slides, a software developing company facing delays of new computers that can still produce new code on the old machines. Even though the processes are not becoming impossible to execute, they become less efficient or reduce in quality, thus leading to effectively less output. Text books usually state that linear production functions are of limited or no practical use (Varian, 2014) . However, it is obviously a sensible approximation for the above mentioned situations. Clearly, it fails to describe the production process of physical goods where inputs are inherently combined and transformed to become an output of a new type. Figure 4 (c) shows a Cobb-Douglas production function for firm i with two different types of input, Π i1 on the x-axis, Π i2 on the y-axis, and the resulting production level, x i on the z-axis,

It is clear that "production" is more efficient if both inputs are available, because the fewer of input one is available, the more of input two is needed to compensate for this loss and keeping the production level constant. In the extreme, an infinite amount of input one would be needed to completely substitute two. This is of course not possible in reality, but indicates that production with one input alone is not possible. In its original formulation it was not intended to describe production based on different intermediate inputs, but on labour (l i ) and capital Douglas, 1976) . Note that if the log is taken on both sides of the equation, this yields a linear production function. The interpretation is that if Π i1 grows (decreases) by 1%, output will grow (decrease) by α i1 %, at the given production level. On the sector level α ik can be set to the technical coefficients, derived from input-output tables (Carvalho and Tahbaz-Salehi, 2019; Miller and Blair, 2009) , i.e., α ik = Π ik /x i . This is equally possible on the firm level.

Interpretation in the context of shock spreading. Usually the production function specifies how much goods and services (output), x i , of type p i ∈ {1, 2, . . . , m} firm i can produce with a given amount of m different intermediate products (inputs) (Π i1 , Π i2 , . . . Π im ) ∈ R m + , manufacturing equipment (capital) c i , and its employees (labour), l i . The production functions are usually used in a long-term perspective, where different inputs (traditionally labour and capital) can be substituted, by for example, using more machines and less labour. In practical terms this implies a change of a firm's production process. Such changes are usually not possible in the short-term and often need at least a few months (lead times for buying new machinery, designing the new business processes, hiring differently qualified employees, train employees on new machines or software, etc.). This is especially true for complicated production processes of complex goods. Thus, for short-term spreading of shocks caused by, for example, the failure of a supplier, we assume that the production processes are constant and we neglect capital and labour within the production function. Similarly, we assume that intermediate inputs can not be substituted in the short-term by simply buying more of a different input. In the long-term, however, certainly, missing inputs can be substituted by other inputs when the production process is sufficiently adapted. Consequently, for the analysis of shock propagation in the production network, we use the production functions to determine how much a firm can still produce in the short-term if a certain input is not available. As discussed in Appendix E, the assumption of replacing a failing supplier by another that can deliver a comparable input, has a strong impact of how shocks are spreading.

Leontief-, linear-, and generalized Leontief production functions. We now show the exact specifications of the production functions used to model the spreading of up-and downstream shocks in the economic systemic risk index, ESRI. Note that the calibrations are made for an observed production network, W , and an industry affiliation vector, p. These two objects determine the output volume of every firm and the volume of respective inputs purchased in the observed production network. We specifically use the linear production function for trades businesses and services companies, and a modified Leontief production function for firms with a physical production processes, to allow for a different treatment of inputs from physical production processes and from services.

We start from observed outputs and inputs of firms determined by the supply network, W , where W ij denotes the volume (value) supplied from firm i to firm j in the respective observation period. The input, k, of firm i is mapped to its suppliers by,

where δ a,b is the Kronecker delta that is equal to 1, if a = b and 0, otherwise. The output, x i = l=1 W il , is the sum over all supply transactions to other firms, l. First, we consider the Leontief production function for company i as

Note that this formulation implies that all suppliers j of i supplying the same input, p j = k (industry classification) are treated as perfect substitutes. The corresponding parameters α k for company i are the technical coefficients,

α ik determines the fraction of output that firm i needs to spend on input k to produce this respective output. Second, we consider the simplified linear production function

with

It is clear that Eq. (D4) is a special case of Eq. (D2) when all product types are the same. For shock propagation this implies that the loss of an input affects the output proportional to the value of the inputs. Note that the linear production function has a natural interpretation for a trades business, where 1/α i − 1 corresponds to the definition of the markup (M arkup = (Sale P rice − Cost)/Cost). In our 2017 dataset, trades businesses (NACE classes 46 and 47) make up for 17% of firms. 2 . Third, we consider a generalized Leontief production function (GL) where only essential inputs are treated in the Leontief sense and non-essential inputs are treated in the linear sense. We denote the set, I es i , that contains all essential input types k, entering production in a Leontief way, and the set I ne i that contains all non-essential input types k entering in a linear way. Then, we define the modified Leontief production function as

The parameter β i is defined as the production level that is attainable with only essential inputs k ∈ I es i , i.e.,

As for the linear production function, the parameter α i is the fraction of output spent on the value of inputs. It ensures that a lack of non-essential inputs interpolates between the full production level x i and β i . It is clear that both the Leontief and the linear production function are special cases when either all inputs are essential or non-essential. Note that here non-essential means that production is not stopped when the input is lacking, but there is still an impact on output, i.e. the input is non-essential yet relevant. This formulation can be easily extended to cover the additional case that inputs are non-essential and non-relevant. In that case a third group of inputs needs to be specified and those inputs simply do not affect the production function. We leave this for future research.

On shorter time horizons the mode of production can only be changed in relatively subtle ways, e.g., by using a given input from a different supplier and even that is sometimes not possible if goods are not highly standardized. However, major modifications to the production process itself are much harder, especially for sophisticated products. Detailed modelling of how different companies can be replaced by others (as suppliers) would require the detailed knowledge of inventory levels, production capacities of suppliers, and detailed product information (e.g. if two suppliers in the same industry can produce the same good). With such additional data available, a consistent modelling of supplier replacement could be implemented with a dynamic rewiring of the network, i.e., companies looking for new suppliers with idle production capacities or inventory to replace a failed supplier. However, it is difficult to model this behaviour in a realistic way, without many additional parameters. Such extended modelling efforts would be computationally significantly more expensive.

Only a few strategies in the firm level supply chain contagion literature tackle this challenge, however, generally all attempts suffer from limited data availability. Ref. (Inoue and Todo, 2019) does consider replaceability only between existing suppliers. (Wu, 2016 ) uses a different strategy and creates a measure of replaceability based on the weight of a given supplier in its costumers' production-related cost. This measure assumes that if a supplier is an important part of a firm's production, then it is harder to find a substitute for it. Another, frequently used solution, is to distinguish between standardized goods (goods with a clear reference price listed in trade publications) and differentiated goods (goods with multidimensional characteristics) based on (Rauch, 1999) . Although this distinction can be used as a categorical variable in econometric estimations, e.g., (Giannetti et al., 2011) it is not informative regarding the extent of replaceability even for standardized products. Ref. (Barrot and Sauvagnat, 2016) uses two other proxies to measure

W ji δ 2,p j + · · · + 1 α im n j=1 W ji δm,p j ) .

Then, for trades businesses the parameters α ik would be chosen such that 1/α ik − 1 corresponds to the markup of product group k for company i. the specificity of suppliers: the level of R&D expenditures and the number of patents held by a firm. Unfortunately, these pieces of information are only relevant for a tiny fraction of companies and not at all applicable to the entire network of Hungarian firms. We propose a different strategy that employs a straightforward, data driven way to construct a short-term supplier replaceability index based on intra-industry market shares. The basic intuition is that a supplier, having a small market share within its industry, on average should be relatively easy to replace by a small increase of the production of its competitors to cover the additional demand for their products. However, a supplier producing a considerable share of the goods in a given industry is more difficult to replace, as it is unlikely that its competitors can increase their production immediately. Having a large market share within an industry category in a country does of course not necessarily mean that a firm is large. It also needs to be taken into account that potential alternative suppliers in the given industry might also have experienced shocks during the contagion processes in the model. To make this approach more realistic, we take the deterioration in their production capacity into account. Here, it is important to distinguish between the two different sources that could cause reductions in the production level of these potential alternative suppliers. On the one hand, one should account for downstream shocks, i.e. shocks coming from the suppliers, which is a truly limiting disruption in their production. On the other hand, one should disregard the upstream shocks experienced by them, because demand shocks coming from customers are not actual restraints on their production (at least from the point of view of replacing reduced output of their competitors). We can also account for the fact that in the case of a system-wide crisis it might not be possible to find alternative suppliers; see Eq. (G14) for the mathematical formulation of the supplier replaceability.

The interpretation of the replaceability of suppliers in our approach can be illustrated by the following example. Assume that a supplier with a 10% market share within its industry (after considering also its competitors' states) is responsible for 50% of the input required by one of its customers in a given product category. If the production level of this supplier drops to 80%, then the production of the customer will decrease by 1% (= 10% × 50% × (1 − 80%)). If we disregarded the possibility of replacing this supplier, the corresponding decrease in the firm's output would be 10%. This way, we are able to replace missing supplies w.r.t. the market conditions; i.e. our replaceability factor reflects not only the fact that the given input in this example can be bought from the remaining 90% of the market, but also acknowledge that this replacement might not be possible entirely, or it entails some costs. In future research the interaction of upstream and downstream shocks needs to be taken into account also for modelling the replaceability. This is important because, if some producers receive upstream shocks, they would have an additional capacity to supply other firms with this, which themselves might be suffering a downstream shock. If these idle supply and idle demand is matched, this would result in less contagion.

We quantify the systemic risk of a (temporally) failing firm to the overall production network -and hence economy and society (e.g. through employment problems) -by the hypothetical direct and indirect up-and downstream effects of this initial failure on the remaining production network. The basic philosophy of our systemic risk measureconsidering the up-and downstream dependencies of firms -is captured in the following recursion scheme. First, we define the observed production network (for a given year) as our initial state, W (0), and calibrate the chosen production functions to this state t = 0, as shown in Appendix D. Second, we recursively simulate how the production levels of the other firms in the network are affected in response to the initial shock (failure of a firm). Third, we observe the production levels of all firms after the recursion reaches a new stable state (shocks stop propagating). We define the economic systemic risk index, ESRI i , of company i as the fraction of the production networks' overall output that is (temporally) impeded in response to the firm's (temporal) failure. The downstream recursion is derived from linking the amount firm i can produce at t + 1 with the available inputs at t,

to the production network, W (t). The amount Π ik (t) firm i uses from input type k in period t can be mapped to its in-links (suppliers), W ij (t) by Π ik (t) = n j=1 W ji (t)δ pj ,k , where δ a,b is the Kronecker delta. Similarly, the upstream recursion is derived from the dependence of firm i's output on the available demand for its products. The demand is mapped to the out-links (buyers), W il , via

The in-links of node i, W ji (t), themselves depend on the production x d j of suppliers j ∈ {1, 2, . . . , n} and subsequently on their in-links (supply). The out-links, W il , depend on the buyers l ∈ {1, 2, . . . , n} production x u l and subsequently their own out-links (demand). Consequently, the production of each firm recursively depends on the other firms' production. To simulate the impacts of firm j's default, we assume an initially stable state at t = 0, corresponding to the observed production network, W (0). The initial output levels are given by the initial state W (0) and x d (0) = x u (0) = x(0) and x i (0) = n l=1 W il (0). We apply an initial shock to the network by assuming the (temporary) failure of firm j at t = 1. We set x d j (1) = x u j (1) = 0, and consequently, the out-links W ji (1) = 0, for all i and the in-links W lj (1) = 0, for all l. We can simulate the shock propagation by updating the production levels of all other firms. We continue updating the production levels iteratively until production levels of firms do not change anymore.

The relative production level at time t depending on the received up-and downstream shocks is given by

respectively. We assume a proportional rationing mechanism, by setting

. This implies that if a firm received a downstream shock it will forward it to it's customers pro rata according to their initial percentage shares. Similarly, if a firm received an upstream shock it will forward it to it's customers pro rata according to their initial percentage shares. Note that this assumption also implies a static production network W . For an overview of the effects of rationing mechanisms on the sector level, see . The amount that firm i can produce at t + 1, given the suppliers' production levels, h d j (t) at t is

The amount firm i produces at t + 1, given the current production level of its buyers, h u l (t) at t is

Note that we assume here that companies do not produce simply on stock if there is no demand. The failure of firm j, h u j (1) = h d j (1) = 0 is propagated through the network by iterating the Eq. (F4) and (F5) until the network reaches a stable state at time T . The stable state is reached at time

where = 10 −2 is chosen as a convergence threshold. Thus, we assume that once all shocks are smaller than , the propagation stops -the production level of the firm is no longer adapted, when shocks are becoming smalland the corresponding time point of convergence is T . The final production level of every firm i is set to h i (T ) = min(h d i (T ), h u i (T )). The ESRI of firm j is now calculated by weighting the relative losses in production of each firm i by its out-strength, s out i ,

The calibration of the recursion algorithm to the specific production functions used in the empirical analysis is shown in Appendix H; the derivation of the recursion equations is shown in the following Appendix G.

We explicitly determine the relations for Eqs. (F4-F5) for the generalized Leontief, the Leontief, and the linear production functions. We use the term update equation and recursion interchangeably.

Downstream recursion. We derive the downstream update equation for the generalized Leontief production function from Eq. (D7)

It is convenient to work with a direct recursion with h d i (t + 1) on the left hand side and h d j (t) on the right hand side. Thus, we divide Eq. (G1) on both sides by

After a few simplifications (see Eq. (G18-G20)) we can write more compactly

The elements of the matrix Λ d1 are defined as

andβ i is simply the relative fraction of production possible with only essential inputs k ∈ I es i . Note that the Leontief production function (Eq. D2) and the linear production function (Eq. D4) are special cases, where either only the first, or the second part of Eq. (G3) is present.

The elements of the Leontief downstream impact matrix, Λ d1 ji , capture the impact that the failure of firm j has on firm i, given that j is supplying a "Leontief input" (essential) to firm i. Note that Λ d1 ji = 1 means that firm j is the only supplier of goods of type p j to firm i. Similarly, the elements of the linear downstream impact matrix, Λ d2 ji , capture the impact the failure of firm j has on firm i if j is supplying a "linear input" (non-essential) to firm i. Note that the case Λ d2 ji = 1 is only possible if firm j is the only supplier of firm i. This difference gives already a glimpse on the fact that the Leontief production function leads to much higher contagion levels than the linear one. Further, the case of Λ d1 ji = 1 becomes more frequent if the product types (industry classifications) become more fine grained. This leads to the intuitive behaviour that companies producing scarce and essential (to many firms) resources, are expected to have a high systemic risk index.

As will become visible in Appendix H, for the ease of numerical implementation, we split the updating of Eq. (G3) into two parts. First, we define the unified downstream impact matrix

where k ∈ I es i contains all inputs k that are essential to firm i and k ∈ I ne i contains all inputs k that are non-essential to firm i. Again, firms having a pure Leontief or linear production function are special cases. Second, we define the relative share of input k, available to firm i at time t as

Consequently, the n × m matrix,Π(t), contains in row i the relative input vector at t for firm i. Relative refers to the initial state t = 0. Hence,Π(t) can be updated by a matrix multiplication, (Λ d ) P , where P ∈ {0, 1} n×m is defined as

Note that -due to the definition of Λ d -for the inputs k ∈ I es i the relative input,Π ik (t), is equal to 1 if no supplier defaulted, while for k ∈ I ne i it is equal to the share of input k's value out of the value of all inputs, i.e., its only exactly 1, if i buys only one single non-essential input, k. Third, we can update the variable h d i by evaluating the relative production function for each firm i h d i (t + 1) = min min

Note that the quantity k∈I ne iΠ ik (t) is 0 if all inputs k ∈ I ne i are not available. In the case of a pure Leontief producer this is always zero andβ i = 1, while for the case of a pure linear producer this quantity is one at time t = 0 andβ i = 0.

Upstream recursion. We follow the same procedure for the upstream update equation based on Eq. (F5). The production of company i in response to demand reductions from its customers is

and we divide both sides by x i (0) to get

The upstream impact matrix is defined as

and its elements, Λ u ji , determine the impact of the failure of buyer j on supplier i. Λ u ji = 1 only occurs if j is the only buyer of i. Note that in comparison to the definition of the downstream update equations, we assume that upstream shocks are independent on the production functions. We assume implicitly that firms keep the proportion of their inputs fixed when output is reduced or increased. Therefore, we ignore that some inputs are "fixed costs" and will or cannot be affected by upstream contagion. Further, we can not take into account contractual obligations, which enforce a supply transaction also in case the demand for the customers product is reduced.

Incorporating exogenous shocks. For analysing an exogenous shock we have to make the initial exogenous shock to every firm explicit in the updating equations Eq. (G10) and Eq. (G9). Simply iterating these equations after the initial failure of firm j, x d j (1) = x u j (1) = 0 =⇒ h d j (1) = h u j (1) = 0, leads to an instant recovery of the initially failed firm because all its suppliers and customers are not affected (yet). 3 In our model the initial failure is an abstract exogenous shock to firm j that could constitute everything from the destruction of the business premises of a firm, by for example a fire, a government mandated closure in a pandemic, to a strike. Technically, such exogenous shocks affect either capital, labour, or both in the production function. As mentioned in the main text, we don't model capital and labour explicitly but simply define an exogenous constraint, ψ i ∈ [0, 1], and set it to the fraction of production that is still possible after the initial shock. For example, ψ i = 0.8 means that 80% of production is still possible after the initial shock. Both shocks to labour or capital can be implemented with the variable ψ i . Equations (G10) and (G9) change in the following ways

Note that instead of introducing the parameter, ψ i , one could also work with incremental updates, h d i (t) − h d i (t − 1), as e.g. in (Bardoscia et al., 2015) . When firms receive an exogenous shock, (1 − ψ i ) < 1 (not a 100% failure) upstream and downstream shocks are propagated on top of the initial shock 1 − ψ i . However, when including ψ i directly into upstream and downstream updates Eqs. (G12-G13), the additionally received upstream and downstream shocks need to be larger than the initial shock, 1 − ψ i , so that firm i propagates them further. Note that a positive shock can be considered if ψ i > 1. As in (Inoue and Todo, 2019) the initial shock could be made time dependent ψ i (t) to also model potential recovery from the initial shock.

Modelling supplier replaceability. So far the replaceability of suppliers was ignored. It is conceivable that in practice it is totally unrealistic for many firms (suppliers) to be not replaceable at all, even on the short-term; see also Appendix E. A straight forward proxy for a firm's replaceability -that can be inferred from available datais its market share. Intuitively, a company with a small market share (within a fine grained industry classification) will be easier to replace than one with a large market share. For example, the lack of inputs caused by a (temporary) failure of a firm with a 5% market share can be most likely compensated to a large degree by its competitors, as long as the produced goods are reasonably standardized. This degree of standardisation varies significantly and is ignored in our approach. On the other hand a temporary failure of a firm with a 50% market share is much more difficult to compensate by competitors, given that there are inventory and production capacity constraints.

Modelling the replaceability of suppliers explicitly would require that a firm, whose supplier defaulted creates a new link with one or more firms of the same industry classification. This link creation depends on the available inventory level, geographical proximity, and available production capacity. To be meaningful, such a network rewiring model needs to be calibrated to appropriate data that is hard to access. Moreover, it becomes computationally significantly more expensive and involved. We therefore choose a simpler approach for calculating the systemic risk index for every firm. We model the replaceability based on firms' market shares with a simple linear factor. In particular, we calculate the replaceability factor at time t

Over time, σ i (t) can only increase and σ j (t) ∈ [0, 1]. It is bounded from below by firms with zero market share. When other suppliers of the same product have reduced outputs too, the supplier becomes harder to replace. Note that since the supply of firm i itself is also part of the total market (denominator), a market share higher than 50% can not be replaced anymore by competitors. This would imply a short-term doubling of their capacities. Note that Λ d ji 1 − h d j (t) is the shock that i receives from j and since σ j (t) ∈ [0, 1], small market shares dampen this shock substantially.

To integrate this factor into the recursion we have to adapt the update equation for the availability of relative inputsΠ ik (t). For k ∈ I es i we now updatẽ

and for k ∈ I ne i we define an artificial product category k that is updated according tõ

Note that if all σ j (t) = 1 and all h d j (t) = 0 for all non-essential input suppliers (p j ∈ I ne i ) the right hand side of Eq. (G16) is equal toβ i and we can write the downstream update more compactly as

The matrix P can be adjusted accordingly.

Simplifications. For completeness, we specify the simplification steps for the first term in Eq. (G2)

and the second term in Eq. (G2)

We denote the relative fraction of production that possible with only essential inputs k ∈ I es i as

(G20)

We finally show the necessary equations to calculate the ESRI for each firm in the production network, given the network W , the industry affiliation vector p and the sets of essential products I es i and the set of non-essential products I ne i . The downstream impact matrix Λ d is computed as

The elements of Λ d1 and Λ d2 are defined as

The upstream impact matrix is calculated as

Next, the initial exogenous shock parameter is set to ψ i = 0 for the failing firm i (for which we compute ESRI i ) and ψ i = 1 for all other firms i. Then, we iterate the following set of equations. First, for each firm the market share at time t is calculated as

Next we calculate the relative amount of essential and non-essential inputs available for each firm i. For essential inputs k ∈ I es i we update the relative share available according tõ

For the non-essential inputs, k ∈ I ne i , we update the relative share available of the auxiliary product category k according tõ

Finally, we update for each firm i its relative production level determined by the received downstream shocks

and its relative production level, determined by the received upstream shocks

The update equations are iterated until the algorithm converges at time

where = 10 −2 is chosen as the convergence threshold. Once all shocks are smaller than the propagation stops and the corresponding time point of convergence is T . The ESRI of firm i is now calculated as the weighted sum of lost output in the production network in response to the (temporary) failure of firm i

Note that the interpretation is conditional on the assumption that neither the supply of inputs nor the demand of the failing firm are replaced. Thus, the ESRI should be interpreted more as an economically motivated metric to rank firms with respect to their systemic importance than as a precise estimate of the actually lost total output. Note that the Hungarian production network does not compromise all sales and purchases of firms based in Hungary, due to the reporting threshold and international transactions. When calculating the up-and downstream impact matrices, Λ d and Λ u , we over-estimate the impacts single links have. The upstream impact from j to i, Λ u ji , associated to link W ij needs to be adjusted for the sales transactions that are not included in the observed network, W . These residual sales transactions are included in the revenue of firm i. Thus, for the empirical application of the algorithm we re-weight the element Λ u ij by the factor s out j /r j (0), where the numerator is the output recorded in the observed production network and r j (0) is the revenue of firm j for the respective year. Similarly, we re-weight Λ d ji by s in j /c j (0), where the numerator is the volume of purchased inputs within the observed production network and c j (0) denotes the material costs in the respective year. We use income statement data of the years 2017 and 2016 of the firms, wherever available to make this adjustment. For those firms where this is not available we apply a factor of 1. We report that ( Omitting this re-weighting would lead to a significant overestimation of ESRI. A potential error arises by re-weighting for firm i Λ d ji by the same factor for all input suppliers j, since we do not know which input types the unobserved quantity c j (0) − s in j compromises. We apply the algorithm to the four production scenarios described in the main text. First, in the linear-only scenario (LIN) we assume all firms have only non-essential inputs, i.e. I ne i = {01, . . . , 99}. Second, in the Leontief-only scenario (LEO) we assume all firms only have essential inputs, i.e. I es i = {01, . . . , 99}. LIN (LEO) is the lower (upper) bound for shock propagation in our model. Third, in the mixed scenario (MIX) we assume all firms within NACE classes 01-45 (physical production) have only essential inputs I es i = {01, . . . , 99} and all firms within NACE classes 46-99 (non-physical production) have only non-essential inputs I ne i = {01, . . . , 99}. Fourth, in our baseline scenario (GL), we assume that for firms within NACE 01-45 the set of essential inputs consists of I es i = {01, . . . , 45} only physical producers and the set of non-essential inputs is I ne i = {46, . . . , 99} all inputs from non physical producers, while for firms with NACE 46-99 we assume they have only non-essential inputs I ne i = {01, . . . , 99}. We assume that firms within NACE 01-45 have a physical production process, while firms within NACE 46-99 trade or provide services.

To calibrate the production functions and to model the spreading of shocks, we use the supply network, W , consisting of all relevant supplier buyer transactions of firms within Hungary. For the current study we obtained data for the years 2017 and 2016. The firm level supplier-buyer connection data is collected by the National Tax and Customs Administration of Hungary as a part of the VAT reporting of firms. We got access to this dataset in an anonymized format through the Central Bank of Hungary. The data contains trade links among Hungarian companies between 2014 and 2017, when the tax content of the transactions between two firms exceeds HUF 1 million (EUR 3,000) in the given year. The cleaning and filtering of the data was based on (Borsos and Stancsics, 2020) . The two most important corrections are described below. Many of the links in the network disappear between the observed periods and new links emerge to a similar extent. This happens mainly because of the presence of many one-off, incidental transactions, which relationships are not particularly relevant from the point of view of supply chain contagion. As these links increase the noise in the data, we filtered the network to contain only long-term supplier connections. We consider a link long-term if there were at least two trade events between the parties and if there is at least 90 days time difference between the first and the last transaction. Under these mild requirements, only 54% of the links are long-term, however, these cover 93% of the aggregate trade volume in the network (in 2017).

In the Hungarian VAT regulation there is no general rule that states if firms belonging to the same ownership based group should report individually or at the group level. To handle the potential distortions arising from this inconsistency, we collapsed the supplier network to the group level in every case based on the ownership data obtained from the OPTEN database. The details of this procedure are described in (Borsos and Stancsics, 2020) .

To obtain not only the transactions between nodes, but also information about the nodes the VAT data set was enriched with anonymized firm level information from the Hungarian corporate tax dataset. This contains the NACE classification on the 4 digit level 4 for 64,053 firms (in 2017) as well as revenue and material cost information. For those firms where no NACE classification is available, we merge it into an artificial 569 th category. In the future, classification methods could be used to predict NACE membership based on the input and output links of a company. Symbol-size corresponds to the strength si, red symbols belong to the plateau, emphasised by the shaded area. We find large and small companies in the plateau, suggesting that firm-size is not a good predictor for very high ESRI. Even though we find that for the bulk of the companies strength explains parts of the variations in ESRI (R 2 = 0.90 and slope βreg = 1.11 for regressing log-ESRI on log-strength, p = 2 · 10 −16 ), for individual companies strength does not serve as a reliable predictor of ESRI. (c) Network of the 32 most systemically risky firms (plateau) for 2016. Node size is proportional to the square root of strength, si. Link (from i to j) colors correspond to the downstream "criticality" i.e. the percentage of j's production should i stop producing, Λ d ij . Red thick (blue, thin) links indicate very large (small) losses of production. Small companies predominately supply to large high-risk companies, thereby inheriting systemic risk. Between the companies in the plateau, almost all supply relations are highly critical (red thick). In this sub-network the default of one firm's production will lead to the default of many others.

227,355 links. The core-periphery structure is still visible with some denser regions (darker areas in the center and around); the network becomes substantially sparser towards the periphery. In the center large nodes, highlighed in yellow, are present. Size is proportional to the square root of node strength. This hints to the fact that denser regions are agglomerations around large firms and clusters of large firms. This should be verified with community detection methods in future work. In FIG. 5 (c) we visualize the graph of all suppliers of the largest 3 nodes (5,945 nodes) and all 26,783 links between them. It is obvious that many suppliers of the largest firms are densely connected. To better understand which companies are forming the plateau region of extremely risky companies in Appendix FIG. 6 (b) we show ESRI i as a function of the strength, s i (firm size within the network). Red color indicates the plateau companies; Symbol size represents the strength, s i . Clearly, as for the year 2017, also in 2016 the ESRI of plateau firms (located in the shaded area) is not changing with size; in the plateau we find large and small companies (note the range of strength of 4 orders of magnitude), suggesting that firm-size is not able to predict extreme ESRI values at all also for the year 2016. For the bulk of the companies (blue), we find an overall strong correlation of log-ESRI and log-strength (R 2 = 0.90 and slope β reg = 1.11 for regressing log-ESRI on log-strength, p = 2 · 10 −16 ). For individual companies, strength does not serve as a reliable predictor of the ESRI, since the spread of the ESRI extends over 4 orders of magnitude. Note that all firms in the plateau are physical producers (NACE 01-45).

Figure 6 (c) shows the network between the plateau nodes. The pattern is similar to 2017. There is a strongly connected component of core nodes that are connected with each other via highly critical supply relationships. In the periphery there are mostly small nodes having only one or two critical supply relations to nodes in the strongly connected component.

We investigate the empirical distribution of the ESRI for the years 2017 and 2016 in more detail. In FIG. 7 (a) we plot the empirical cumulative distribution function -indicating the probability P(ESRI i > x) of observing an ESRI value of larger than x -in log-log scale, for the four scenarios, GL, MIX, LIN, and LEO. The three scenarios including shares of Leontief production functions behave differently than the LIN scenario. These three scenarios decay over a wide range of values as a power law (linearly in log-log), followed by a sharp cut-off (the plateau). The linear-only (LIN) scenario does not show this behaviour. We estimate the power-law exponents for the three scenarios and indicate the corresponding slopes with a red dashed line for LEO, a grey long-dashed line for MIX and a blue solid line for the GL. For the estimation we use the maximum likelihood estimator (MLE)

see (Hanel et al., 2017) or Eq. (3.1) in (Clauset et al., 2009) . Since, the linear decay does not extent over the whole range of the ESRI values we restrict the estimation intervals. For LEO the estimate gives a slopeα LEO = 0.5 for the values ESRI i ∈ [1.5·10 −6 , 10 −1 ] covering 53% of the observations. For the MIX the estimate gives a slopeα MIX = 0.63 for the values ESRI i ∈ [5 · 10 −6 , 3 · 10 −2 ], covering 29% of the observations. For the GL the estimate gives a slope ofα GL = 0.67 for the values ESRI i ∈ [7 · 10 −6 , 3 · 10 −2 ], covering 24% of the observations. The MLE in Eq. (L1) corresponds to the probability density function, while the plotted empirical ESRI distributions and the estimated slopesα GL ,α MIX ,α LEO in the figure corresponds to the cumulative distribution function that has an exponent α − 1. Figure 7 (b) shows the results for the year 2016. The patterns are practically identical. For GL the estimate of the slope changes slightly from 0.67 to 0.68. Table I contains the number of observations above several threshold values of ESRI in 2017. The value 0.41 is the threshold for the plateau in the LEO scenario (66 firms) and 0.22 is the threshold for the plateau in GL (32 firms), and MIX (47 firms). In all scenarios only a few firms have a high ESRI of more than one percent. Same for 2016 in Table II . We conduct a more extensive regression analysis to investigate if there are firm level factors that explain ESRI i better than total strength, s i . We estimate 5 additional regression models with the dependent variable in log-scale. Details of the model specification and results are shown in Table III . Each column corresponds to one model fit. First, we regress systemic risk on log-in-strength, log s in i , (first column) and on log-out-strength, log s out i , (second column) separately and jointly (third column). Then, we regress log-systemic risk on log-total-strength (fourth column), log s i , and on log-revenue (fith column). Additionally, we control for the number of customers, the number of suppliers in the network, market shares of firms in their respective NACE 4 sector and add fixed effects for the NACE 4 industry affiliation of firms. Note that market share is important because it determines the degree by which firms can be replaced as a supplier. Note that for in-strength and out-strength not all observations can be included, because there are firms without in-or out-links in the production network dataset.

We see that all independent variables are highly significant, which is not surprising given the large number of observations. Note that log-out-strength seems to be a better explanatory variable (R 2 = 0.82) than log-in-strength (R 2 = 0.73) for log-ESRI. The likely reason is that due to the non-linearity of the downstream contagion (for firms with essential inputs) downstream effects weigh larger than upstream effects. Interestingly, explanatory power (R 2 = 0.85) increases only slightly when adding both log-in-and log-out-strength. The difference in sample size makes a direct comparisons difficult. Interestingly, the explanatory power of log-total-strength (R 2 = 0.91) seems to be a better explanatory variable than log-in-strength and log-out-strength jointly. This can be attributed to the difference in the observations used in the model estimation. Note that the model containing log-total-strength consists of all 91,515 observations while the model containing log-in-and log-out-strength contains 35,058 observations. The ESRI of firms having both in-and out-links seems to be more difficult to explain than the one for firms having only in-or out-links. Log-revenue can explain less than either of the other four models. This is not surprising, because revenue captures also sales transactions that are not present in the Hungarian production network and this leads to a distorted picture of firm size within the network. For example, a firm with a large revenue can export almost its entire production and thus cause relatively small downstream contagion, while the firm can import most of its inputs and thus causes only little upstream contagion. The additional control variable market share and the industry fixed effects increased the R 2 only marginally.

There is a clear indication that size proxies like strength seem to explain large parts in the variation of ESRI. However, as seen in main text FIG. 2 (b) on the individual firm level for a given level of strength there exists still an extremely large variation in ESRI.

We investigate the changes of the systemic risk index, ESRI i , and firm strength, s i from 2016 to 2017. For a comprehensive view we present respective bi-variate density plots in log-log scale in FIG. 8 (a) and (b) . It is clear that . 8 (b) . Again most values are found on the diagonal indicated by red and dark blue colors. Outliers are visible (yellow). Strength of smaller firms shows much larger fluctuations than the strength of of larger firms, seen by a diminishing variance around the diagonal for large values. The correlation of firm strength of 2016 and 2017 in linear scale is ρ ESRI 16 , ESRI 17 = 0.86 (p-value 2 10 −16 ) and correlation in log scale is ρ log(ESRI 16 ), log(ESRI 17 ) = 0.88 (p-value 2 10 −16 ). This confirms the intuition obtained from the visual inspection. There are three main observations. First, the bulk of values is clustered at the center with small changes in both variables. Second, there is a group of firms exhibiting a positive relation between changes in strength and ESRI. Third, there is a group of firms exhibiting large variation in ESRI, but almost no variation in strength. Overall this indicates again that changes in strength seems to have some influence on changes in ESRI, but for many firms the changes in strength are not predictive for changes in ESRI. It seems that it matters more with whom firms form buyer-supplier relations than how large the sum of these relations is. For completeness Here we explain deeper why sector level aggregation of the firm level production network leads to distortions in the picture of shock propagation. If a sector contains at least two firms, A and B, with linearly independent input sector vectors or customer sector vectors, it is possible to construct three firm level shock scenarios that have the same size on the sector level, but are in fact heterogeneous shocks on the firm level (single companies are affect differently). Consequently, these three shock scenarios affect the direct input-and customer-sectors in different ways and hence actual shocks will propagate very differently for each of the scenarios. Note that shocks only propagate differently when analyzed on the firm level, but when the production network is aggregated to the sector level the cascades triggered by these three initial shocks are the same.

We present a concrete example where three shock scenarios are indistinguishable in terms of sector size affected, but have very different effects on how other sectors are affected by them. The implied assumption of a classical sector shock of size x% is that all companies i within the sector are affected by the same percentage shock of x i % = x%. In our example the sector shock is the baseline scenario. We simulate a 18% shock to each company in NACE sector 2611 (Manufacturing of electronic components). The percentage shock is measured in percent of the sector's strength, s 2611 = j : pj =2611 s j , that is initially affected. Then, we construct two firm level shocks that affect also 18% of the strength s 2611 of sector 2611, but are not uniformly distributed across the sector's firms. In the first scenario we apply a shock to a single company -called A for anonymity reasons-of 100% (de facto a temporary failure). Note that s A /s 2611 = 0.18. In the same way we apply a shock of 59% to firm B (we assume that 59% of the firms inputs and outputs are not bought and supplied). Note that 0.59s A /s 2611 = 0.18. Thus all three shocks are the same on the sector level in the classical sense, but affect different companies within the sector to a different degree. We know from FIG. 9 (b) and main text FIG. 3 (b) that firm A and B have different customer-and input-sector vectors, consequently, first order shocks spread to different sectors. This leads to a distinct spreading of shocks within the whole network and we expect other sectors to receive different indirect shocks for the three scenarios.

In FIG. 9 (a) we show how all the other 567 NACE 4 sectors are affected, measured in percent of output lost, i.e., for each NACE sector (on the x-axis) the y-axis shows the percentage of the shock this NACE sector received in response to the initial shocks from the three scenarios uniform (black line) shock to firm A (red) and shock to firm B (blue). We sorted the NACE 4 sectors on the x-axis in ascending order with respect to the baseline scenario (the FIG. 9 (a) Effects of different shock scenarios to NACE sector 2611 (Manufacture of electronic components). Shocks are of the same size but affect different companies in the other 567 NACE 4 sectors. The black line shows the baseline shock scenario where a 18% shock is applied to all companies in sector 2611. Sectors are sorted in ascending order w.r.t. the baseline scenario. The red curve corresponds to a 100% shock to firm A (corresponds to an 18% shock on sector). The blue curve represents a 59% shock to firm B (corresponds to an 18% shock on sector). so that the three shock scenarios are all of the same size on the sector level. The choice of the initially shocked companies have a drastically different effect on how all the other sectors are affected. (b) Customer sector vectors (at NACE 4 level) for 69 firms in sector 2611. Even though all companies belong to the same sector, their customer sector vectors are substantially different. We highlight the customer sectors of firm A (red, 5 distinct customer sectors) and firm B (blue highlight, 4 customer sectors). They differ substantially and have only 1 common customer sector. uniform sector shock), i.e., to the very left there is the NACE category that is affected least by the 18% shock to all firms in sector 2611. We see that the effects of the shock to firm A are in general larger than the shock to firm B, but this is not the case for every single sector. The blue line is below the baseline scenario in most cases (shock to firm A affects sectors less), while in several cases the red line is above the baseline scenario (shock to firm B affects these sectors stronger than the baseline). Note that the red (blue) line is above the black line only in 37% (6%) of the cases. This means that looking at the sector level would overestimate the shocks received by the NACE 4 sectors in 63% (94 %) of the cases and underestimate it in the other cases if the true shock was a shock to firm B (firm A). The correlation between the percentage of output lost per sector, caused by the initial shock scenarios, is given in Table  IV . We see that the two firm level shocks have a low, but significant correlation of 0.11 (p-value p = 0.0085). Both firm level scenarios are positively correlated with the baseline shock scenario. In FIG. 9 (b) we show for all companies in NACE sector 2611 the customer sectors at the NACE 4-digit level (one row corresponds to one firm, a column to a NACE 4 class). Companies (rows) are sorted with respect to the similarity of their customer sectors. Even though all companies belong to the same fine grained industry sector, their customer sectors vary substantially. The customer sectors of firm A (red shading) and firm B (blue shading) have a small overlap seen in a Jaccard index of 0.13. Depending on which of two companies fail, very different customer sectors are affected. Among the 51 firms (in NACE 2611 that have NACE 4 classified inputs) 85% have no pairwise customer sector overlap (Jaccard index of zero). This means that when choosing two firms randomly within the sector the probability of having no common customer sector is 85%. Consequently, for theses cases not only different companies, but entirely different economic sectors will be affected, depending on which of the two randomly selected firms suffers the initial shock. The most common customer sector is supplied to by 30% of the companies in 2611 and the most common input sector is a supplier to 35% of the companies in 2611. Note that 19 firms have no input sectors (empty rows).

Executive order on america's supply chains

Meat-Shortage Risk Climbs with 25% of U.S. Pork Capacity Offline

Unfolding the hidden structure of the Hungarian multi-layer firm network

International Conference on Computational Science

The Belgian production network

Archiv fur Sozialwissenschaft und Sozialpolitik 60

Understanding National Accounts

Heterogeneous firms and the micro origins of aggregate fluctuations

The Review of Economic Studies

Input-Output Analysis: Foundations and Extensions

Consumption Tax Trends 2020

Modeling simultaneous supply and demand shocks in input-output networks

Production Networks and Epidemic Spreading: How to Restart the UK Economy?

The Economist, 2021, How vaccines are made, and why it is hard

forthcoming, A Complex Systems Perspective on Macroprudential Regulation

Intermediate Microeconomics: A Modern Approach: Ninth International Student Edition

Essays on the interface between finance and technology

the International Society for Inventory Research

This work was supported in part by the OeNB Hochschuljubiläumsfund P17795 the Austrian Science Promotion Agency FFG under 857136, the Austrian Science Fund FWF under P29252, H2020 SoBigData-PlusPlus grant agreement ID 871042.

To give a more detailed view on the network structure, we show three additional aspects of it in FIG. 5. (a) shows a section of the network with 5,541 nodes and all 27,793 links between them. It is based on the same randomly sampled 1,500 nodes used for the visualization in FIG. 1 (c) in the main text. In FIG. 1 (c) we sample 1,500 random nodes and consider all their in-links, yielding 6,113 nodes, and visualizing the giant component consisting of 4070 nodes. Here we consider all links between the 6,113 nodes -not just the in-links of the 1,500 nodes-and visualize the resulting giant component consisting of 5,541 nodes. The giant component becomes substantially larger, and needless to say, denser. This makes the structure more difficult to visualize. The center of the network is denser than the periphery. However, by considering all links, some nodes on the periphery have higher degrees (nodes with high degree on the outer part of the circle). This is most likely due to the nature of considering a random sample. In the overall network they are most likely not in the periphery, but more central. In FIG. 5 (b) we show the giant component of 86,470 nodes and