Commodity-specific triads in the Dutch inter-industry production network

Triadic motifs are the smallest building blocks of higher-order interactions in complex networks and can be detected as over-occurrences with respect to null models with only pair-wise interactions. Recently, the motif structure of production networks has attracted attention in light of its possible role in the propagation of economic shocks. However, its characterization at the level of individual commodities is still poorly understood. Here we analyze both binary and weighted triadic motifs in the Dutch inter-industry production network disaggregated at the level of 187 commodity groups, which Statistics Netherlands reconstructed from National Accounts registers, surveys and known empirical data. We introduce appropriate null models that filter out node heterogeneity and the strong effects of link reciprocity and find that, while the aggregate network that overlays all products is characterized by a multitude of triadic motifs, most single-product layers feature no significant motif, and roughly 85% of the layers feature only two motifs or less. This result paves the way for identifying a simple ‘triadic fingerprint’ of each commodity and for reconstructing most product-specific networks from partial information in a pairwise fashion by controlling for their reciprocity structure. We discuss how these results can help statistical bureaus identify fine-grained information in structural analyses of interest for policymakers.

In the last decade, the increasing availability of data at the industry and firm level led to a vast number of studies analyzing the system of customer-supplier trade relationships -the production network -among industries [1][2][3][4][5][6] or firms  and their impact on country-level macroeconomic statistics [40].
The heterogeneity encoded in the production network structure plays an essential role in amplifying economic growth [34] and in the propagation of shocks [1,12] related to exogenous events, such as Hurricane Sandy [29], the Great East Asian Earthquake [11,27], the Covid-19 pandemic [6,17,28], or endogenous events such as the 2008 financial crisis [41,42].
Even in the time of globalization -characterized by highly interconnected global supply chains -domestic production networks are still relevant.In fact, it has been shown that for a small country as Belgium, while almost all firms directly or indirectly import and export to foreign firms, these exchanges represent the minority of domestic firms' total revenues [16].
While aggregated information about single firms is contained in most National Statistical Institutes' repositories, reliable data on input/output relationships is available only for a small number of countries.For instance, the Compustat dataset contains the major customers of the publicly listed firms in the USA [8].The FactSet Revere dataset contains major customers of publicly listed firms at a global level, with a focus on the USA, Europe, and Asia [30].Two datasets are commercially available in Japan, namely, the dataset collected by Tokyo Shoko Research Ltd. (TSR) [11] and the one collected by Teikoku DataBank Inc. (TDB) [35].They are characterized by a high coverage of Japanese firms but with a limited amount of commercial partners.Other domestic datasets contain transaction values among VAT-liable firms: this is the case for countries such as Brazil [14], Belgium [15], Hungary [17], Ecuador [7], Kenya [19], Turkey [20], Spain [21], Rwanda and Uganda [22], West Bengal [23]; or contain transaction values among the totality of registered domestic firms such as in the case of Dominican Republic [18] and Costa Rica [43].
However, in production networks, user firms connect to supplier firms to buy goods for their own production.Customer-supplier relationships are, hence, characterized by an intrinsic product granularity that is usually neglected.The importance of product-specific information has been highlighted, for instance, in a rare study that utilizes surveys with limited data for Japanese automotive firms [44].Generally, in the economic theory of industries and firms, the problem of product granularity is 'solved' artificially, by assuming that industries/firms supply a single product [1,4].This is an oversimplification that often conflicts with reality: indeed, a single firm can possess more than a production pipeline and is capable of supplying multiple products (e.g.Samsung, a Telecommunication company, sells also household appliances, and multinational companies such as Amazon and Google supply a large number of different products).
Recently, Statistics Netherlands (CBS) reconstructed from national statistics two multi-layer production network datasets for domestic intermediate trade of Dutch firms for 2012 [25] and 2018 [10], with each layer corresponding to a different product exchanged by a firm for its own production process, as illustratively depicted in Fig. 1(a).The 2012 dataset has been recently used to prove the complementarity structure of production networks [33] by inspecting the number of cycles of order 3 and 4 compared to a null model taking into consideration the in-degree and out-degree distributions.We use the improved version for 2018 and construct an inter-industry network that will be presented in the next section.
In this study, we focus on triadic motifs and anti-motifs that are over-occurrences and under-occurrences of different patterns of directed triadic connections, respectively.They are represented in Fig. 1(b).Triadic and tetradic connections are known as the building blocks of complex networks [45], playing the role of functional modules or evolutionary signs in biological networks [46,47], homophily-driven connections in social networks [48], complementarity-driven structures in production networks [33,36], their change in time being interpreted as self-organizing processes in the World Trade Web (WTW) [49,50], and early-warning signals of topological collapse in inter-bank networks [51,52] and stock market networks [42].It has been proven that for the majority of (available) real-world networks, the triadic structure is maximally random [53] and by fixing it their global structure is statistically determined [54].
In contrast, research on weighted motifs and anti-motifs is still underdeveloped.To our knowledge, only one study involves trade volumes circulating on triadic subgraphs, using a probabilistic model based on random walks on the WTW [55].
Motif detection strictly depends not only on the properties of the real network but also on the randomization method used for the computation of random expectations.In Network Science literature, various methods have been advanced for network randomization, primarily edge-stub methods, edge-swapping methods, and Maximum-Entropy methods, we focus on the latter.Randomization methods based on Entropy Maximization [56][57][58] build Graph Probability distributions that are maximally random by construction.Available global or node-specific data are encoded as constraints in the optimization procedure, and their corresponding Lagrange Multipliers are computed by Maximum Likelihood Estimation (MLE) [59].This theoretical framework has been proven to successfully reconstruct economic and financial systems [60][61][62][63], statistically predicting both the topology and the weights of the WTW [64][65][66], in an integrated [67], or conditional fashion [68], with only structural constraints, or informing the models with economic factors [69][70][71], statistically predicting banks' risk exposures [72], and most recently, statistically reconstructing payment flows among Dutch firms that were clients of ABN Amro Bank or ING Bank, constraining their industry-specific production functions [26].These methods have been proven to give the best insurance of unbiasedness with respect to missing data, as seen by independent testing [73][74][75][76].Two studies using Maximum-Entropy modeling are especially worthy of note for motif detection: a theoretical study where the authors develop null models for triadic motif detections and compute z-scores of triadic occurrences analytically [77], and an applied study where triadic motifs and their time evolution are used as early warnings of topological collapse during the 2008 financial crisis [51].
Our contribution goes in this direction, using Maximum-Entropy methods constraining degree distributions and strength distributions -in their directed form and taking into account their reciprocal nature -to characterize triadic connections and the total money circulating on them for different product layers of the reconstructed Dutch production network.An analysis of this kind can give better insight into how much product-level granularity is needed in production network datasets and how the links and weights of a production network are organized for different products.Once product layer patterns have been detected, National Bureau officials -having experience in the domestic trade of that single commodity -can infer if such motifs and anti-motifs are due to commodity-specific characteristics, market imbalances, or represent structures aided by laws.If unbalances and anomalies are detected, executive government agencies can use this as input to eventually advance policy laws to nudge a more convenient redistribution of connections and trade volumes.

The CBS Production Network
The CBS production network for 2018 [10] improves on the 2012 version by integrating more auxiliary micro and industry-level data.Before going into detail it is helpful to explain which are the industry classifications and product classifications used by Statistics Netherlands.Industries are classified using the Dutch Standard Industrial Classifications, in brief SBI, which are equivalent to the European Standard Classification NACE rev.2 in the first two digits, although the subsequent digits can differ.Statistics Netherlands has industry data on two different levels, the SBI4 level, containing 132 industries, and the SBI5 level, corresponding to 888 industries.Regarding the CPA product classification, Statistics Netherlands uses a modified version of the original European CPA, mainly at 4 and 6 digits.In the data, we retrieved 192 commodities for 4 digits and 623 commodities at the 6-digit level.
Firm-level data is obtained from the Statistical Business Register (SBR) for 2018 for over 1 700 000 firms.The SBR contains values about net turnover, geographical location, business id, and business sector at the SBI5 classification level.After cleaning for micro-firms with annual net turnover below 10 000e, around 900 000 firms remain, accounting for 99.5% of the Dutch economy output in 2018.Details regarding the breakdown of output and input at the commodity level are primarily available at the industry level and for a limited number of firms.The Dutch National Supply-Use tables provide data on inter-industry and intra-industry intermediate input/output transactions for various commodities, classified at the Dutch CPA 4-digit level with industries categorized according to the SBI4 level.
While industry-wide transactions are validated, estimating output and input for individual firms per commodity and matching suppliers with users within commodity layers remains a challenge.To estimate supply per firm, domestic turnover, calculated as VAT turnover minus export turnover, is employed as a distributional key.Firms are assumed to supply in proportion to the ratio of their domestic turnover to the overall industry turnover.Additional adjustments are made for wholesale and retail trade firms to account only for domestic turnover associated with actual production.
Estimating use per firm from VAT turnover data involves determining the ratio between intermediate use and turnover.This ratio is estimated using SBS survey data.The breakdown of supply/use per firm at the commodity level is available for a relatively large number of firms through surveys conducted by Structural Business Statistics (SBS) for commercial firms, Prodcom for manufacturing firms, and estimates generated by National Accounts for non-commercial firms.
Specifically, SBS provides a breakdown of sales and intermediate purchases into ten to twenty commodity categories for small firms and at the CPA-classification level for large firms.Prodcom conducts a similar survey.SBS categories are then mapped into CPA commodities by National Account experts.
For firms not covered by the aforementioned surveys, the breakdown in commodities of intermediate supply/use is estimated using the distribution of the industries as a whole from the Supply-Use tables.This approximation can result in implausible values of annual supply and use.To address this issue, thresholds are imposed, setting supply values below 2 000e and also use values below 1 000e to zero.Finally, an Iterative Proportional Fitting (IPF) procedure is implemented to ensure consistency with industry-level Tables.
Once supply and use per firm per commodity are obtained, their out-degree distribution is estimated using stylized facts from Japanese firms [39], connecting out-degrees with firm sizes through a power-law function, while their indegree distribution is estimated assuming a power-law connection to firm-specific input at the commodity-level, an assumption that is consistent with recent studies on Dutch inter-firm payments [26].
Once in-degree and out-degree distributions per commodity are estimated for each firm, suppliers and users are matched according to a deterministic procedure that takes into account (1) a company score, encoding their net turnover, (2) a distance score, that takes into account their mutual distance, (3) the presence of a link between respective industries in the Supply-Use tables, (4) the presence of the observed relationship in the Dun & Bradstreet dataset, i.e. a dataset containing the list of the users of the 500 largest suppliers in the Dutch Economy.After the computation of the related 'link score', users in each commodity layer are ordered according to their purchase volumes.The top user, then, selects the best X suppliers and establish a connection with them, where X represents its commodity-specific in-degree.The procedure continues from the second-highest purchasing volume user to the last until no available links remain and degree distributions are reproduced.Network weights are then distributed across generated links according to a power-law distribution.Finally, the resulting weighted inter-firm network at the 650 commodity level (National CPA level 6) is compared to the Supply-Use tables (National CPA level 4) and consequent adjustments are made to weights and links.Further details can be found in [10].
We aggregate the inter-firm network at the commodity level, passing from 623 commodities (CPA level 6) to 192 commodities (CPA level 4, compatible with Supply-Use tables).Then, we aggregate firms in industries at the SBI5 level, taking their business sector ids from the SBR.For the topic of interest, the self-loops implied by intra-industry trade are not important and can be removed from the dataset without adversely affecting the subsequent analysis.After cleaning for intra-industry trade, we obtain a multi-layer inter-industry production network containing linkages and weights for 862 industries (nodes) and 187 commodity groups (layers).
The firm-level reconstructed dataset is not without limitations.One source of error arises from the breakdown provided by SBS and Prodcom surveys, particularly regarding the documented intermediate purchases and sales.The purchases may include imports, and the sales may also include sales for final consumptive use.Another source of error stems from the distributions and assumptions made for firm out-degree and in-degree distributions.While these assumptions are supported by stylized facts from Japanese firms (for out-degrees) and payment data from a large sample of Dutch firms (for in-degrees), it cannot be assumed that the parameters used in the reconstruction are universally applicable or representative of 'true values'.Finally, the matching procedure results in a deterministic network where the 'best' users have priority in connecting with their more closely aligned suppliers.This algorithm cannot account for noisy behavior or real-world uncertainties.In fact, for the 2012 version, with similar assumptions on degree distributions and the same matching algorithm, it has been demonstrated that these assumptions lead to biases in core network statistics such as the number of links in commodity layers [37], when compared with the ground-truth provided by a known sample of firm-to-firm connections collected by Dun & Bradstreet (for 2012).While aggregation at the SBI5 level is bound to reduce the biases that arose at the firm-level, it is still not clear how much the results are impacted by the propagation of these errors.Further discussion on limitations is provided at the beginning of the Discussion Section.

Network randomization methods
The main goal of Network randomization methods is the generation of a statistical ensemble of networks, which are maximally random given available data.In our case, we randomize each product layer of our industry-multilayer network separately using Maximum-Entropy methods.The available data -encoded as constraints in the Entropy maximization -consists of the supplier's(user's) tendency to supply(use) a specific commodity and its output(input).The obtained statistical ensemble of networks represents the possible realizations of the system taking into account suppliers' and users' tendencies.After the generation of the synthetic ensemble of networks it is possible to extract metrics of interest as ensemble averages.
The null models we take into account are the Directed Binary Configuration Model (DBCM) [65] and the Reciprocal Binary Configuration Model (RBCM) [77] for the estimation of network links, and the Conditional Reconstruction Method A (CReM A ) [68] and the newly developed Conditionally Reciprocal Weighted Configuration Model (CRWCM) for the conditional estimation of network weights.The DBCM corresponds to the model that maximizes the Shannon Entropy attached to the distribution of possible binary adjacency matrices, given that in-degree and out-degree distributions are constrained on average.The RBCM is also used for estimation of links by maximizing the Shannon Entropy attached to the distribution of possible binary adjacency matrices, but makes use of additional information, namely the non-reciprocated out-degree, in-degree and the reciprocated degree distributions.These metrics are originated distinguishing links that are reciprocated from the ones that are not and summing on them.Turning our attention to weighted networks, the CReM A is the Maximum-Entropy model that maximizes the conditional Shannon Entropy attached to the distribution of weighted networks, given the realization of the adjacency matrix A. The constraints used in the conditional optimization are the out-strength and in-strength distributions, corresponding to sum of weights going from and to a node, respectively.The CRWCM, instead, is an augmented version of CReM A , which can take better account of reciprocation by constraining the out-strength and in-strength distributions for  I. Description of the distribution of statistics such as the number of active industries N , the number of links L, the total weight Wtot, the topological reciprocity rt and the weighted reciprocity rw across commodity layers of the inter-industry network.
reciprocated and non-reciprocated links.Both CReM A and CRWCM are estimated using an annealed approach, following the articles [68,80], and consequently coupled with the relative binary model.Specifically directionality is encoded in the DBCM+CReM A model, also denoted as the directed model, while directional and reciprocal information is encoded in the RBCM+CRWCM model, denoted as the reciprocated model.For further information and the mathematical generation of link and weight distributions, please refer to the Method section.

Measuring empirical reciprocity statistics
The presence of data on product granularity gives us the opportunity to study heterogeneity across commodity layers.Let us consider in Table I the number of layer-active industries N , the number of links L, the total weight W tot , and reciprocity measures such as the topological reciprocity r t , defined as the ratio of reciprocated links to L, i.e.
and its weighted counterpart r w , defined as the ratio of total weight on reciprocated links to W , i.e.
The median for N is 149, meaning that for around 50% of commodity layers there are less than 149 active industries (as suppliers or users).At the same time, 25% of commodity layers have less than 62 industries, and another 25% have more than 544 industries.Consequently, industries are specialized among a small number of business activities for half of the commodity groups but, a small, and not negligible, number of layers is characterized by a high number of active industries and hence of industry heterogeneity.Some examples are suppliers of plastic goods that are sold to users with heterogeneous specializations, for instance, Bread, Beer, Cereals, Fish, etc. Also the distributions regarding the number of commodity-specific links L and the related total weight W tot have wide distributions, with a minimum with few digits, respectively 3 and 0.95 (in millions of euro), and a maximum in 5 digits, respectively 15198 and 23767, implying a high degree of heterogeneity in network structure across commodity layers.
Passing from the commodity global statistics to r t and r w , we see a high degree of heterogeneity also in this case, namely a minimum value of 0 stands for layers where no link is reciprocated, i.e. users and suppliers represent two distinct sets of nodes (bipartite graph).Instead, in the majority of the commodities (above 75%) there is a notnull reciprocity.In fact, the median is respectively 0.05 and 0.08.There is also the presence of a small number of commodities (below 10%) which are characterized by a large reciprocity, with a maximum of 0.78 for both r t and r w .
Reciprocity can arise for different reasons: (1) the aggregation from firms to industries or (2) the aggregation of products.To mention the first case, consider two firms A and B in the industry i and other two firms C and D in industry j. Suppose firm A supplies to firm D, while firm C supplies to firm B, in the same commodity layer.Once the firms are aggregated in the related industries, a reciprocated link emerges between them, even if reciprocity is not present at the firm level.The second case follows from the fact that if each commodity layer represents a unique product, that could be represented by the finest CPA product classification (with around 5000 products), and we take into account only intermediate supply and use, it is not reasonable to think that firms are at the same time suppliers and users (of that specific product).Instead, in case of product aggregation, firms may be suppliers of a product inside that commodity group and also users of another product inside that same commodity group.Let us now move to the analysis of triads.We define triadic occurrences N m , the number of times a specific m-subgraph appears and triadic fluxes F m , the total amount of money circulating on each m-subgraph.In Fig. 2, we depict their values normalizing by their sum across the m-types.The normalized N m and F m can be considered as the relative importance of a specific type of triadic subgraph in the network.The aggregated network (depicted in blue), where the weights of all commodity groups are summed, and three commodity layers, namely 'Cereals' (in green), 'Gas/Hot Water/City Heating (in orange) and 'Agricultural Services' (in pink) are displayed.
In the aggregated network, the structures that occur relatively more are m = 1, represented by a supplier connected to two users and m = 13, the totally reciprocated cyclical triad.While m = 13 is probably due to product aggregation, the predominance of m = 1 is a signal of structural dependency on a limited number of suppliers.However, when normalized F m are investigated, m = 13 still contain the majority of the volumes.A similar profile, in the binary case, is given by the Agricultural Services, with the predominance of m = 1 and m = 13.At the same time a relatively smaller amount of money is concentrated on m = 13 with respect to the aggregated case, while m = 1 and m = 11 carry a greater amount of money.During the product disaggregation weights on m = 13 in the aggregated network are redistributed on other subgraphs, especially m = 1.In 'Cereals' and 'Gas/Hot Water/City Heating' these differences are even larger, with a relevant increase of triadic occurrences and fluxes on m = 1, further increasing the dependency of the network on a limited amount of suppliers.Note that when counting the different triads in Fig. 2 they are not nested, i.e. a subgraph of type m = 8 requires two reciprocated links and hence does not contain two subgraphs of type m = 1, which contain only non-reciprocated links.Consequently, the number and fluxes over all triadic subgraphs are structurally independent across different types.

Binary Motif Analysis
We analyze the number of occurrences N m of all the possible triadic connected subgraphs, depicted in Fig. 1(b).To quantify their deviations to randomized expectations, we define the binary z-score of subgraph m where N m (A * ) is the number of occurrences of the m-type subgraph in the empirical adjacency matrix, ⟨N m ⟩ is its model-induced expected number of occurrences, and σ [N m ] is the model-induced standard deviation.An analytical procedure [77] has been developed to compute the binary z-scores for the binary case.However, the assumption on the confidence intervals -represented as the interval (−3, 3) -holds true only if the ensemble distribution of N m is Normal for each m.For all the commodities, m-types, and binary null models, we test the assumption using a Shapiro Test [79].According to the test, N m ensemble distributions are in a large proportion not normal at the 5% confidence level.Consequently, we must use a numeric approach.Networks are sampled according to the DBCM recipe by (1) computing the induced connection probability p ij;DBCM and (2) establishing a link between industry i and j if and only if a uniformly distributed random number u ij ∈ U (0, 1) is below p ij;DBCM .The analogous recipe for RBCM requires (1) computing the set of connection probabilities for non-reciprocated connection between i and j, namely p → ij , p ← ij and p ̸ ↔ ij , and reciprocated connection p ↔ ij , generate a uniform random variable u ij ∈ (0, 1) and ( 2) establishing the appropriate links in the dyad in the following way: • no links from i to j and from j to i otherwise.
In both cases, we generate a realization of A and extract the N m statistic.⟨N m ⟩ and σ [N m ], are the average and standard deviation of N m extracted from the ensemble distribution of 500 realizations of A. After having computed z [N m ], we also extract the 2.5-th and 97.5-th percentiles from the ensemble distribution of N m over all models and we standardize them using Eq. ( 3) by replacing the empirical N m with the percentile.Such measures will serve as the 95% CI for the z-score.
The results for the aggregated inter-industry network are in Fig. 3(a).The z-scores computed with respect to the DBCM are depicted in blue on the left panel, while the z-scores computed with respect to the RBCM are depicted in orange on the right panel.The corresponding confidence intervals at the 5% percent are depicted with the same color (blue or orange) but in slight transparency.The majority of N m are not reproduced by the randomized methods, i.e. the z-scores are outside the confidence intervals.Specifically, only N 8 is reproduced by the DBCM, while both N 1 and N 9 are reproduced by the RBCM.Discounting reciprocal information does not only increase the number of triads that are statistically well described, but potentially changes their type, implying a qualitatively different z-score profile.At the same time, in the aggregated picture, m = 1 and m = 9 are seen as described by a null model implementing reciprocity, i.e. neither high dependency on suppliers (m = 1), nor unstable feedback loops (m = 9), where industries supply to each other in a cyclical fashion, are revealed.The aggregated network, is hence, characterized by a multitude of structures that are not well described by the null model and are due to additional three-node correlations but is relatively resilient to supply shocks and cyclical input/output.By disaggregating from the aggregated monolayer to the multi-commodity network, the majority of commodity-layers have triadic structures which are statistically reproduced by the reciprocal null model.Only 1 or 2 motifs or anti-motifs are present for the majority of the remaining commodities, a result indicating that beneath the aggregated picture, commodity groups are characterized by a small number of commodity-specific motifs and anti-motifs.
In Fig. 3(b-d) z-score profiles for three commodity layers are displayed, namely Cereals, Electrical Components, and the Construction of Tunnels, Waterways, and Roads.RBCM well describes all subgraph occurrences (z Nm is within CI), while the DBCM signals the presence of anti-motifs for m = 10, m = 11 and m = 12 for Cereals, and anti-motif m = 12 and motif m = 13 for the Construction layer.In Fig. 3(e-f) two z-score profiles are displayed -namely for Bread & other Bakery Products and Gasoline -for which RBCM signals the presence of at least a motif or anti-motif.A motif m = 12 is present for the former layer while an anti-motif for m = 4 is present for the latter.Notice that for Bread the DBCM does not signal any motif or anti-motif, implying that deviations can emerge by introducing information on the reciprocal structure.Moreover, subgraph m = 9 in Bread and the majority of subgraphs in the Gasoline commodity layer are characterized by a degenerate Confidence Interval: in all of the generated synthetic networks N m=9 correspond to the empirical N * 9 with null variance, i.e. the constraints imposed on the ensemble totally describe the specific m-type motif, a matter which can arise regardless of the lack of statistics in the related N m .Finally, in Fig. 3(g) the z-profile for the commodity layer Beer/Malt is considered.The DBCM signals a large number of motifs, specifically for m = 2, m = 10, and m = 11, and anti-motifs for m = 3 and m = 8.In contrast, the RBCM signals a lone motif m = 3 and an anti-motif m = 6.In Fig. 4(a), the empirical counter cumulative distribution for the number of deviating binary triads is shown.Introducing reciprocal structure information reduces the number of motifs and anti-motifs present across commodities.For instance, the percentage of commodities with at least a motif or anti-motif is 61% when compared to the DBCM, and 48% when compared to the RBCM, while the percentage of commodities having at least two motifs or anti-motifs is 46% when compared to the DBCM and 27% when compared to the RBCM.
Lastly, we identify the occurrence of m-type of motifs and anti-motifs across commodities by introducing two quantities, c h (m) and c l (m).c h (m) represents the number of commodities having a motif of type m while c l (m) represents the same measure for anti-motifs.The addition of the reciprocal structure reduces the number of commodity-specific motifs for each subgraph type, with the exception of motif m = 6 as depicted in Fig. 4(b), and the number of antimotifs for each type, with the exception of anti-motif m = 8 as depicted in Fig. 4(c).The reciprocal null model, hence, reveals a higher number of commodities that are relatively more vulnerable to demand shock due to bankruptcy of industries of type k in triadic formations m = 6, while it reveals an increased resilience to supply/demand shocks originating from bankruptcy of industries of type j in formations m = 8.

Weighted Motif Analysis
While the bankruptcy of an entire industry is unrealistic, a shock due to a decrease in the flow of goods among industries can propagate along the supply chain, with side effects on the real economy.This implies that not only binary information is important for shock propagation but also weighted information, namely the amount of money circulating on connected structures.
Consider the triadic flux F m on motif m, defined as the total money circulating on triadic subgraphs of type m.We characterize the deviation of empirical F m to null models by defining the weighted z-scores as where ⟨F m ⟩ is the model-induced average amount of money circulating on motif m and σ [F m ] represents the modelinduced standard deviation over the ensemble of network realizations.The theoretical benchmark (or null model) is built by using a combination of binary and conditional weighted models, depending on the wanted constraints.If we deem reciprocal information of negligible importance we should use the combination of models given by DBCM, for the sampling of the binary adjacency matrix, and the CReM A , constraining the out-strength and in-strength sequences.If we deem reciprocal information necessary, a combination of the RBCM and CRWCM should be used.We compare here the two to establish the importance of the addition of reciprocity information for the detection of weighted motifs.
In operative terms, using a two-step model such as the DBCM+CReM A reduces to (1) establishing a link between industries i and j when a uniform random number u ij ∈ U (0, 1) is such that u ij ≤ p ij;DBCM , (2) if i and j are connected, sampling w ij by using the inverse transform sampling method technique, i.e., we generate a uniformly distributed random variable η ij ∈ U (0, 1) such that then we invert the relationship finding the weight v ij to load on the link (i, j).
The network sampling for the RBCM+CRWCM follows the same concepts with two major differences: (1) a link is established using the RBCM recipe and (2) the dyadic conditional weight probability q CReM A (w ij |a ij = 1) is substituted with q CRW CM (w ij |a ij = 1) in the inverse transform sampling.
In Fig. 5(a) the z-score profile for the aggregated network with a single representative commodity is depicted using the directed (in blue on the left panel) or the reciprocal models (in orange on the right panel).There is a large number of motifs and anti-motifs when the benchmark model is directed, only F 3 does not deviate significantly.
When reciprocity information is considered, the picture changes: only three motifs, namely m = 6, m = 11, and m = 13, are identified, and four anti-motifs, namely m = 3, m = 8, m = 10, and m = 12, are found when the reciprocal null model is employed.This model's enhanced accuracy unveils a higher-than-expected volume of financial activity on sub-types characterized by a single exclusive user and two suppliers utilizing each other's products (m = 6), two users supplying to each other while employing a product from the same supplier (m = 11), and entirely cyclical triads (m = 13).In contrast, a lower-than-expected level of financial activity transpires in open triads with two reciprocated ties (m = 8), one reciprocated link and one exclusive user (m = 3), or in closed triads of type m = 10 and m = 12.While it might be contended that the heightened concentration of funds on m = 13 is attributable to aggregation bias, it is crucial to recognize that aggregation solely accounts for the increased monetary worth of the particular sub-type in absolute terms, not for the weighted motif obtained after adjusting for the statistical null model.It should be noticed that the emergence of these specific motifs cannot be easily explained without delving into greater detail, given the representative commodity scheme, while the picture cannot be merely reduced to a higher activity on open triads and a lower activity on closed triads.
Similarly to the binary case, passing from the aggregated network to the disaggregated product-level layers, it is possible to identify a small number of commodity-specific weighted motifs and anti-motifs.
In Fig. 5(b-d) three commodity layers are depicted for which no motifs and anti-motifs are present when z-scores are computed using the reciprocal model.They are 'Seeds', 'Metal components for Doors & Windows' and 'Airline Services'.In the 'Seeds' layer, the directed model signals the presence of an anti-motif for m = 5.In the second layer, no deviations are registered by both null models but CIs are of different nature, in fact, the reciprocal model allows a more restricted range of z-scores with respect to the directed model for m = 9.In the 'Airline Services' layer, for both models, no deviations are present and three CIs are degenerate for m = 5, m = 9, and m = 10.In Fig. 5(e-f) the z-scores relative to the commodity groups 'Coffee/Tea' and 'Textile raw materials and products' are depicted, for which 1 motif is present by using the reciprocal model.For both the directed and reciprocal models there is a weighted motif m = 2 in the 'Coffee/Tea' layer.In contrast, in the Textile products layer the directed model signals an anti-motif for m = 2, while the reciprocal model signals a motif for m = 1.If Fig. 5(g) the z-score profile for the commodity layer 'Shipping Services' is shown: the directed model signals a large number of anti-motifs, specifically for m = 5, m = 7 and m = 12, while it registers a motif for m = 11.The reciprocal model, instead, registers a motif for m = 4 and anti-motifs for m = 5 and m = 12.Different commodity layers call for different motifs and anti-motifs which are due to their specific characteristics.In this paper, we refrain from characterizing every single commodity layer, but a specific and thorough analysis is possible by visualizing the number of triadic sub-types, the z-score profile for N m and their weighted analogs.
The empirical counter cumulative distribution ECCDF(# deviating W∆) for the number of deviating weighted triads is depicted in Fig. 6(a).The number of deviating triadic fluxes is steadily lower using the reciprocal model.F m are maximally random for 49% when the directed model benchmark is used and for 55% according to the reciprocal model.The reduction of the number of motifs is however not as significant as in the binary case.
In Fig. 6(b-c) we plot the weighted analogous of c h (m) and c l (m).Reciprocal information decreases the occurrence of all types of anti-motifs across commodities, with the exception of m = 8.Instead, the profile induced by c h (m) is significantly different using the two null models.For instance, according to the directed model, F 1 is almost always well predicted, instead, it is the most occurring motif according to the reciprocal model.At the same time, reciprocity unveils the dependency of more than 40 commodity layers on the supply of a limited amount of suppliers, which in this case control the market.In fact, the high presence of m = 1 weighted motif signals the vulnerability of the industry-industry network to supply shocks provoked by a reduction of supply volumes.

DISCUSSION
The study of triadic motifs on production networks is still in its infancy due to a scarcity of reliable data.In the existing literature, only binary triadic motifs on one production network, the Japanese one, have been characterized for a single representative commodity [36], while the Hungarian dataset has been analyzed only for triadic occurrences without recurring to a null model [81].The Japanese study revealed a simple but significant pattern: open triadic subgraphs are over-represented while closed triadic subgraphs are under-represented.This phenomenon was attributed to complementarity, where economic actors connect in tetradic structures -better explained by open triads -due to complementary needs [33].
Our findings corroborate the notion that an analysis based on a single representative commodity is insufficient to fully characterize a production network.Product-level data is essential for disaggregating the network into layers that are characterized by commodity-specific binary motifs and anti-motifs.Moreover, we found that the majority of layers exhibit maximally random triadic structures when the reciprocal structure is considered.
At the level of binary motifs, we detected that cyclical reciprocated triadic subgraphs, which are dominant in the aggregated network, break up in the disaggregated product layers, where open triangles become dominant, especially m = 1.However, using the RBCM as a benchmark, we proved that m = 1 is always well described.Conversely, the completely cyclical triads, even if partially broken in the disaggregated layers, are often over-represented compared to the benchmark estimate.In general, constraining the reciprocation capacity of industries -by constraining the reciprocated degrees -is of the foremost importance when characterizing triadic motifs, as explained by the better accuracy and the decrease in binary triadic motifs and anti-motifs when using RBCM as a benchmark compared to DBCM.
We also characterized weighted motifs and anti-motifs, defined as the amount of money circulating on triadic subgraphs, with a novel model which constrains strengths, decomposing them according to the character of the corresponding links.This type of analysis is totally novel in the context of production networks, and rarely seen with benchmark models [55].We find a non-trivial result already when analyzing the aggregated network, subgraphs that are well explained in binary terms -their occurrence is well described by the statistical ensemble induced by the DBCM or RBCM -can be not well described in weighted terms, meaning that even if a binary triadic subgraph has the expected occurrence it can accommodate an unexpected concentration of money.Furthermore, we identified a high presence of m = 1 weighted motifs across commodity layers, a signal of commodity-specific dependency on a limited number of suppliers, which control the market.This implies that a large number of layers are vulnerable to supply shocks, which can arise due to a decrease in supplied volumes (and not only to the supplier's bankruptcy as in the binary case).
Changing the benchmark from a directed to a reciprocal model significantly changes the identity of motifs and anti-motifs across commodities.Hence, it is essential to take into account the type of the corresponding link in which weights are sampled by constraining reciprocated and non-reciprocated strengths.
Overall, our results indicate that product-level information is strictly necessary to identify triadic structures and fluxes in production networks.We hope that our study can encourage Statistics Bureaus around the world to implement policies and techniques to reveal or reconstruct a reliable product heterogeneity for firm-level transaction data.Our analysis also shows that most commodity-specific layers can be reconstructed via null models that incorporate reciprocity while maintaining dyads independent.For these layers, network reconstruction methods of the type introduced in [26], if extended to incorporate reciprocity, are likely to perform well in replicating the properties of the entire layers starting from partial, node-specific information.Most other layers show at most one or a couple of deviating triadic motifs that are unexplained by the null model.For these layers, additional information is needed to achieve a good reconstruction.Once a rigorous product analysis has been performed, experts in the single commodity can interpret why such triadic formations over-occur or under-occur, accommodating an excessive or insufficient amount of trade volume, unveiling the detailed structure of the commodity-specific production networks.
In order to suggest improvements for further research, we conclude by noticing that our study is subject to two main limitations.First, an industry-level analysis inherently yields results that differ from those obtained from firm-level studies and underestimates the risk associated with exogenous and endogenous shocks [40].Second, the dataset analyzed pertains to industries at the SBI5 level, a classification intermediate between firms and SBI4-level industries.Analyzing the dataset at the firm level was not feasible due to potential biases arising from the deterministic imputation and degree distribution assumptions, which are exacerbated when dealing with highly granular data.Conversely, analyzing industries at the SBI4 level, which encompasses a maximum of 132 industries, would imply that for a substantial number of commodities, very few industries are active.Consequently, the null model would trivially replicate, in a statistical sense, the triadic structures for the majority of commodity layers due to a lack of relevant observations.However, the same biases anticipated at the firm level can arise, even if mitigated, by selecting SBI5level industries.This could potentially lead to biases in our analysis, especially in the type of motifs and anti-motifs found for each commodity.However, in order to validate all of our 'fingerprints' we would need fully empirical data for industries at the SBI5 level for each of the 187 commodities, an information that is not available in any country until now, to the best of our knowledge.

Binary Null Models
For binary-directed graphs, the Maximum Entropy formalism prescribes the maximization of the Graph Entropy functional S[P (A)] subject to the normalization of the Graph Probability P (A) and to the constraints on network properties C * α , i.e.
A∈A P (A) = 1 hence maximizing the unbiasedness of the resulting P (A) given available data.Solving the optimization problem, we obtain the canonical P (A) where H(A) is denoted as the Graph Hamiltonian and is defined as In this section, we focus on the binary reconstruction methods taking into account local properties.
The Directed Binary Configuration Model In the Directed Binary Configuration Model (DBCM), we choose as local properties the the out-degree (k out i ) and the in-degree (k in i ) representing the number of industries industry i sells to and the number of industries industry i buys from respectively.
Out-degrees and in-degrees can be defined mathematically in terms of the adjacency matrix A = (a ij ) as Solving the constrained Entropy maximization we obtain the Graph Probability P (A) in Eq. ( 8) where The Graph Probability P (A) can be re-written as the product of Bernoulli trials where p ij = P (a ij 1) denotes the probability of connection of supplier i with user j and is equal to with x out i ≡ e −α out i and x in i ≡ e −α in i .By Maximum Log-Likelihood Estimation (MLE) on the log-likelihood L = ln(P (A)) we obtain the Lagrange parameters α out i and α in i ∀i, a procedure equivalent to solving a system of 2N coupled equations where N is the number of industries in the network and ⟨k out i ⟩ and ⟨k in i ⟩ denote the ensemble averages of out-degrees and in-degrees respectively.
The Reciprocal Binary Configuration Model In the Reciprocal Binary Configuration Model (RBCM), we decompose the degree according to the reciprocal nature of the connection at hand, namely in non-reciprocated out-degree k → i , non-reciprocated in-degree k ← i and reciprocated degree k ↔ i .Those measures can be defined mathematically in terms of the adjacency matrix A = (a ij ) as Solving the Constrained Maximization Entropy problem, we obtain the Graph Probability P (A) as in Eq. ( 8) with Graph Hamiltonian given by The model-induced Graph Probability P (A) is the product of Bernoulli trials of mutually exclusive events with where are the exponentiated Lagrange multipliers tuning for the nonreciprocated out-degree, non-reciprocated in-degree and reciprocated degree respectively.The Lagrange multipliers α → i , α ← i and α ↔ i are found using MLE on the Log-likelihood L = ln(P a procedure equivalent to solving the system of 3N coupled equations reading i.e., equating the reciprocated and non-reciprocated degrees to their ensemble averages.

Conditional Weighted Null Models
When inspecting network weights, the numeric character of the involved trade volumes restricts the basket of available models.If the weights are discrete-valued, the constrained Entropy maximization leads to a family of geometric distributions [66,70,78].In contrast, continuous values lead to a family of exponential probability distributions when the constraints arise from node-specific properties [68,71].We treat the conditional problem, which is well defined only after deciding the form of the binary adjacency matrix A.
The conditional Graph Entropy S[Q(W |A)], measuring the uncertainty attached to the probability of having a weighted adjacency matrix W compatible with a given realization of the binary adjacency matrix A, i.e.

S[Q(
is maximized given the normalization of the conditional weighted probability density function Q(W |A) and the constraints C α (W ) where the set of C * α represent known node-specific properties.From this constrained conditional maximization we obtain Q(W |A), as where W A stands for the ensemble of realizations of W compatible with A (with weights sampled only on connected dyads a ij = 1) and the Graph Hamiltonian H(W ) is defined as Parameters β α are estimated using MLE on the log-likelihood function L W reading where Z ⃗ β,A is the conditional partition function and its computation is possible only if total information about A is available.However, estimating parameters on the empirical topology A neglects its intrinsic random variability when it is sampled using a binary model, such as DBCM or RBCM.This problem is solved in Network Science literature by defining the generalized log-likelihood G ⃗ β [68,80] where P (A) is the Graph Probability induced by the binary model.In the following, we mainly deploy the estimation based on G ⃗ β for weighted models.Using the framework mentioned above, we can solve the conditional maximum Entropy problem taking into account weighted local properties.
The CReM A When randomizing the weighted adjacency matrix W , trade marginals such as the out-strength s out i and the in-strength s in i -representing the total output or total input of industry i -are usually constrained [66,67].The out-strength s out i and the in-strength s in i sequences are defined as the marginals of the weighted adjacency matrix W , namely Solving the constrained conditional Entropy maximization leads to a conditional cumulative function Q(W |A) as in Eq. ( 22) where with a conditional Graph distribution i.e. the product of dyadic exponential distributions in w ij conditional on the establishment of the link a ij and regulated by the node-specific Lagrange parameters β out i and β in i ∀i.By using Generalized Log-likelihood Estimation (GLE), we find the Lagrange parameters -a procedure that equates to slightly changing the dyadic conditional probability by substituting a ij with a dyadic term f ij such that f ij = ⟨a ij ⟩, i.e., f ij is the ensemble average of a ij and Maximizing G ⃗ β we obtain a system of 2N coupled equations reading and find {β in i , β out i } for each industry.The CRWCM model In order to take into account reciprocity, we develop a novel model denoted as Conditionally Reciprocal Weighted Configuration Model (CRWCM), that considers the different nature of links on which weights are sampled, namely reciprocated and non-reciprocated links.This choice leads to the definition of four trade marginals for each supplier/user, namely • the non-reciprocated out-strength s → i which measures the output of supplier i to users from which it does not buy, defined in terms of W as • the non-reciprocated in-strength s ← i , which measures the input of industry i from suppliers to which it does not itself supply, defined as • the reciprocated out-strength s ↔,out i , measuring the output of supplier i to users from which it also purchases, reading • and the reciprocated in-strength s ↔,in i , measuring the input of user i from suppliers to which it also supplies, defined as Solving the constrained conditional Maximum Entropy problem, we obtain the Conditional Weighted Graph Probability in Eq. ( 22) where the Graph Hamiltonian is given by leading to q ij (w|a ij = 1) = (36) where q ij (w|a ij ) for the single dyad depends on the possible states of w ij , namely Rephrasing the vector {a → ij , a ← ij , a ↔ ij , a ̸ ↔ ij } of a ij -states into the vector of their ensemble averages for the non-reciprocated sub-problem and for the reciprocated sub-problem.

FIG. 1 .
FIG. 1.(a) Graphical representation of the Dutch multi-layer production network.For illustrative purposes, we represent three industries/firms i, j and k as nodes, all placed on three commodity group layers, namely (from top to bottom): Cereals, Beer/Malt, Bread and other bakery products.The connections between the same three nodes are different in the different layers.(b) The possible 13 types of connected triadic subgraphs.Each triple of Industries/Firms can trade different products by forming, on each commodity-specific layer, either one of the 13 possible connected subgraphs or one of the remaining subgraphs where at least one node is disconnected (not shown).

FIG. 2 .
FIG. 2. Normalized Triadic Occurrences (a) and Fluxes (b): the Aggregated Network (− • −) presents a high occurrence of subgraphs m = 1 and m = 13, representing open-Vs and completely reciprocated triads, respectively.The latter covers most of the total amount of money traded.The Cereals commodity layer (− • −), with a high occurrence of subgraph m = 1.A relatively high amount of money is distributed across m = 1, m = 4 and m = 6.Gas/Hot Water/City Heating layer (− • −) with a predominant occurrence and flux in subgraph m = 1.Agricultural Services layer (− • −), with a highly heterogenous spectrum of occurrences and fluxes.Completely cyclical triads have a high occurrence in the aggregated network, but break apart when passing to single commodity layers as G.H.C and Cereals, if not for rare cases such as Agricultural Services.In single commodity layers m = 1 receives the highest concentration of money, signalling a large amount of money flows over structures that greatly depend on a limited number of suppliers, which control the market.

FIG. 3 .
FIG. 3. Triadic binary motif analysis: DBCM (•) vs RBCM (•).(a) Analysis of the aggregated network with a single representative commodity.Numerous motifs and anti-motifs are present using DBCM and RBCM as null models.(b-d) Commodity groups where RBCM reproduces all the triadic structures, and they are, respectively, Cereals, Electrical Components, and the Construction of Tunnels, Waterways, and Roads.(e-f) Commodity groups with one network motif, namely Bread and Gasoline.(g) Commodity group with two network motifs, namely Beer/Malt.The CIs are computed by extracting the 2.5-th and 97.5-th percentile from an ensemble distribution of 500 graphs.The numerous motifs and anti-motifs in the aggregated network can be seen as the aggregation of commodity groups presenting very few characteristic patterns.

FIG. 5 .
FIG. 5. Triadic weighted motif analysis: DBCM+CReMA (•) vs RBCM+CRWCM (•).(a) Analysis of the aggregated network with a single representative commodity.A large number of motifs and anti-motifs are present when using DBCM+CReMA, while three motifs are present when using the RBCM+CRWCM.(b-d) Commodity groups where RBCM+CRWCM reproduces all the triadic structures, and they are, respectively, Seeds, Metal Components for Doors & Windows, and Airline Services.(e-f) Commodity groups with one network motif, namely Coffee/Tea and Textile raw materials and products.(g) Commodity group with two network motifs, namely Shipping Services.The CIs are computed by extracting the 2.5-th and 97.5-th percentile from an ensemble distribution of 500 graphs.Passing from the aggregated network to the disaggregated product layers unveils the presence of a few commodity-specific motifs and anti-motifs.
j̸ =i;aij =1 q ij (w|a = 1 ij ⟩ depends on the binary model of choice, we can use GLE for the estimation of the 4N parameters.The resulting generalized log-likelihood is separable in a reciprocal and non-reciprocal component, i.e., G ⃗ β = G → ⃗ β + G ↔ ⃗ β (see Appendix B for details).The Lagrange parameters ⃗ β are retrieved by maximizing G ⃗ β , which equates to solving two uncoupled systems of 2N coupled equations reading