Supply Chain Network Reconstruction

Network reconstruction methods are designed to estimate a complete network based on a partial observation of the underlying network. This can include individual nodes or connections, as well as global observations. For example, the number of connections for each node or simply the distribution of connections per node. We refer here to a directed network where the nodes are individual businesses, and the network connections indicate the supply or receipt of goods and services between two businesses.

Given the missing information, it is usually not possible to perfectly reconstruct the network. Instead, the aim is to generate an ensemble \(\mathbb{G}\) of all possible networks that satisfy the observed data. Specifically, to generate a probability matrix \(\small{P=[{p_{ij}]}_{i,j=1}^N}\) that encodes the distribution of all possible networks. Here, each \(p_{ij}\) describes the probability \(\small{0<p_{ij}<1}\) of a connection between \(i\) and \(j\). A Bernoulli trial of each \(p_{ij}\) results in a graph \(\small G\in\mathbb{G} \) described by an adjacency matrix \(\small{A=[{a_{ij}]}_{i,j=1}^N}\) where each element \(a_{ij}=1\) if a link exists, otherwise \(a_{ij}=0\). After reconstructing a graph \(\small G\in\mathbb{G} \), we then need to ascertain the weight \(w_{ij}\) or value of goods and services being traded between businesses. Therefore, we construct a weighted adjacency matrix \(\small{W=[{w_{ij}]}_{i,j=1}^N}\) to describe the weight and direction of goods and services from node \(i\) to node \(j\).

For the reconstruction to be statistically defensible, it should minimise bias with respect to all available information. Maximising entropy ensures that the link probabilities are 'maximally random' and therefore captures the uncertainty inherent in the estimate. Moreover, this approach reproduces the values of the information constraints on average, ensuring the outcome is unbiased. In this regard, entropy maximisation can be considered as a more general (non-parametric) form of Bayes’ theorem.

To maximise entropy, first consider that the reconstructed network contains information \(\small{I\left(G\right)=-\text{ln}P\left(G\right)}\) about the underlying network (Cover & Thomas, 2006; Squartini et al., 2018). Where a reconstructed network \(\small{G}\) is unlikely \(\small{P\left(G\right)\simeq0}\), the information presents an almost infinite amount of 'surprise' (Cover & Thomas, 2006; Squartini et al., 2018). Where the reconstructed network is certain \(\small{P\left(G\right)=1}\), no new information about the underlying network is presented. Iterating over the ensemble of all possible networks results in the Shannon entropy equation:

\(S=\left\langle I\right\rangle=\sum_{G\in\mathbb{G}} P\left(G\right)I\left(G\right)=\sum_{G\in\mathbb{G}}{-P\left(G\right)\text{ln} P\left(G\right)}\)

(Squartini et al., 2018)

From here, the aim is to maximise entropy subject to \(m\) constraints \(C_m\), each representing known information of the network. In general the constraints take the form:

\(\sum_{G\in\mathbb{G}}{P\left(G\right)C_m(G)=\left\langle C_m\right\rangle}\)

where the first constraint acts as a normalising condition,  \( \sum_{G\in\mathbb{G}}{P\left(G\right)=1} \).

Here \(C_m(G)\) is the value of the \(m^{th}\) constraint for graph \(G\) and \(\small \langle C_m \rangle\) is the observed (or expected) value for constraint \(m\). Constraining entropy by the expected value ensures the subsequent distribution preserves our observation in expectation (on average). Therefore, the Lagrangian maximising entropy with respect to our \(m\) constraints follows:

\(\displaystyle \mathcal{L}\left[P\right]=S-\lambda_0\left[\sum_{G\in\mathbb{G}}{P\left(G\right)-1}\right]-\sum_{m=1}^{M}{\lambda_m\left[\sum_{G\in\mathbb{G}}{P\left(G\right)C_m(G)-\left\langle C_m\right\rangle}\right]}\)

Solving the partial derivatives results in the exponential random graph model distribution:

\(\sum_{G\in\mathbb{G}}{P\left(G\right)} =\displaystyle \frac{e^{-\sum_{m=1}^M \lambda_m \cdot C_m(G)}} {\sum_{G\in\mathbb{G}} \exp (-\sum_{m=1}^M \lambda_m \cdot C_m(G))} = \frac{e^{-H(G)}}{Z} \)

(Squartini et al., 2018)

In statistical mechanics terms, the Hamiltonian \(\small{H\left(G\right)=\sum_{m=1}^{M}{\lambda_mC_m\left(G\right)}}\) describes the total energy of the system, and the partition function \(\small{Z=\sum_{G\in\mathbb{G}} \exp\left(-\sum_{m=1}^{M}{\lambda_mC_m\left(G\right)}\right)}\) acts as a normalising constant (Squartini et al., 2018).

If the degree or number of connections \(k\) for each node is known, then one of the constraints becomes \(\sum_{G\in\mathbb{G}}\sum_{i}{P\left(G\right)k_i\left(G\right)=\left\langle k\right\rangle}\ \)where \(k_i\left(G\right)\) is the degree of node \(i\) for graph \(G\). Given the degree \(k_i=\sum_{j} a_{ij}\) is the number of connections for each node and the network is directed, each node has an in-degree \(k_i^{in}\) and an out-degree \(k_i^{out}\). So, the Hamiltonian takes the form:

\(H\left(G\right)=\sum_{i}\sum_{j\neq i}{\lambda_i^{out}a_{ij}+\lambda_j^{in}a_{ji}}\ =\sum_{i}\sum_{j\neq i}{\lambda_i^{out}k_i^{out}+\lambda_j^{in}k_j^{in}}\)

(Squartini et al., 2018)

Now that we have defined the Hamiltonian at the node-level, we can also define the probability distribution at the node-level. Consider that the probability \(p_{ij}\) of a connection between \(i\) and \(j\) is equal to the expected value of a connection existing \(\small \left\langle a_{ij} \right\rangle\) (Rachkov, Pjipers, & Garlashceilli, 2021). Therefore, the network can be split into pairs of directed links:

\(\displaystyle{p_{ij}=\frac{e^{-\sum_{i}{\lambda_i^{\text{out}}k_i^{\text{out}}}}e^{-\lambda_j^{\text{in}}k_j^{\text{in}}}}{1+e^{-\sum_{i}{\lambda_i^{\text{out}}k_i^{\text{out}}}}e^{-\lambda_j^{\text{in}}k_j^{\text{in}}}}}\)

For simplicity, we use the convention that \(x_j^{in}=e^{-\lambda_j^{in}k_j^{in}} \) and \(x_i^{out}=e^{-\sum_{i}{\lambda_i^{out}k_i^{out}}}\). Therefore:

\(\displaystyle{p_{ij}=\frac{x_i^{\text{out}}x_j^{\text{in}}}{1+x_i^{\text{out}}x_j^{\text{in}}}}\)

Recall that each \(p_{ij}\) is an element in a probability matrix \(P\) that encodes the probability of selecting a realisation \(G\) of all possible graphs \(\mathbb{G}\):

\(P=\left[\begin{matrix}p_{ij}^\ &\cdots&p_{iN}^\ \\\vdots&\ddots&\vdots\\p_{Nj}^\ &\cdots&p_{NN}^\ \\\end{matrix}\right]\)

The degree, or number of connections per node, is usually unknown given the confidentiality of financial flows. However, Caldarelli et al. (2002) show that the link probability can be derived from a node's 'fitness', as opposed to any topological property. Here, fitness denotes a node's influence in the network, meaning high fitness nodes attract more connections. This method has been demonstrated to hold in international trade networks (Almog, Squartini, & Garlaschelli, 2015; Garlaschelli & Loffredo, 2008), interbank networks (De Masi et al., 2006; Garlaschelli & Loffredo, 2004; Bollobas et al., 2003), and financial securities networks (Cimini et al., 2015; Garlaschelli et al., 2005). Therefore, a 'fitness ansatz' \(\small{z}\) correlated with the Lagrangian multipliers is used to approximate the nodes' degrees. 

The theory of node fitness has also been observed in highly detailed domestic trade networks in Japan (Bernard, Moxnes and Saito, 2015) and Belgium (e.g. Bernard, Dhyne, Magerman, Manova and Moxnes, 2018; Dhyne and Duprez, 2017). These examples show that:

1. Larger firms generally have more trade connections.
2. Larger firms tend to transact with other large firms.
3. Firms tend to transact with geographically closer firms.

Additionally, the presence of the power law in many aspects such as the connection distribution (Bernard et al., 2018; Bernard, Moxnes and Saito, 2015), distance distribution (Bernard, Moxnes and Saito, 2015; Dhyne and Duprez, 2017) and weight distribution (Bernard et al., 2018) affirm the structural properties of the network.

Unlike the degrees, the weight of incoming and outgoing flows (total sales and purchases) for each businesses is generally known (e.g. from publicly disclosed balance sheets or taxation data). This aggregated weight or strength is defined for each node as \(\small{\hat{s}_i^{out}=\sum_{j=1}^{N} w_{ij}}\) for the out-strength, and  \(\small{\hat{s}_j^{in}=\sum_{i=1}^{N} w_{ij}}\) for the in-strength. For businesses, this can be considered as the total value of purchases (out-strength) and sales (in-strength) to other businesses. Therefore the link probability can be defined as:

\(\displaystyle{p_{ij}=\frac{x_i^{\text{out}}x_j^{\text{in}}}{1+x_i^{\text{out}}x_j^{\text{in}}}\approx\frac{zs_i^{\text{out}}s_j^{\text{in}}}{1+zs_i^{\text{out}}s_j^{\text{in}}}}\)

Notice that under this formulation the sum of incoming and outgoing probabilities for each node is equal to the expected in-degree \(\small{\sum_{j} p_{ij}=\langle k_i^{\text{out}} \rangle}\) and out-degree \(\small{\sum_{i} p_{ij}=\langle k_j^{\text{in}} \rangle}\). That is, the row-sum and column-sum of the \(\small{P}\) matrix describe the expected number of links for each node. It follows too, that the probability matrix sums to the expected number of links \(\small{\langle L \rangle}\) in the network:

\(\displaystyle{\sum_{i}\sum_{j\neq i} p_{ij}=\langle L \rangle}\)

This implies that if the total number of links (or even a sample of links) in the network is known, then \(z\) can be calculated and shown to meet the expected value of the Lagrangian multipliers (Squartini et al., 2018). Furthermore, the reconstruction takes on a more intuitive process, wherein the expected number of links for each firm is simply the total number of links in the market, apportioned by each firm’s market share.

Finally, Squartini et al. (2018) demonstrate the validity of incorporating a gravity model, which employs a distance function \(\small{f(d_{ij})}\) to penalise connections between geographically distant nodes, thereby accounting for logistics costs observed in the Japanese and Belgium trade networks (Bernard, Moxnes and Saito, 2015; Dhyne and Duprez, 2017):

\( \displaystyle p_{ij}=\frac{z\hat{s}_i^{out}\hat{s}_j^{in}e^{f(d_{ij})}}{1+z\hat{s}_i^{out}\hat{s}_j^{in}e^{f(d_{ij})}} \)

The Dutch Central Bureau of Statistics (CBS) has demonstrated the potential for applying network reconstruction to an entire economic network (Rachkov, Pijpers, & Garlaschelli, 2021; Hooijmaaijers & Buiten, 2019). Drawing from the work of Squartini et al. (2018), the CBS adapts the maximum entropy framework to the system of national accounts. Specifically, Rachkov, Pijpers, & Garlaschelli (2021) make use of supply-use tables to determine commodity flows and input-output tables to determine industry flows:

\( \displaystyle p_{ij}^{[\alpha]}=\frac{(z^{[\alpha]}\hat{s}_i^{out,[\alpha]}\hat{s}_j^{in,[\alpha]}I_{ij})d_{ij}} {1+(z^{[\alpha]}\hat{s}_i^{out,[\alpha]}\hat{s}_j^{in,[\alpha]}I_{ij})d_{ij}}\)

where:

• the link probabilities of each commodity group \(\alpha\) are estimated independently;
• the out-strength \(\small{\textstyle \hat{s}_j^{out,[\alpha]}}\) represents the volume of goods and services of commodity group \(\textstyle \alpha\) that firm \(\textstyle j\) supplies;
• the in-strength \(\small{\textstyle \hat{s}_i^{in,[\alpha]}}\) represents the volume of goods and services of commodity group \(\alpha\) that firm \(i\) purchases;

•  \(\small c_i\) refers to the industry classification of business \(i\);

• \(\small I_{ij}\) is a binary indicator that takes the value \(\small 1\) where industry \(\textstyle c_i\) trades with industry \(\textstyle c_j\), and \(0\) otherwise; and
• \(\small d_{ij} \) is the sum of the absolute difference in location (i.e., longitude and latitude coordinates) between two businesses.

Note that the in-strength and out-strength are defined for each product market (or commodity group) as:

\(\hat{s}_j^{in,[\alpha]} =\frac{\mathrm{Net \ turnover \ of \ firm \ }j} {\mathrm{Total\ net \ turnover \ of \ industry \ }c} \small{D_{c_j}^{in,[\alpha]}}\)       and          \(\hat{s}_i^{out,[\alpha]} =\frac{\mathrm{Net \ turnover \ of \ firm \ }i} {\mathrm{Total\ net \ turnover \ of \ industry\ }c} \small{D_{c_i}^{out,[\alpha]}}\) 

where:

•  \(\small D_{c_j}^{in,[\alpha]} \) represents the value of commodity group \(\alpha\) used in industry \(c\); and
• \(\small D_{c_i}^{out,[\alpha]} \)  represents the value of commodity group \(\alpha\) supplied by industry \(c\).

The simplest method for generating weights follows the maximum entropy approach. Maximising the Shannon entropy \(\small S=\sum_{i=1}^{N}\sum_{j=1}^{N} w_{ij}\text{ ln }w_{ij}\) leads to \(\small w_{ij}^{ME}=\frac{s_i^{out}s_j^{in}}{W}\ \)\(\forall i,j\) where \(\small W=\sum_{j=1}^{N}s_j^{in}\equiv\sum_{i=1}^{N}s_i^{out}\) is the total network weight. Here, the weights are conditional on the existence of a link. That is, each weight matrix is conditional on the graph (adjacency matrix) drawn from the probability matrix \(\small P\). The marginal weights and marginal degrees are met on average (Squartini et al., 2018).

An alternative approach described by Parisi, Squartini, & Garlaschelli (2020) takes the probability matrix \(\small P\) as prior information and maximises the Shannon entropy, constrained by the observed strength. The conditional entropy maximisation produces a distribution of weighted connections described by the tensor of coefficients \(\beta_{ij}\) that can be sampled:

\(\displaystyle\beta_{ij}=\frac{p_{ij}}{w_{ij}^{ME}}\ \ \ \ \ \forall i\neq j\)

(Parisi, Squartini, & Garlaschelli, 2020)

Despite the size of the tensor, the method is computationally feasible given that each parameter can be calculated independently, enabling parallelisation (Parisi, Squartini, & Garlaschelli, 2020).

Observations of the degree (or degree distribution) are required to calibrate and validate the reconstruction. Although bank transaction data would be ideal, practical limitations mean that this source would require significant investment to ensure it is usable for statistical purposes. Similarly, consignment data from freight transportation or package tracking information may prove to be viable, but currently faces similar practical limitations (aside from the obvious absence of digital goods and services). Instead, the ABS is considering manageable and scalable alternatives.

Administrative Data 

The introduction of single-touch payroll reporting has proliferated the use of common accounting software for businesses with more than 19 employees (Australian Tax Office, 2020). As a result, there has been a concurrent uptake in e-invoicing (i.e., using accounting software to send digital invoices). As such, a small number of accounting software companies provide business accounting software services (including e-invoices) to the majority of businesses. Although the coverage is skewed towards small businesses, the data covers pairs of trading businesses (buyers and suppliers). Therefore, large businesses will be present as suppliers to small businesses.

E-invoicing is a key feature of the Government’s modern digital economy strategy (Prime Minister of Australia, 2021; Department of Treasury, 2020). Although on hold due to the COVID-19 pandemic, it is reasonable to expect that e-invoicing will eventually be rolled out nation-wide (Prime Minister & Cabinet, 2021). Therefore, in the long term, e-invoicing data with near complete coverage could be suitable for statistical purposes.

Survey Data

Existing surveys could be amended to include information on trading relationships. To preserve privacy and reduce administrative burden, these questions would only need to capture aggregate relational data (Breza et al., 2020). For example, additional questions could include: 'How many businesses in industry \(c\) do you supply commodity \(\alpha\) to?', and 'What is the total value?' Therefore, in principle, existing surveys can be tailored to inform reconstruction without identifying individual businesses.