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

Marzio Di Vece<EMAIL_ADDRESS>IMT School for Advanced Studies
Lucca, P.zza San Francesco 19, 55100 Lucca (Italy)  Lorentz Institute for
Theoretical Physics, Leiden University, Niels Bohrweg 2, 2333CA Leiden (The
Netherlands)  Frank P. Pijpers  Statistics Netherlands, Henri Faasdreef 312,
2492 JP Den Haag (the Netherlands)  Korteweg - de Vries Institute for
Mathematics, University of Amsterdam, Amsterdam (the Netherlands)  Diego
Garlaschelli  IMT School for Advanced Studies Lucca, P.zza San Francesco 19,
55100 Lucca (Italy)  Lorentz Institute for Theoretical Physics, Leiden
University, Niels Bohrweg 2, 2333CA Leiden (The Netherlands)  INdAM-GNAMPA
Istituto Nazionale di Alta Matematica (Italy)

###### Abstract

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 analyse both
binary and weighted triadic motifs in the Dutch inter-industry production
network disaggregated at the level of 187 commodity groups, using data from
Statistics Netherlands. 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 80% 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.

###### pacs:

89.75.Fb; 02.50.Tt; 89.65.Gh

##  I Introduction

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

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 [7, 8, 9, 10, 11, 12, 14, 13, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]
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 [17, 28, 6], 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]. 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.
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) produced 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 presence of product granularity makes it an invaluable source for
the analysis of commodity-specific structural patterns.

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. Firms are matched according to a deterministic
procedure, which is shown to decrease the dataset quality by inducing a bias
in the network density and the degree distribution, as proved in [37] for a
sample of known links of the production network collected by Dun & Bradstreet.
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
[44], playing the role of functional modules or evolutionary signs in
biological networks [45, 46], homophily-driven connections in social networks
[47], complementarity-driven structures in production networks [36, 33], their
change in time being interpreted as self-organizing processes in the World
Trade Web (WTW) [48, 49], and early-warning signals of topological collapse in
inter-bank networks [50, 51] and stock market networks [42]. It has been
proven that for the majority of (available) real-world networks, the triadic
structure is maximally random [52] and by fixing it their global structure is
statistically determined [53].

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 [54].

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 [55, 56, 57] 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) [58]. This theoretical framework has been proven
to successfully reconstruct economic and financial systems [59, 60, 61, 62],
accurately predicting both the topology and the weights of the WTW [63, 64,
65], in an integrated [66], or conditional fashion [67], with only structural
constraints, or informing the models with economic factors [68, 69, 70],
accurately predicting banks’ risk exposures [71], and most recently,
accurately reconstructing payment flows among Dutch firms that were clients of
ABN Amro Bank or ING Bank, constraining their industry-specific production
functions [26].

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 [72], and an applied study where triadic motifs and their time
evolution are used as early warnings of topological collapse during the 2008
financial crisis [50].

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 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, they can
eventually advance policy laws to nudge a more convenient redistribution of
connections and trade volumes.

The rest of the paper is structured as follows. In Sec. II, we describe the
inference procedure and the limitations of the CBS dataset and the pre-
processing performed for the generation of the inter-industry production
network. In Sec. III, we describe the Maximum-Entropy randomization methods
used for triadic pattern detection in their binary (in Sec. III.1.) and
weighted declinations (in Sec. III.2.). Then, in Sec. IV. we perform a
descriptive analysis of the empirical inter-industry network, followed by a
Binary Motif Analysis in Sec. IV.A. and a Weighted Motif Analysis in Sec.
IV.B. In Sec. V., we discuss our findings.

##  II 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. Firm-level data is
obtained from the General Business Register (ABR) for 2018, containing data
for over $1\,700\,000$ firms. After cleaning for micro-firms with annual net
turnover below $10\,000$€, around $900\,000$ firms remain, accounting for
$99.5\%$ of the Dutch economy output in 2018. The breakdown in commodities is
extracted from the Structural Business Statistics (SBS) survey for commercial
industries, from the Prodcom survey for manufacturing industries, and
estimated by National Accounts for non-commercial industries. In most cases,
the commodity breakdown data is available for the industry as a whole and not
for individual organizations within industries. A breakdown of intermediate
supply and use per firm follows using intermediate purchases as the
distributional key, and then an ulterior breakdown in commodity groups is
performed by integrating data from the SBS and Prodcom surveys. The resulting
dataset is then compared to the industry-level supply/use tables at the SBI4
classification, and appropriate rescaling of supply and use per firm is
performed by Iterative Proportion Fitting. Once supply and use per firm per
commodity are obtained, their in-degree distribution is estimated using
stylized facts from [39] Japanese firms. Suppliers and users are then matched
according to a deterministic procedure that takes into account (1) their trade
capacity, encoded by their net turnover, (2) their mutual distance, (3) the
presence of a link between respective industries, (4) the observed
relationship in the Dun & Bradstreet dataset, which contains the customers of
the largest 500 suppliers in the Dutch Economy. Finally, the resulting inter-
firm network at the $650$ commodity level is compared to the inter-industry
(known) network at the $250$ commodity level and consequent adjustments are
done to weights and links.
Even if a consistent number of auxiliary micro-data is available, the
intensive imputation procedure and the known biases contained in the previous
version [37] makes it desirable to reduce the bias by appropriately
transforming the dataset. Consequently, we take advantage of the tested
coherence between the inter-firm network for $192$ commodities and the inter-
industry network as the key point of our pre-processing. Specifically, we
aggregate the $650$ commodity groups into $192$, coherently with industry-
level known tables, and then we aggregate firms according to their SBI5
Standard Industry Classification extracted from the ABR for 2018. Passing from
the SBI4 to the SBI5 classification leads to a greater granularity in industry
resolution, increasing the number of industries from 132 to 870. Finally, we
clean for intra-industry trade and obtain a multi-layer inter-industry
production network containing linkages and weights for $862$ industries
(nodes) and $187$ commodity groups (layers). For the topic of interest,
triadic motifs, the self-loops implied by intra-industry trade are not
important and can be removed from the dataset without adversely affecting the
subsequent analysis. Therefore, although intra-industry trade is certainly
relevant for both intensive and extensive margins (weights and links
respectively) they are ignored for the rest of this paper.
The resulting dataset, while being the most reliable and detailed multi-
product inter-industry domestic production network for intermediate use in
existence to our knowledge, has clear limitations: (1) the lack of (known)
firm granularity, which is a necessary ingredient to unveil more detailed
network anomalies, especially in intra-industry supply/use relationships, (2)
the possible presence of imputation bias in the disaggregation from the SBI4
to the SBI5 Industry Classification. Nonetheless, we believe that our findings
shed light on the role of product granularity in the modeling of production
networks and the characterization of their triadic structures.

##  III 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 which give the best insurance
of unbiasedness with respect to missing data, as proven by independent testing
[73, 74, 75, 76]. Statistical measures of interest are then extracted as
ensemble averages.

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. If the empirical networks’ statistics show
significant deviations from the model-induced ensemble averages, they are a
sign of higher-order correlations not explained by the structural constraints.

Reliable data is available for $187$ commodity groups and not for the totality
of specific commodities, implying that a commodity group can contain distinct,
but similar, specific products. A supplier of one of the products belonging to
the group can be simultaneously a user of a different product in the same
group, leading to possible reciprocated links among industries in the
commodity group layer. This leads us, in Sec. III.1, to also treat models that
encode not only the tendencies to supply and use but also the tendency to
reciprocate their supply/use relationship in the single commodity group.
Finally, in Sec. III.2, we employ models to predict the expected trade volumes
among industries in their usual declination or distinguish them according to
the reciprocal nature of the corresponding links.

###  III.1 Binary Null Models

For binary-directed graphs, the Maximum Entropy formalism prescribes the
maximization of the Graph Entropy functional $S[P(A)]$

$S[P(A)]=-\sum_{A\in\mathbf{A}}P(A)\ln P(A)$ (1)

subject to the normalization of the Graph Probability $P(A)$ and to the
constraints on network properties $C_{\alpha}^{*}$, i.e.

$\begin{cases}\sum_{A\in\mathbf{A}}P(A)&=1\\\
\sum_{A\in\mathbf{A}}P(A)C_{\alpha}(A)&=C_{\alpha}^{*},\quad\forall\alpha,\end{cases}$
(2)

hence maximizing the unbiasedness of the resulting $P(A)$ given available
data. Solving the optimization problem, we obtain the canonical $P(A)$

$\displaystyle P(A)$
$\displaystyle=\dfrac{e^{-\sum_{\alpha}\theta_{\alpha}C_{\alpha}^{*}(A)}}{\sum_{A\in\mathbf{A}}e^{-\sum_{\alpha}\theta_{\alpha}C_{\alpha}^{*}(A)}}=$
$\displaystyle=\dfrac{e^{-H(A)}}{\sum_{A\in\mathbf{A}}e^{-H(A)}}$ (3)

where $H(A)$ is denoted as the Graph Hamiltonian and is defined as

$H(A)\equiv\sum_{\alpha}\theta_{\alpha}C_{\alpha}(A)^{*}.$ (4)

In this section, we focus on the binary reconstruction methods taking into
account local properties.

####  III.1.1 The Directed Binary Configuration Model

In the Directed Binary Configuration Model (DBCM), we choose as local
properties the the out-degree $(k_{i}^{out})$ and the in-degree $(k_{i}^{in})$
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

$\begin{cases}k_{i}^{out}&=\sum_{j\neq i}a_{ij}\\\ k_{i}^{in}&=\sum_{j\neq
i}a_{ji}.\end{cases}$ (5)

Solving the constrained Entropy maximization we obtain the Graph Probability
$P(A)$ in Eq. (III.1) where

$H(A)=\sum_{i}\alpha_{i}^{out}k_{i}^{out}+\alpha_{i}^{in}k_{i}^{in}.$ (6)

The Graph Probability $P(A)$ can be re-written as the product of Bernoulli
trials

$P(A)=\prod_{i,j\neq i}(p_{ij})^{a_{ij}}(1-p_{ij})^{1-a_{ij}}$ (7)

where $p_{ij}=P(a_{ij}=1)$ denotes the probability of connection of supplier
$i$ with user $j$ and is equal to

$p_{ij}=\dfrac{x_{i}^{out}x_{j}^{in}}{1+x_{i}^{out}x_{j}^{in}}$ (8)

with $x_{i}^{out}\equiv e^{-\alpha_{i}^{out}}$ and $x_{i}^{in}\equiv
e^{-\alpha_{i}^{in}}$. By Maximum Log-Likelihood Estimation (MLE) on the log-
likelihood $\mathcal{L}=\ln(P(A))$ we obtain the Lagrange parameters
$\alpha_{i}^{out}$ and $\alpha_{i}^{in}$ $\forall i$, a procedure equivalent
to solving a system of $2N$ coupled equations

$\begin{cases}k_{i}^{out,*}&=\langle k_{i}^{out}\rangle=\sum_{j\neq
i}p_{ij}\\\ k_{i}^{in,*}&=\langle k_{i}^{in}\rangle=\sum_{j\neq
i}p_{ji}.\end{cases}$ (9)

where $N$ is the number of industries in the network and $\langle
k_{i}^{out}\rangle$ and $\langle k_{i}^{in}\rangle$ denote the ensemble
averages of out-degrees and in-degrees respectively.

####  III.1.2 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}^{\rightarrow}$, non-reciprocated in-degree
$k_{i}^{\leftarrow}$ and reciprocated degree $k_{i}^{\leftrightarrow}$. Those
measures can be defined mathematically in terms of the adjacency matrix
$A=(a_{ij})$ as

$\begin{cases}k_{i}^{\rightarrow}&=\sum_{j\neq i}a_{ij}(1-a_{ji})=\sum_{j\neq
i}a_{ij}^{\rightarrow}\\\ k_{i}^{\leftarrow}&=\sum_{j\neq
i}a_{ji}(1-a_{ij})=\sum_{j\neq i}a_{ij}^{\leftarrow}\\\
k_{i}^{\leftrightarrow}&=\sum_{j\neq i}a_{ij}a_{ji}=\sum_{j\neq
i}a_{ij}^{\leftrightarrow}.\end{cases}$ (10)

Solving the Constrained Maximization Entropy problem, we obtain the Graph
Probability $P(A)$ as in Eq. (III.1) with Graph Hamiltonian given by

$H(A)=\sum_{i}\alpha_{i}^{\rightarrow}k_{i}^{\rightarrow}+\alpha_{i}^{\leftarrow}k_{i}^{\leftarrow}+\alpha_{i}^{\leftrightarrow}k_{i}^{\leftrightarrow}.$
(11)

The model-induced Graph Probability $P(A)$ is the product of Bernoulli trials
of mutually exclusive events

$P(A)=\prod_{j<i}\left(p_{ij}^{\rightarrow}\right)^{a_{ij}^{\rightarrow}}\left(p_{ij}^{\leftarrow}\right)^{a_{ij}^{\leftarrow}}\left(p_{ij}^{\leftrightarrow}\right)^{a_{ij}^{\leftrightarrow}}\left(p_{ij}^{\nleftrightarrow}\right)^{a_{ij}^{\nleftrightarrow}}$
(12)

with

$\begin{cases}p_{ij}^{\rightarrow}&=\dfrac{x_{i}^{\rightarrow}x_{j}^{\leftarrow}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}\\\
p_{ij}^{\leftarrow}&=\dfrac{x_{i}^{\leftarrow}x_{j}^{\rightarrow}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}\\\
p_{ij}^{\leftrightarrow}&=\dfrac{x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}\\\
p_{ij}^{\nleftrightarrow}&=\left[1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}\right]^{-1}.\end{cases}$
(13)

where $x_{i}^{\rightarrow}\equiv e^{-\alpha_{i}^{\rightarrow}}$,
$x_{i}^{\leftarrow}\equiv e^{-\alpha_{i}^{\leftarrow}}$ and
$x_{i}^{\leftrightarrow}\equiv e^{-\alpha_{i}^{\leftrightarrow}}$ are the
exponentiated Lagrange multipliers tuning for the non-reciprocated out-degree,
non-reciprocated in-degree and reciprocated degree respectively. The Lagrange
multipliers $\alpha_{i}^{\rightarrow}$, $\alpha_{i}^{\leftarrow}$ and
$\alpha_{i}^{\leftrightarrow}$ are found using MLE on the Log-likelihood
$\mathcal{L}=\ln(P(A))$, a procedure equivalent to solving the system of $3N$
coupled equations reading

$\begin{cases}k_{i}^{\rightarrow}&=\langle
k_{i}^{\rightarrow}\rangle=\sum_{j\neq i}p_{ij}^{\rightarrow}\\\
k_{i}^{\leftarrow}&=\langle k_{i}^{\leftarrow}\rangle=\sum_{j\neq
i}p_{ij}^{\leftarrow}\\\ k_{i}^{\leftrightarrow}&=\langle
k_{i}^{\leftrightarrow}\rangle=\sum_{j\neq
i}p_{ij}^{\leftrightarrow},\end{cases}$ (14)

i.e., equating the reciprocated and non-reciprocated degrees to their ensemble
averages.

###  III.2 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 [69, 65, 77]. In contrast, continuous values lead to a family of
exponential probability distributions when the constraints explicate node-
specific properties [67, 70]. 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(W|A)]=-\sum_{A\in\mathbf{A}}P(A)\int_{W_{A}}Q(W|A)\ln Q(W|A)dW$ (15)

is maximized given the normalization of the conditional weighted probability
density function $Q(W|A)$ and the constraints $C_{\alpha}(W)$

$\begin{cases}\int_{W_{A}}Q(W|A)dW&=1\\\
\sum_{A}P(A)\int_{W_{A}}Q(W|A)C_{\alpha}(W)dW&=C_{\alpha}^{*},\quad\forall\alpha.\end{cases}$
(16)

where the set of $C_{\alpha}^{*}$ represent known node-specific properties.
From this constrained conditional maximization we obtain $Q(W|A)$, as

$Q(W|A)=\begin{cases}\dfrac{e^{-H(W)}}{\int_{W_{A}}e^{-H(W)}dW_{A}}\quad&W\in
W_{A}\\\ 0&W\notin W_{A}\end{cases}$ (17)

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

$H(W)\equiv\sum_{\alpha}\beta_{\alpha}C_{\alpha}(W).$ (18)

Parameters $\beta_{\alpha}$ are estimated using MLE on the log-likelihood
function $\mathcal{L}_{W}$ reading

$\mathcal{L}_{W|A}=-H_{\vec{\beta}}(W)-\ln(Z_{\vec{\beta},A})$ (19)

where $Z_{\vec{\beta},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 $\mathcal{G_{\vec{\beta}}}$ [67]

$\mathcal{G_{\vec{\beta}}}=-H_{\vec{\beta}}(\langle
W\rangle)-\sum_{A\in\mathbf{A}}P(A)\ln(Z_{\vec{\beta},A})$ (20)

where $P(A)$ is the Graph Probability induced by the binary model. In the
following, we mainly deploy the estimation based on
$\mathcal{G_{\vec{\beta}}}$ for weighted models. Using the framework mentioned
above, we can solve the conditional maximum Entropy problem taking into
account weighted local properties.

####  III.2.1 The CReMa

When randomizing the weighted adjacency matrix $W$, trade marginals such as
the out-strength $s_{i}^{out}$ and the in-strength $s_{i}^{in}$ \- explicating
the total output or total input of industry $i$ \- are usually constrained
[65, 66]. The out-strength $s_{i}^{out}$ and the in-strength $s_{i}^{in}$
sequences are defined as the marginals of the weighted adjacency matrix $W$,
namely

$\begin{cases}s_{i}^{out}&=\sum_{j\neq i}w_{ij}\\\ s_{i}^{in}&=\sum_{j\neq
i}w_{ji}.\end{cases}$ (21)

Solving the constrained conditional Entropy maximization leads to a
conditional cumulative function $Q(W|A)$ as in Eq. (17) where

$H(W)=\sum_{i}\beta_{i}^{out}s_{i}^{out}+\beta_{i}^{in}s_{i}^{in}$ (22)

with a conditional Graph distribution

$\displaystyle Q(W|A)$ $\displaystyle=\prod_{i,j\neq
i;a_{ij}=1}q_{ij}(w|a=1)=$ $\displaystyle=\prod_{i,j\neq
i;a_{ij}=1}\left[\left(\beta_{i}^{out}+\beta_{j}^{in}\right)e^{-(\beta_{i}^{out}+\beta_{j}^{in})w_{ij}}\right]^{a_{ij}}$
(23)

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 $\beta_{i}^{out}$ and $\beta_{i}^{in}$ $\forall 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}=\langle a_{ij}\rangle$, i.e., $f_{ij}$ is the ensemble
average of $a_{ij}$ and

$q_{ij}(w_{ij}|a_{ij}=1)=\left[\left(\beta_{i}^{out}+\beta_{j}^{in}\right)e^{-(\beta_{i}^{out}+\beta_{j}^{in})}\right]^{f_{ij}}.$
(24)

Maximizing $G_{\vec{\beta}}$ we obtain a system of $2N$ coupled equations
reading

$\begin{cases}s_{i}^{out}&=\sum_{j\neq
i}\dfrac{f_{ij}}{\beta_{i}^{out}+\beta_{j}^{in}}=\langle s_{i}^{out}\rangle\\\
s_{i}^{in}&=\sum_{j\neq
i}\dfrac{f_{ji}}{\beta_{i}^{in}+\beta_{j}^{out}}=\langle
s_{i}^{in}\rangle\end{cases}$ (25)

and find $\\{\beta_{i}^{in},\beta_{i}^{out}\\}$ for each industry.

####  III.2.2 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}^{\rightarrow}$ which measures the
output of supplier $i$ to users from which it does not buy, defined in terms
of $W$ as

$s_{i}^{\rightarrow}=\sum_{j\neq i}a_{ij}^{\rightarrow}w_{ij}=\sum_{j\neq
i}w_{ij}^{\rightarrow}$ (26)

  * •

the non-reciprocated in-strength $s_{i}^{\leftarrow}$, which measures the
input of industry $i$ from suppliers to which it does not supply, defined as

$s_{i}^{\leftarrow}=\sum_{j\neq i}a_{ij}^{\leftarrow}w_{ji}=\sum_{j\neq
i}w_{ij}^{\leftarrow}$ (27)

  * •

the reciprocated out-strength $s_{i}^{\leftrightarrow,out}$, measuring the
output of supplier $i$ to users from which it also buys from, reading

$s_{i}^{\leftrightarrow,out}=\sum_{j\neq
i}a_{ij}^{\leftrightarrow}w_{ij}=\sum_{j\neq i}w_{ij}^{\leftrightarrow,out}$
(28)

  * •

and the reciprocated in-strength $s_{i}^{\leftrightarrow,in}$, measuring the
input of user $i$ from suppliers to which it also supplies to, defined as

$s_{i}^{\leftrightarrow,in}=\sum_{j\neq
i}a_{ij}^{\leftrightarrow}w_{ji}=\sum_{j\neq i}w_{ij}^{\leftrightarrow,in}$
(29)

Solving the constrained conditional Maximum Entropy problem, we obtain the
Conditional Weighted Graph Probability in Eq. (17) where the Graph Hamiltonian
is given by

$H(W)=\sum_{i}\beta_{i}^{\rightarrow}s_{i}^{\rightarrow}+\beta_{i}^{\leftarrow}s_{i}^{\leftarrow}+\beta_{i}^{\leftrightarrow,out}s_{i}^{\leftrightarrow,out}+\beta_{i}^{\leftrightarrow,in}s_{i}^{\leftrightarrow,in}$
(30)

leading to

$\displaystyle Q(W|A)$ $\displaystyle=\prod_{j\neq
i,a_{ij}=1}q_{ij}(w|a_{ij}=1)=$

where $q_{ij}(w|a_{ij})$ for the single dyad depends on the possible states of
$w_{ij}$, namely

$\begin{cases}(\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow})e^{-(\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow})w_{ij}^{\rightarrow}}&\text{for
$w_{ij}^{\rightarrow}>0$}\\\
(\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in})e^{-(\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in})w_{ij}^{\leftrightarrow,out}}&\text{for
$w_{ij}^{\leftrightarrow,out}>0$}\\\ 0&\text{for $w_{ij}=0$}.\end{cases}$ (32)

Rephrasing the vector
$\\{a_{ij}^{\rightarrow},a_{ij}^{\leftarrow},a_{ij}^{\leftrightarrow},a_{ij}^{\not{\leftrightarrow}}\\}$
of $a_{ij}$-states into the vector of their ensemble averages
$\\{f_{ij}^{\rightarrow},f_{ij}^{\leftarrow},f_{ij}^{\leftrightarrow},f_{ij}^{\not{\leftrightarrow}}\\}$,
where $f_{ij}^{(\cdot)}=\langle a_{ij}^{(\cdot)}\rangle$ 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.,
$\mathcal{G}_{\vec{\beta}}=\mathcal{G}_{\vec{\beta}}^{\rightarrow}+\mathcal{G}_{\vec{\beta}}^{\leftrightarrow}$
(see Appendix B for details). The Lagrange parameters $\vec{\beta}$ are
retrieved by maximizing $\mathcal{G}_{\vec{\beta}}$, which equates to solving
two uncoupled systems of $2N$ coupled equations reading

$\begin{cases}s_{i}^{\rightarrow}&=\sum_{j\neq
i}\dfrac{f_{ij}^{\rightarrow}}{\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow}}=\langle
s_{i}^{\rightarrow}\rangle\\\ s_{i}^{\leftarrow}&=\sum_{j\neq
i}\dfrac{f_{ij}^{\leftarrow}}{\beta_{i}^{\leftarrow}+\beta_{j}^{\rightarrow}}=\langle
s_{i}^{\leftarrow}\rangle\end{cases}$ (33)

for the non-reciprocated sub-problem and

$\begin{cases}s_{i}^{\leftrightarrow,out}&=\sum_{j\neq
i}\dfrac{f_{ij}^{\leftrightarrow}}{\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in}}=\langle
s_{i}^{\leftrightarrow,out}\rangle\\\ s_{i}^{\leftrightarrow,in}&=\sum_{j\neq
i}\dfrac{f_{ij}^{\leftrightarrow}}{\beta_{i}^{\leftrightarrow,in}+\beta_{j}^{\leftrightarrow,out}}=\langle
s_{i}^{\leftrightarrow,in}\rangle\end{cases}$ (34)

for the reciprocated sub-problem.

Figure 2: Triadic Occurrences and Fluxes: (a) the Aggregate 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. (b) 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$. (c) Gas/Hot Water/City Heating layer with a
predominant occurrence and flux in subgraph $m=1$. (d) Agricultural Services
layer, with a highly heterogenous spectrum of occurrences and fluxes.
Completely cyclical triads have a high occurrence in the aggregate picture (a)
but break apart when passing to single commodity layers (b and c) if not for
rare cases (d).

##  IV Results

Product granularity gives us the opportunity to study heterogeneity across
commodity layers. Let us consider the number of layer-active industries $N$
and reciprocity measures such as the topological reciprocity $r_{t}$, defined
as the ratio of reciprocated links to $L$, i.e.

$r_{t}=\dfrac{L^{\leftrightarrow}}{L}=\dfrac{\sum_{i,j\neq
i}a_{ij}^{\leftrightarrow}}{\sum_{i,j\neq i}a_{ij}}.$ (35)

and its weighted counterpart $r_{w}$, defined as the ratio of total weight on
reciprocated links to $W$, i.e.

$r_{w}=\dfrac{W_{tot}^{\leftrightarrow}}{W_{tot}}=\dfrac{\sum_{i,j\neq
i}w_{ij}^{\leftrightarrow,out}}{\sum_{i,j\neq i}w_{ij}}.$ (36)

For around $50\%$ of the layers, $N$ is below 100. For another $30\%$ of
layers $N$ is between 100 and 600, while for the last $20\%$, it is between
600 and 818. There are few reciprocated connections ($r_{t}<0.05$) for half of
the commodity layers, while large reciprocity $r_{t}>0.4$ is present for
around $8\%$ of the layers. Analogously the weights on reciprocated links
represent a minority share for half of the commodity groups ($r_{w}<0.08$),
while goods are traded in large part on reciprocated links ($r_{w}>0.5$) for
$10\%$ of layers. 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 large industry heterogeneity. Some
examples are suppliers of plastic goods who sell to users with heterogeneous
specializations, for instance, Bread, Beer, Cereals, Fish, etc. In addition,
there is high binary reciprocity only for a small number of layers but in
those connections a high concentration of money is present.

Let us now move to the analysis of triads. Consider in Fig. 2 triadic
occurrences $N_{m}$, defined as the number of times a specific m-subgraph
appears and triadic fluxes $F_{m}$, defined as the total amount of money
circulating on each m-subgraph. They are shown in Fig. 2(a) for the aggregated
network, with a single representative commodity, and in Fig. 2(b-d) for three
commodity layers, namely ‘Cereals’, ‘Gas/ Hot Water /City heating’ and
‘Agricultural Services’. In the Aggregated Network, occurrences of m-subgraphs
of type $m=1$ and $m=13$ are predominant and the majority share of money
circulates on $m=13$. Disaggregating in product layers leads to a decrease of
occurrences of type $m=13$ in both Cereals and Gas/ Hot Water /City heating,
with type $m=1$ remaining predominant, i.e. completely cyclical triads break
up in favor of open triangles. However, $N_{13}$ is still high in a small
number of commodity layers. This is the case for Agricultural Services, where
$m=13$ is the second most frequent subgraph. Conversely, its weighted
counterpart $F_{13}$ is lower than $F_{1}$ and $F_{11}$, i.e. in single
product layers even if $m=13$ does not decrease in occurrence, it decreases in
terms of money circulation.

Figure 3: Triadic binary motif analysis: DBCM
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bullet}$)
vs RBCM
(${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}$).
(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. Figure 4: Comparison DBCM
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bullet}$)
vs. RBCM
(${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}$):
(a) Empirical Counter Cumulative Distribution Function $ECCDF$ representing
the number of commodities that have at least $n_{m}$ deviating triadic motifs.
(b) Number of commodities $c_{h}(m)$ having a $m$-type motif. (c) Number of
commodities $c_{l}(m)$ having a $m$-type anti-motif. RBCM explains more
triadic structures than DBCM, as shown in the difference of their
$ECCDF(n_{m})$. Passing from DBCM to RBCM reduces the number of $m$ motifs
across commodities, with the exception of $m=6$, and anti-motifs, with the
exception of $m=8$.

###  IV.1 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$

$z\left[N_{m}\right]=\dfrac{N_{m}(A^{*})-\langle
N_{m}\rangle}{\sigma\left[N_{m}\right]}$ (37)

where $N_{m}(A^{*})$ is the number of occurrences of the $m$-type subgraph in
the empirical adjacency matrix, $\langle N_{m}\rangle$ is its model-induced
expected number of occurrences, and $\sigma\left[N_{m}\right]$ is the model-
induced standard deviation.

An analytical procedure [72] 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
[78]. 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}\in 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}^{\rightarrow}$, $p_{ij}^{\leftarrow}$ and
$p_{ij}^{\not{\leftrightarrow}}$, and reciprocated connection
$p_{ij}^{\leftrightarrow}$, generate a uniform random variable
$u_{ij}\in(0,1)$ and (2) establishing the appropriate links in the dyad in the
following way:

  * •

a non-reciprocated link from $i$ to $j$ if $u_{ij}\leq p_{ij}^{\rightarrow}$;

  * •

a non-reciprocated link from $j$ to $i$ if
$u_{ij}\in(p_{ij}^{\rightarrow},p_{ij}^{\rightarrow}+p_{ij}^{\leftarrow}]$;

  * •

a reciprocated link from $i$ to $j$ (and from $j$ to $i$) if
$u_{ij}\in(p_{ij}^{\rightarrow}+p_{ij}^{\leftarrow},p_{ij}^{\rightarrow}+p_{ij}^{\leftarrow}+p_{ij}^{\leftrightarrow}]$;

  * •

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. $\langle N_{m}\rangle$ and $\sigma\left[N_{m}\right]$, are the
average and standard deviation of $N_{m}$ extracted from the ensemble
distribution of $500$ realizations of $A$. After having computed
$z\left[N_{m}\right]$, we also extract the $2.5$-th and $97$-th percentiles
from the ensemble distribution of $N_{m}$ over all models and we standardize
them using Eq. (37) 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
red on the right panel. The corresponding confidence intervals at the $5\%$
percent are depicted with the same color (blue or red) 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, i.e. the information on the reciprocal structure
qualitatively changes the z-score profile.

By disaggregating from the aggregated monolayer to the multi-commodity
network, $N_{m}$ are well reproduced once you consider the tendency to
reciprocate the trade relationship for most commodities. 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_{N_{m}}$ 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 - all the ensemble values of $N_{m}$ correspond to the
empirical $N_{m}^{*}$ \- CIs are not symmetrical and the distribution of the
corresponding $N_{m}$ is trivially not Normal. 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 probability that a commodity has at least $n_{m}$
motifs or anti-motifs 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
$60\%$ 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 $22\%$ 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 anti-motifs for each type, with the exception of anti-
motif $m=8$ as depicted in Fig. 4(c).

###  IV.2 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 the triadic subgraph of type $m$. We characterize the deviation
of empirical $F_{m}$ to null models by defining the weighted z-scores as

$z\left[F_{m}\right]=\dfrac{F_{m}(W^{*})-\langle
F_{m}\rangle}{\sigma\left[F_{m}\right]}$ (38)

where $\langle F_{m}\rangle$ is the model-induced average amount of money
circulating on motif $m$ and $\sigma\left[F_{m}\right]$ represents the model-
induced standard deviation over the ensemble of network realizations.

Figure 5:  Triadic weighted motif analysis: DBCM+CReMa
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bullet}$)
vs RBCM+CRWCM
(${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}$).
(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. Figure 6:  Comparison DBCM+CReMa
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bullet}$)
vs. RBCM+CRWCM
(${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bullet}$):
(a) Empirical Counter Cumulative Distribution Function $ECCDF$ representing
the number of commodities that have at least $n_{m}$ deviating triadic fluxes.
(b) The number of commodities $c_{h}(m)$ having a $m$-type motif. (c) The
number of commodities $c_{l}(m)$ having a $m$-type anti-motif. RBCM+CRWCM
explains slightly more triadic fluxes than DBCM+CReMa, as shown in the
difference of their $ECCDF(n_{m})$. Passing from the directed to the
reciprocal model reduces the number of anti-motifs, with the exception of
$m=8$. In contrast, it changes qualitatively the motif profile, with a slight
dominance of $m=11$-type motifs when the directed model is used and a clear
dominance of $m=1$-type motifs when the reciprocal model is used.

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 CReMa, constraining the out-strength and in-strength
sequences. We denote this model as $DBCM+CReMa$ or directed model. If we deem
reciprocal information necessary, a combination of the RBCM and CRWCM should
be used. We denote this model as $RBCM+CRWCM$ or reciprocated model. 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+CReMa reduces to
(1) establishing a link between industries $i$ and $j$ when a uniform random
number $u_{ij}\in U(0,1)$ is such that $u_{ij}\leq 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 $\eta_{ij}\in U(0,1)$ such that

$F(v_{ij})=\int_{0}^{v_{ij}}q_{CReMa}(w_{ij}|a_{ij}=1)dw_{ij}=\eta_{ij},$ (39)

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_{CReMa}(w_{ij}|a_{ij}=1)$ is
substituted with $q_{CRWCM}(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 red 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. The picture is different when taking
reciprocity information in input: there are only three motifs, namely $m=1$,
$m=6$, and $m=13$ when the reciprocal null model is used.

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

The empirical counter cumulative distribution $ECCDF(n_{m})$ indicates the
probability that in a commodity group, there are at least $n_{m}$ weighted
motifs or anti-motifs and 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 $48\%$ when the directed model benchmark is used and for
$56\%$ 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.

###  IV.3 The NuMeTriS Python package

As an additional result, we release a Python package named ‘NuMeTriS - Null
Models for Triadic Structures’ and containing solvers and routines for triadic
motif analysis for the mentioned models, namely the DBCM, the RBCM and the
mixture models DBCM+CReMa and RBCM+CRWCM. The package is available at the
following URL: https://github.com/MarsMDK/NuMeTriS.

##  V Conclusions

The study of triadic motifs on production networks is still in its infancy due
to a lack of reliable data. In the literature, only one production network has
been characterized, the Japanese one, for a single representative commodity
[36]. In that study a simple message emerged, open triadic subgraphs are over-
represented while closed triadic subgraphs are under-represented. This
phenomenon was explained as due to complementarity, i.e. economic actors
connect in tetradic structures - better explained by open triads - because of
complementary needs [33].

Our production network consists of industries and not firms and our results
are bound to be different. However, even at the level of industries, we find
that an analysis based on a single representative commodity is not enough to
characterize a production network. Product level data is essential to
disaggregate the network into layers that are characterized by commodity-
specific binary motifs and anti-motifs. Moreover, we detected that the
majority of layers are characterized by maximally random triadic structures
when the reciprocal structure is taken into account.

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 [54]. 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.
Moreover, 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 firm-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.

## Acknowledgements

This work is supported by the European Union - NextGenerationEU - National
Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR),
project ‘SoBigData.it - Strengthening the Italian RI for Social Mining and Big
Data Analytics’ - Grant IR0000013 (n. 3264, 28/12/2021). This work has been
also supported by the project ‘Network analysis of economic and financial
resilience’, Italian DM n. 289, 25-03-2021 (PRO3 Scuole) CUP D67G22000130001
and by the PNRR-M4C2-Investimento 1.3, Partenariato Esteso PE00000013-“FAIR-
Future Artificial Intelligence Research”-Spoke 1 “Human-centered AI”, funded
by the European Commission under the NextGeneration EU programme. DG
acknowledges support from the Dutch Econophysics Foundation (Stichting
Econophysics, Leiden, the Netherlands) and the Netherlands Organization for
Scientific Research (NWO/OCW). MDV and DG acknowledge support from the
‘Programma di Attività Integrata’ (PAI) project ‘Prosociality, Cognition and
Peer Effects’ (Pro.Co.P.E.), funded by IMT School for Advanced Studies Lucca.

## Appendix A - Binary Null Models

### The Directed Binary Configuration Model

The Directed Binary Configuration Model (DBCM) is the maximum-entropy model
where the out-degree and in-degree sequences are constrained. The
corresponding Graph Hamiltonian is

$H(A)=\sum_{i}\left(\alpha_{i}^{out}k_{i}^{out}+\alpha_{i}^{in}k_{i}^{in}\right)=\sum_{i,j\neq
i}\left(\alpha_{i}^{out}+\alpha_{j}^{in}\right)a_{ij}.$ (40)

The partition function reads

$\displaystyle Z(A)$
$\displaystyle=\sum_{A}e^{-H(A)}=\sum_{A}e^{-\sum_{i,j\neq
i}\left(\alpha_{i}^{out}+\alpha_{j}^{in}\right)a_{ij}}=$
$\displaystyle=\prod_{i,j\neq
i}\sum_{a_{ij}=0,1}\left(x_{i}^{out}x_{j}^{in}\right)^{a_{ij}}=\prod_{i,j\neq
i}\left(1+x_{i}^{out}x_{j}^{in}\right)$ (41)

where $x_{i}^{y}=e^{-\alpha_{i}^{y}}$ with $\cdot=\\{out,in\\}$.

After computing the partition function the binary graph distribution $P(A)$ is

$\displaystyle P(A)$ $\displaystyle=\dfrac{e^{-H(A)}}{Z(A)}=$
$\displaystyle=\prod_{i,j\neq
i}\dfrac{\left(x_{i}^{out}x_{j}^{in}\right)^{a_{ij}}}{1+x_{i}^{out}x_{j}^{in}}$
(42)

It is possible to define a Log-Likelihood from $P(A)$ as

$\displaystyle\mathcal{L}$ $\displaystyle=\ln P(A)=$
$\displaystyle=-H(A)-\ln(Z(A))=$
$\displaystyle=-\sum_{i}\left(\alpha_{i}^{out}k_{i}^{out}+\alpha_{i}^{in}k_{i}^{in}\right)-\sum_{i,j\neq
i}\ln(1+x_{i}^{out}x_{j}^{in}).$ (43)

Parameters are then estimated using MLE on the log-likelihood function which
consists in solving the following set of equations for the node-specific
parameters $\alpha_{i}^{out},\alpha_{i}^{in}$.

$\begin{cases}\dfrac{\partial\mathcal{L}}{\partial\alpha_{i}^{out}}&=-k_{i}^{out}+\sum_{j\neq
i}\left(\dfrac{x_{i}^{out}x_{j}^{in}}{1+x_{i}^{out}x_{j}^{in}}\right);\\\
\dfrac{\partial\mathcal{L}}{\partial\alpha_{i}^{in}}&=-k_{i}^{in}+\sum_{j\neq
i}\left(\dfrac{x_{i}^{in}x_{j}^{out}}{1+x_{i}^{in}x_{j}^{out}}\right).\end{cases}$
(44)

### The Reciprocal Binary Configuration Model

The Reciprocal Binary Configuration Model (RBCM) is the maximum-entropy model
constraining the reciprocated and non-reciprocated degree sequences. The
corresponding Graph Hamiltonian reads

$\displaystyle H(A)$
$\displaystyle=\sum_{i}\left(\alpha_{i}^{\rightarrow}k_{i}^{\rightarrow}+\alpha_{i}^{\leftarrow}k_{i}^{\leftarrow}+\alpha_{i}^{\leftrightarrow}k_{i}^{\leftrightarrow}\right)=$
$\displaystyle=\sum_{i,j<i}\left(\alpha_{i}^{\rightarrow}+\alpha_{j}^{\leftarrow}\right)a_{ij}^{\rightarrow}+\left(\alpha_{i}^{\leftarrow}+\alpha_{j}^{\rightarrow}\right)a_{ij}^{\leftarrow}+$
$\displaystyle+\left(\alpha_{i}^{\leftrightarrow}+\alpha_{j}^{\leftrightarrow}\right)a_{ij}^{\leftrightarrow}$
(45)

The model-induced partition function is

$\displaystyle Z(A)$ $\displaystyle=\sum_{A}e^{-H(A)}=$
$\displaystyle=\sum_{A}\prod_{i,j<i}\left(x_{i}^{\rightarrow}x_{j}^{\leftarrow}\right)^{a_{ij}^{\rightarrow}}\left(x_{i}^{\leftarrow}x_{j}^{\rightarrow}\right)^{a_{ij}^{\leftarrow}}\left(x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}\right)^{a_{ij}^{\leftrightarrow}}=$
$\displaystyle=\prod_{i,j<i}\sum_{\\{a_{ij}\\}}\left(x_{i}^{\rightarrow}x_{j}^{\leftarrow}\right)^{a_{ij}^{\rightarrow}}\left(x_{i}^{\leftarrow}x_{j}^{\rightarrow}\right)^{a_{ij}^{\leftarrow}}\left(x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}\right)^{a_{ij}^{\leftrightarrow}}=$
$\displaystyle=\prod_{i,j<i}\left(1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}\right)$
(46)

where $x_{i}^{y}=e^{-\alpha_{i}^{y}}$ with
$y=\\{\rightarrow,\leftarrow,\leftrightarrow\\}$, and the last equality arises
because by definition the events described by different arrows - e.g.
$a_{ij}^{\rightarrow}$, $a_{ij}^{\leftarrow}$, $a_{ij}^{\leftrightarrow}$ and
$a_{ij}^{\not{\leftrightarrow}}$ \- are mutually exclusive, i.e. only one of
those terms is equal to $1$ while all the others are equal to zero for each
dyad (i,j).

After computing the partition function the binary Graph Probability $P(A)$ is
obtained as follows

$\displaystyle P(A)$ $\displaystyle=\dfrac{e^{-H(A)}}{Z(A)}=$
$\displaystyle=\prod_{i,j<i}\dfrac{\left(x_{i}^{\rightarrow}x_{j}^{\leftarrow}\right)^{a_{ij}^{\rightarrow}}\left(x_{i}^{\leftarrow}x_{j}^{\rightarrow}\right)^{a_{ij}^{\leftarrow}}\left(x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}\right)^{a_{ij}^{\leftrightarrow}}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}$
(47)

From the partition function, we define the log-likelihood function
$\mathcal{L}$

$\displaystyle\mathcal{L}$ $\displaystyle=\ln P(A)=$
$\displaystyle=-H(A)-\ln(Z(A))=$
$\displaystyle=-\sum_{i}\left(\alpha_{i}^{\rightarrow}k_{i}^{\rightarrow}+\alpha_{i}^{\leftarrow}k_{i}^{\leftarrow}\alpha_{i}^{\leftrightarrow}k_{i}^{\leftrightarrow}\right)-$
$\displaystyle-\sum_{i,j<i}\ln(1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow})$
(48)

and we solve for the node-specific set of parameters
$\\{\alpha_{i}^{\rightarrow},\alpha_{i}^{\leftarrow},\alpha_{i}^{\leftrightarrow}\\}$
by using the MLE framework

$\begin{cases}\dfrac{\partial\mathcal{L}}{\partial\alpha_{i}^{\rightarrow}}&=-k_{i}^{\rightarrow}+\sum_{j\neq
i}\left(\dfrac{x_{i}^{\rightarrow}x_{j}^{\leftarrow}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}\right);\\\
\dfrac{\partial\mathcal{L}}{\partial\alpha_{i}^{\leftarrow}}&=-k_{i}^{\leftarrow}+\sum_{j\neq
i}\left(\dfrac{x_{i}^{\leftarrow}x_{j}^{\rightarrow}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}\right);\\\
\dfrac{\partial\mathcal{L}}{\partial\alpha_{i}^{\leftrightarrow}}&=-k_{i}^{\leftrightarrow}+\sum_{j\neq
i}\left(\dfrac{x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}{1+x_{i}^{\rightarrow}x_{j}^{\leftarrow}+x_{i}^{\leftarrow}x_{j}^{\rightarrow}+x_{i}^{\leftrightarrow}x_{j}^{\leftrightarrow}}\right).\end{cases}$
(49)

## Appendix B - Conditional Weighted Null Models

### Conditional Reconstruction Method A

The Conditional Reconstruction Method A (CReMa) is the conditional maximum-
entropy model constraining the out-strength and in-strength sequences.

The Graph Hamiltonian reads

$\displaystyle H(W)$
$\displaystyle=\sum_{i}\left(\beta_{i}^{out}s_{i}^{out}+\beta_{i}^{in}s_{i}^{in}\right)$
$\displaystyle=\sum_{i,j\neq
i}\left(\beta_{i}^{out}+\beta_{j}^{in}\right)w_{ij}.$ (50)

It induces a partition function of the following form

$Z(W_{A})=\int_{W_{A}}\prod_{i,j\neq
i}e^{-\left(\beta_{i}^{out}+\beta_{j}^{in}\right)w_{ij}}dw_{ij}=\left(\dfrac{1}{\beta_{i}^{out}+\beta_{j}^{in}}\right)^{a_{ij}}.$
(51)

The conditional graph probability function $Q(W|A)$ is then defined as

$Q(W|A)=\prod_{i,j\neq
i}\left[(\beta_{i}^{out}+\beta_{j}^{in})e^{-(\beta_{i}^{out}+\beta_{j}^{in})w_{ij}}\right]^{a_{ij}}$
(52)

and the corresponding log-likelihood $\mathcal{L}$ is

$\displaystyle\mathcal{L}$
$\displaystyle=-\sum_{i}\left(\beta_{i}^{out}s_{i}^{out}+\beta_{i}^{in}s_{i}^{in}\right)+\sum_{i,j\neq
i}a_{ij}\ln(\beta_{i}^{out}+\beta_{j}^{in})$ (53)

In order to take into account the random variability of the binary adjacency
matrix $A$ we average the log-likelihood over the ensemble realizations and
obtain the generalized log-likelihood $\mathcal{G}$, defined as

$\displaystyle\mathcal{G}$
$\displaystyle=-\sum_{i}\left(\beta_{i}^{out}s_{i}^{out}+\beta_{i}^{in}s_{i}^{in}\right)+\sum_{i,j\neq
i}p_{ij}\ln(\beta_{i}^{out}+\beta_{j}^{in})$ (54)

The node-specific parameters $\\{\beta_{i}^{out},\beta_{i}^{in}\\}$ are tuned
according to the first order equation for GLE, i.e.

$\begin{cases}\dfrac{\partial\mathcal{G}}{\partial\beta_{i}^{out}}&=-s_{i}^{out}+\sum_{j\neq
i}\dfrac{p_{ij}}{\beta_{i}^{out}+\beta_{j}^{in}}\\\
\dfrac{\partial\mathcal{G}}{\partial\beta_{i}^{in}}&=-s_{i}^{in}+\sum_{j\neq
i}\dfrac{p_{ji}}{\beta_{i}^{in}+\beta_{j}^{out}}\end{cases}$ (55)

which amounts to equating the empirical and model-induced out-strength and in-
strength sequences.

### Conditionally Reciprocal Weighted Configuration Model

The Conditionally Reciprocal Weighted Configuration Model (or CRWCM) is a
novel conditional maximum-entropy model that constrains out-strength and in-
strength dividing them according to the character of the reciprocity in the
underlying links. The constraints are non-reciprocated out-strengths and in-
strength - $s_{i}^{\rightarrow}$ and $s_{i}^{\leftarrow}$ \- and reciprocated
out-strengths and in-strengths, namely $s_{i}^{\leftrightarrow,out}$ and
$s_{i}^{\leftrightarrow,in}$.

The corresponding Graph Hamiltonian reads

$\displaystyle H(W)$
$\displaystyle=\sum_{i}\left(\beta_{i}^{\rightarrow}s_{i}^{\rightarrow}+\beta_{i}^{\leftarrow}s_{i}^{\leftarrow}\right)+$
$\displaystyle+\left(\beta_{i}^{\leftrightarrow,out}s_{i}^{\leftrightarrow,out}+\beta_{i}^{\leftrightarrow,in}s_{i}^{\leftrightarrow,in}\right)=$
$\displaystyle=\sum_{i,j\neq i}h(w_{ij})$ (56)

where

$h(w_{ij})=\left(\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow}\right)a_{ij}^{\rightarrow}w_{ij}+\left(\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in}\right)a_{ij}^{\leftrightarrow}w_{ij}.$
(57)

The Hamiltonian induces a conditional partition function defined as

$\displaystyle Z(W_{A})$ $\displaystyle=\int_{W_{A}}e^{-\sum_{i,j\neq
i}h(w_{ij})}dw_{ij}=$ $\displaystyle=\prod_{i,j\neq
i}\int_{0}^{\infty}e^{-h(w_{ij})}dw_{ij}=\prod_{i,j\neq i}Z_{ij|A}$ (58)

where the dyadic-specific conditional partition function $Z_{ij|A}$ is

$Z_{ij|A}=\left(\dfrac{1}{\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow}}\right)^{a_{ij}^{\rightarrow}}\left(\dfrac{1}{\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in}}\right)^{a_{ij}^{\leftrightarrow}}$
(59)

so that $Z_{ij|A}=(\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow})^{-1}$ if
$a_{ij}^{\rightarrow}=1$ and
$Z_{ij|A}=(\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in})^{-1}$
if $a_{ij}^{\leftrightarrow}=1$.

The conditional probability distribution for the weighted network is

$Q(W|A)=\prod_{i,j\neq i}\dfrac{e^{-h(w_{ij})}}{Z_{ij|A}}.$ (60)

It induces a log-likelihood $\mathcal{L}$ of the form

$\mathcal{L}=\sum_{i,j\neq i}\left[-h(w_{ij})+\ln(Z_{ij|A})\right]$ (61)

We take into account the random variability of $A$, induced by the choice of
the binary model, by formulating the average log-likelihood, i.e. the
generalized log-likelihood $\mathcal{G}$. It is defined as

$\mathcal{G}=\sum_{i,j\neq
i}\left[-h(w_{ij})+\ln(Z_{ij|P(A)})\right]=\mathcal{G}^{\rightarrow}+\mathcal{G}^{\leftrightarrow}$
(62)

where

$\begin{cases}\mathcal{G}^{\rightarrow}&=\sum_{i,j\neq
i}\left[-(\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow})w_{ij}+p_{ij}^{\rightarrow}\ln(\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow})\right]\\\
\mathcal{G}^{\leftrightarrow}&=\sum_{i,j\neq
i}-(\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in})w_{ij}+\\\
&+p_{ij}^{\leftrightarrow}\ln(\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in}),\end{cases}$
(63)

i.e. $\mathcal{G}$ can be decoupled into a non-reciprocated component
$\mathcal{G}^{\rightarrow}$ and a reciprocated component
$\mathcal{G}^{\leftrightarrow}$. This implies that the resulting GLE problem
can be divided into two subproblems. The non-reciprocated problem equates to
solving a system of $2N$ coupled equations, namely

$\begin{cases}\dfrac{\partial\mathcal{G}}{\partial\beta_{i}^{\rightarrow}}&=-s_{i}^{\rightarrow}+\sum_{j\neq
i}\dfrac{p_{ij}^{\rightarrow}}{\beta_{i}^{\rightarrow}+\beta_{j}^{\leftarrow}}\\\
\dfrac{\partial\mathcal{G}}{\partial\beta_{i}^{\leftarrow}}&=-s_{i}^{\leftarrow}+\sum_{j\neq
i}\dfrac{p_{ij}^{\leftarrow}}{\beta_{i}^{\leftarrow}+\beta_{j}^{\rightarrow}}.\end{cases}$
(64)

The reciprocated subproblem equates to solving the following set of $2N$
coupled equations

$\begin{cases}\dfrac{\partial\mathcal{G}}{\partial\beta_{i}^{\leftrightarrow,out}}&=-s_{i}^{\leftrightarrow,out}+\sum_{j\neq
i}\dfrac{p_{ij}^{\leftrightarrow}}{\beta_{i}^{\leftrightarrow,out}+\beta_{j}^{\leftrightarrow,in}}\\\
\dfrac{\partial\mathcal{G}}{\partial\beta_{i}^{\leftrightarrow,in}}&=-s_{i}^{\leftrightarrow,in}+\sum_{j\neq
i}\dfrac{p_{ij}^{\leftrightarrow}}{\beta_{i}^{\leftrightarrow,in}+\beta_{j}^{\leftrightarrow,out}}.\end{cases}$
(65)

## References

  * [1] D. Acemoglu, V. M. Carvalho, A. Ozdaglar, and A.Tahbaz-Salehi, _Econometrica_ 80, 1977 (2012).
  * [2] D. Aobdia, J. Caskey, and N. B. Ozel, _Review of Accounting Studies_ 19, 1191 (2014).
  * [3] E. Atalay, _American Economic Journal: Macroeconomics_ 9, 254 (2017).
  * [4] H. Bouakez, E. Cardia, and F. J. Ruge-Murcia, _International Economic Review_ 50, 1243 (2009).
  * [5] A. Brintrup, P. Wichmann, P. Woodall, D. McFarlane, E. Nicks, and W. Krechel, _Complexity_ 2018, e9104387 (2018).
  * [6] A. Pichler and J. D. Farmer, _Economic Systems Research_ 34, 273 (2022).
  * [7] A. Bacilieri, A. Borsos, P. Astudillo-Estévez and F. Lafond, INET Oxford Working Paper No. 2023-08. (2023)
  * [8] E. Atalay, A. Hortaçsu, J. Roberts, and C. Syverson, _Proceedings of the National Academy of Sciences_ 108, 5199 (2011).
  * [9] A. B. Bernard, A. Moxnes, and Y. U. Saito, _Journal of Political Economy_ 127, 639 (2019).
  * [10] G. Buiten, E. Jonge, G. Mooijen, S. Hooijmaaijers, and P. Bogaart, OECD Conference, New Analytical Tools and Techniques for Economic Policymaking (2021).
  * [11] V. M. Carvalho, M. Nirei, Y. U. Saito, and A. Tahbaz- Salehi, _The Quarterly Journal of Economics_ 136, 1255 (2021).
  * [12] V. M. Carvalho and A. Tahbaz-Salehi, _Annual Review of Economics_ 11, 635 (2019).
  * [13] L. Cohen and A. Frazzini, _The Journal of Finance_ 63 (2008).
  * [14] L. Mungo, F. Lafond, P. Astudillo-Estévez and J. D. Farmer, _Journal of Economic Dynamics and Control_ 148, 104607 (2023).
  * [15] E. Dhyne, G. Magerman, and S. Rubìnova, The Belgian production network 2002-2012 , Working Paper 288 (NBB Working Paper, 2015).
  * [16] E. Dhyne, A. K. Kikkawa, M. Mogstad, and F. Tintelnot, _The Review of Economic Studies_ 88, 643 (2021).
  * [17] C. Diem, A. Borsos, T. Reisch, J. Kertész, and S. Thurner, _Scientific Reports_ 12, 7719 (2022).
  * [18] M. Cardoza, F. Grigoli, N. Pierri, and C. Ruane, IMF Working paper, No. 20/205 (2020).
  * [19] P. W. Chacha, B. Kirui, and V. Wiedemann, Mapping Kenya’s Production Network Transaction by Transaction. Oxford WP
  * [20] B. Demir, B. Javorcik, T. K. Michalski, and E. Ors, _The Review of Economics and Statistics_ , pages 1–46, (2022).
  * [21] J. L. Peydr´o, G. Jim´enez, H. Kenan, E. Moral-Benito, and F. Vega-Redondo, CEPR Discussion Paper (2020)..
  * [22] R. Newfarmer, J. Page, and F. Tarp, Industries without Smokestacks: Industrialization in Africa Reconsidered (Oxford, 2018)
  * [23] A. Kumar, A. S. Chakrabarti, A. Chakraborti, T. Nandi, _Physica A: Statistical Mechanics and its Applications_ 568, 125714 (2021).
  * [24] H. Goto, H. Takayasu, and M. Takayasu, _PLoS ONE_ 12, e0185712 (2017).
  * [25] S. Hooijmaaijers and G. Buiten, OECD Conference, New Analytical Tools and Techniques for Economic Policy-making (2019).
  * [26] L. N. Ialongo, C. de Valk, E. Marchese, F. Jansen, H. Zmarrou, T. Squartini, and D. Garlaschelli, _Scientific Reports_ 12, 11847 (2022).
  * [27] H. Inoue and Y. Todo, _Nature Sustainability_ 2, 841 (2019).
  * [28] H. Inoue and Y. Todo, _PLoS ONE_ 15, e0239251 (2020).
  * [29] Y. Kashiwagi, Y. Todo, and P. Matous, _Review of International Economics_ 29, 1186 (2021).
  * [30] M. D. König, A. Levchenko, T. Rogers, and F. Zilibotti, _Proceedings of the National Academy of Sciences_ 119, e2203730119 (2022).
  * [31] E. E. Kosasih and A. Brintrup, _International Journal of Production Research_ 60, 5380 (2022).
  * [32] J. Maluck, R. V. Donner, H. Takayasu, and M. Takayasu, _Journal of Statistical Mechanics: Theory and Experiment_ 2017, 053404 (2017).
  * [33] C. E. S. Mattsson, F. W. Takes, E. M. Heemskerk, C. Diks, G. Buiten, A. Faber, and P. M. A. Sloot, _Frontiers in Big Data_ 4 (2021).
  * [34] J. McNerney, C. Savoie, F. Caravelli, V. M. Carvalho, and J. D. Farmer, _Proceedings of the National Academy of Sciences_ 119, e2106031118 (2022).
  * [35] T. Mizuno, W. Souma, and T. Watanabe, _PLoS ONE_ 9, e100712 (2014).
  * [36] T. Ohnishi, H. Takayasu, and M. Takayasu, _Journal of Economic Interaction and Coordination_ 5, 171 (2010).
  * [37] A. Rachkov, F. Pijpers, and D. Garlaschelli, CBS Technical Reports 10.13140/RG.2.2.31861.29925 (2021).
  * [38] M. Taschereau-Dumouchel, 2017 Meeting Papers, 700, Society for Economic Dynamics (2017).
  * [39] H. Watanabe, H. Takayasu, and M. Takayasu, _Physica A: Statistical Mechanics and its Applications_ 392, 741 (2013).
  * [40] C. Diem, A. Borsos, T. Reisch, J. Kertész and S. Thurner, _arXiv:2302.11451_ (2023).
  * [41] J. Maluck and R. V. Donner, _PLoS ONE_ 10, e0133310 (2015).
  * [42] Z. Wang, S. Liu, C. Han, S. Huang, X. Gao, R. Tang, and Z. Di, _Frontiers in Physics_ 9 (2022).
  * [43] A. Alfaro-Ureña, M. Fuentes, I. Manelici, and J. Vasquez, Research Paper Series, Banco Central De Costa Rica (2018).
  * [44] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, _Science_ 298, 824 (2002).
  * [45] S. S. Shen-Orr, R. Milo, S. Mangan, and U. Alon, _Nature Genetics_ 31, 64 (2002).
  * [46] A. Stivala and A. Lomi, _Applied Network Science_ 6, 1 (2021).
  * [47] A. Asikainen, G. Iñiguez, J. Ureña-Carrión, K. Kaski, and M. Kivelä, _Science Advances_ 6, eaax7310 (2020).
  * [48] T. Squartini and D. Garlaschelli, Self-Organizing Systems 7166, 24 (2012).
  * [49] A. Maratea, A. Petrosino, and M. Manzo, _Procedia Computer Science_ 98, 479 (2016).
  * [50] T. Squartini, I. van Lelyveld, and D. Garlaschelli, _Scientific Reports_ 3, 3357 (2013).
  * [51] T. Squartini and D. Garlaschelli, _Journal of Complex Networks_ 3, 1 (2015).
  * [52] P. Colomer-de Simón, M. Serrano, M. G. Beiró, J. I. Alvarez-Hamelin, and M. Boguñá, _Scientific Reports_ 3, 2517 (2013).
  * [53] A. Jamakovic, P. Mahadevan, A. Vahdat, M. Boguñá, and D. Krioukov, _arXiv.0908.1143_ (2009).
  * [54] F. Picciolo, F. Ruzzenenti, P. Holme, and R. Mastrandrea, _New Journal of Physics_ 24, 053056 (2022).
  * [55] E. T. Jaynes, _Physical Review_ 106, 620 (1957).
  * [56] E. T. Jaynes, _Physical Review_ 108, 171 (1957).
  * [57] E. T. Jaynes, _Proceedings of the IEEE_ 70, 939 (1982).
  * [58] D. Garlaschelli and M. I. Loffredo, _Physical Review E_ 78, 015101 (2008).
  * [59] M. Bardoscia, P. Barucca, S. Battiston, F. Caccioli, G. Cimini, D. Garlaschelli, F. Saracco, T. Squartini, and G. Caldarelli, _Nature Reviews Physics_ 3, 490–507 (2021).
  * [60] G. Cimini, T. Squartini, F. Saracco, D. Garlaschelli, A. Gabrielli, and G. Caldarelli, _Nature Reviews Physics_ 1, 58 (2019).
  * [61] G. Cimini, R. Mastrandrea, and T. Squartini, Reconstructing Networks, Elements in the Structure and Dynamics of Complex Networks (Cambridge University Press, 2021).
  * [62] T. Squartini and D. Garlaschelli, Maximum-Entropy Networks: Pattern Detection, Network Reconstruction and Graph Combinatorics, SpringerBriefs in Complexity (Springer Cham, 2017).
  * [63] D. Garlaschelli and M. I. Loffredo, _Physical Review Letters_ 93, 188701 (2004).
  * [64] T. Squartini, G. Fagiolo, and D. Garlaschelli, _Physical Review E_ 84, 046117 (2011).
  * [65] T. Squartini, G. Fagiolo, and D. Garlaschelli, _Physical Review E_ 84, 046118 (2011).
  * [66] R. Mastrandrea, T. Squartini, G. Fagiolo, and D. Garlaschelli, _New Journal of Physics_ 16, 043022 (2014).
  * [67] F. Parisi, T. Squartini, and D. Garlaschelli, emphNew Journal of Physics 22, 053053 (2020).
  * [68] A. Almog, R. Bird, and D. Garlaschelli, _Frontiers in Physics_ 7, 10.3389/fphy.2019.00055 (2019).
  * [69] M. Di Vece, D. Garlaschelli, and T. Squartini, _Physical Review Research_ 4, 033105 (2022).
  * [70] M. Di Vece, D. Garlaschelli, and T. Squartini, _Chaos, Solitons & Fractals_ 166, 112958 (2023).
  * [71] G. Cimini, T. Squartini, D. Garlaschelli, and A. Gabrielli, _Scientific Reports_ 5, 15758 (2015).
  * [72] T. Squartini and D. Garlaschelli, _New Journal of Physics_ 13, 083001 (2011).
  * [73] K. Anand, I. van Lelyveld, Banai, S. Friedrich, R. Garratt, G. Hałaj, J. Fique, I. Hansen, S. M. Jaramillo, H. Lee, J. L. Molina-Borboa, S. Nobili, S. Rajan, D. Salakhova, T. C. Silva, L. Silvestri, and S. R. S. d. Souza, _Journal of Financial Stability_ 35, 107 (2018).
  * [74] M. Lebacher, S. Cook, N. Klein, and G. Kauermann, _Journal of Network Theory in Finance_ 5, 29 (2019).
  * [75] A. Ramadiah, F. Caccioli, and D. Fricke, _Journal of Economic Dynamics and Control_ 111, 103817 (2020).
  * [76] P. Mazzarisi and F. Lillo, in Econophysics and Sociophysics: Recent Progress and Future Directions, New Economic Windows, (Springer International Publishing, Cham, 2017) pp. 201–215.
  * [77] T. Squartini, F. Picciolo, and F. Ruzzenenti, _Scientific Reports_ 3, 2729 (2013).
  * [78] S. S. Shapiro and M. B. Wilk, _Biometrika_ 52, 591 (1965).