\section*{Results}

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 \ref{SI:data} and \cite{borsos2020unfolding}. Figure \ref{fig:transaction} (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 {\em chains}, but on a tightly interwoven {\em 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 \ref{SI:network_vis}. 

We next calculate the ${\rm 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. \ref{fig:sr_profile} (a) in linear-log scale for 2017. For 2016, see  FIG. \ref{fig:sr_profile2016} (a) in Appendix \ref{SI:2016results}. 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 \ref{SI:table_large_obs2017} in Appendix \ref{SI:powerlaw}. The situation is similar for the MIX scenario, where 47 companies belong to the plateau with values around $0.23$. For the (unrealistic) reference case, LEO with {\em all} companies of Leontief type (red), for 66 firms we find much higher systemic risk levels of about $0.42$.

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 $\alpha^\text{GL}=0.67$ for GL,  $\alpha^\text{MIX}=0.63$ for MIX, and $\alpha^\text{LEO}=0.50$. For details of the fits, see Appendix \ref{SI:powerlaw},  FIG.  \ref{SI_fig:powerlaw} (a). The shape of the rank-ordered ESRI distribution (plateau and power-law tail) is similar to what was found in Fig 4. in \cite{moran2019may}. As expected, the reference case LIN (green) where {\em 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.  \ref{SI_fig:powerlaw} (b) in Appendix \ref{SI:powerlaw}.

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. \ref{fig:sr_profile} (b) we show ${\rm 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 $\beta_{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 \ref{SI:regression_extended}, confirms that generally firm level quantities fail to explain ESRI that is a network-based measure. Note the relation to the Hulten theorem~\cite{Hulten1978Theorem}, 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 \citep{baqaee2019macroeconomic}. 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 (LIN). To illustrate this point, in FIG. \ref{fig:sr_profile} (d), we show the network spanned by the 32 plateau firms. Node size corresponds to the square root of strength,  $\sqrt{s_i}$; link colors correspond to the direct down-stream impact of the supplier on the buyer node (criticality), $\Lambda^\text{d}_{ij}$, (for the definition, see Appendix \ref{SI:propagation}). Red links indicate highly critical supplier-buyer relations; $\Lambda^\text{d}_{ij}=1$ means that 100\% of $j$'s production fails if $i$ fails. Blue links mark non-critical supplier-buyer relations; $\Lambda^\text{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. 

For systemic risk monitoring of companies it is crucial to track the change of their systemic riskiness (ESRI) over time. To estimate the extent of the typical yearly fluctuations, we compare the ESRI values for the years 2016 and 2017 in FIG. \ref{fig:sr_profile} (c). Size corresponds to firms' log-strength. Of the $68,254$ firms that appear in both years, the values of ESRI in 2017 and 2016 are correlated ($\rho ({\rm ESRI}^{16},{\rm ESRI}^{17})=0.78$,  $p=2.2\: 10^{-16}$). The extent of fluctuations is clearly seen in log-log representation in FIG. \ref{fig:sri_16_17_density} (a) in Appendix \ref{SI_annualchange}. For the logged variables we observe an even higher correlation ($\rho ({\rm log(ESRI)}^{16},{\rm log(ESRI)}^{17})=0.85$, $p=2.2\: 10^{-16}$). Figure \ref{fig:sri_16_17_density} (c) shows that relative changes in strength explain the relative changes in ESRI only partially.

In FIG. \ref{fig:sr_profile} (c) temporal changes of companies in the plateau between 2016 and 2017 are seen. The points in the horizontal shaded region are those companies in the plateau in 2017 that had (much) less ESR in 2016. The points in the vertical shaded area are those companies that are in the plateau in 2016, but not 2017. These reduced their ESRI. From the 32 plateau companies in 2016, 20 remain in the plateau (intersection of shaded regions). Of the 12 companies leaving the plateau from 2016 to 2017 some remain very risky (light red) others become less risky (blue). 12 companies enter the plateau in 2017 (also firms not contained in the 2016 dataset). 
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. 

\begin{figure*}[tb]
	\centering
	\includegraphics[width=0.75 \columnwidth]{figures/firm_minus_sector_shock_nace2611.pdf}
	\includegraphics[width=0.75 \columnwidth]{figures/input_overlap_adj2611_n4_sort.pdf}
	\caption{ 
		Importance of firm level analysis. 
		(a) Effect on all 568 NACE 4-digit industry sectors, following a general 18\% shock of the entire NACE sector 2611 (Manufacture of electronic components). 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).
	}
	\label{fig:firm-level}
\end{figure*}

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. \ref{fig:firm-level} (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 {\em initial shock} is exogenously applied and triggers a spreading of shocks in the production network. The {\em 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^\text{2611}=\sum_{i|p_i = 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, \emph{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\cdot 10^{-6}$). In Appendix \ref{SI:firm_level_advantages} 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. \ref{fig:firm-level} (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 \ref{SI:firm_level_advantages}. 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. 


\section*{Discussion}

Based on the reconstruction of all relevant buyer-supplier 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 $\overline{{\rm ESRI}_i} = 0.00018$ (mean). The median is $1.7 \cdot 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 critical 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 \cite{battiston2012debtrank, markose2012systemic, thurner2013debtrank}. 

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 down-stream 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 \cite{pichler2021modeling}. 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 \ref{SI:supplier_replaceability}. 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. 

{\bf 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 \citep{poledna2016elimination}). 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 \cite{Solomon2021} or the USA \citep{biden2021}. The possibilities for systemic risk monitoring as presented in this paper should be kept in mind, when designing new regulation. 

{\bf 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. 

{\bf Avoiding risk concentration.}
In \cite{poledna2016elimination} 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). 

{\bf 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.

{\bf 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 \cite{geodis2017}. 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 \cite{diem2020minimal, pichler2021systemic}. 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 \cite{moran2019may}, 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.


\section*{Material and Methods}


	The economic systemic risk index for firm $i$ is calculated as
	\begin{equation} \label{eq:matmethods_esri}
		{\rm ESRI}_{i} = \sum_{j=1}^{n} \frac{s^\text{out}_j }{\sum_{l=1}^n s^\text{out}_l }\big(1-h_j(T) \big) \quad ,
	\end{equation} 
	where $h_j(T) = \frac{x_j(T)}{x_j(0)}$ is the fraction of the output remaining of firm $j$ at time $T$, where $T$ is the convergence time, see Appendix \ref{SI:ESRI}. ESRI$_i$, is interpreted as the fraction of output that is lost in the entire production network in response to the failure of individual company $i$, given that its supply and demand can not be replaced. The vector $h(T)=\min(h^d(T),h^u(T))$ is the result of propagating the shock from firm $i$'s initial default downstream along the out-links by updating
	\begin{eqnarray}\label{eq:matmethods_downstream}
		x_l^{\text{d}}(t+1) & = &   f_j\Bigg(\sum_{j=1}^{n} W_{ji} h_j^\text{d}(t) \delta_{p_j,1},  \dots, \\ 
		& & \qquad \sum_{j=1}^{n} W_{ji}h_j^\text{d}(t) \delta_{p_j,m}\Bigg) \quad , \notag
	\end{eqnarray}
	and upstream along the in-links by updating 
	\begin{eqnarray}\label{eq:matmethods_upstream}
		x_l^{\text{u}}(t+1) & = &    \sum_{j=1}^{n} W_{lj}h_j^\text{u}(t) \quad. 
	\end{eqnarray} 
	For details of the calibration to the GL, see Appendix \ref{SI:ESRI}, Appendix \ref{SI:derivation}, for the implementation, see Appendix \ref{SI:propagation}.
	






