# Revealing production networks from firm growth dynamics

Luca Mungo  Mathematical Institute and Institute for New Economic Thinking at
the Oxford Martin School, University of Oxford, Oxford, United Kingdom
<EMAIL_ADDRESS>José Moran  Mathematical Institute and Institute for
New Economic Thinking at the Oxford Martin School, University of Oxford,
Oxford, United Kingdom  Complexity Science Hub Vienna, Josefstädter Straße 39,
A-1080, Austria

###### Abstract

We study the correlation structure of firm growth rates. We show that most
firms are correlated because of their exposure to a common factor but that
firms linked through the supply chain exhibit a stronger correlation on
average than firms that are not. Removing this common factor significantly
reduces the average correlation between two firms with no relationship in the
supply chain while maintaining a significant correlation between two firms
that are linked. We then demonstrate how this observation can be used to
reconstruct the topology of a supply chain network using Gaussian Markov
Models.

## Introduction

Fifty years ago, Wassily Leontief was awarded the Nobel prize in Economics for
his development of the input-output method and its application to important
economic problems.111 See e.g. https://www.nobelprize.org/prizes/economic-
sciences/1973/summary/. Interestingly, Leontief took inspiration from the
tableau économique Quesnay (1766) of the physician-turned-economist Quesnay, a
member of the physiocratic school of economic thought, which saw the economy
much like a human body. To Quesnay, mapping the relationships within the
economy was equivalent to studying human anatomy. . His input-output framework
Leontief (1936) views industries as nodes in a network of physical and
monetary flows. Conservation laws for these flows lead, at economic
equilibrium, to linear systems of equations linking the production of
different industries, whose solutions show how differences in the output of an
industry impact the output of any other economic sector.

They were used to determine, for example, how much one should invest in each
sector of an economy in order to increase the production of a given
sector.222This was indeed not an easy problem: to increase the production of
steel, it is necessary to increase the production of coal, but coal extraction
requires having steel. Input-output analysis provided tools to solve this
conundrum.. It was in particular an important tool for central planners in the
decades following the Second World War Bollard (2019).

Later on, input-output analysis was used to understand the origins of
macroeconomic fluctuations, with the seminal paper of Long and Plosser Long
and Plosser (1983), where the input-output network amplifies small shocks that
can lead to system-wide crises. However, most of these analyses are conducted
at a very coarse-grained level, in the sense that they attempt to model the
different sectors of the economy rather than modelling more granular
constituents: there are 405 industries in U.S. Bureau of Economic Analysis’
most disaggregated input-output tables, while there are approximately 200
million firms worldwide. This is an unsettling remark, as recent literature
Acemoglu  _et al._ (2012); Carvalho  _et al._ (2020); Diem  _et al._ (2022)
shows that fine-grained production networks play an important role in the
propagation of shocks, and that aggregating firms into sectors can lead to a
misestimation of risk and distress propagation. Detailed firm-level data will
also be crucial to the coming of age of agent-based modelling, a promising
approach to studying out-of-equilibrium macro-economic phenomena Dessertaine
_et al._ (2022), that recently matched the forecasting accuracy of more
traditional methods Poledna  _et al._ (2023); Hommes  _et al._ (2022); Pichler
_et al._ (2020).

Firm-level production data is thus very useful, but is also scarce Bacilieri
_et al._ (2022): the few datasets that are available only cover certain
countries or certain categories of companies, leaving most of the global
production network inaccessible. To tackle this problem, recent efforts have
attempted to reconstruct the production network, inferring the topology of the
network using only partial, aggregate or related data. For instance, Reisch
_et al._ (2022) uses mobile phone data to reconstruct the supply chain network
of an undisclosed European country, while Brintrup  _et al._ (2018) and
Kosasih and Brintrup (2021) pioneered machine learning for link prediction in
supply chains, leveraging topological features computed by hand or distilled
automatically through Graph Neural Networks. A similar approach was used in
Mungo  _et al._ (2022) to predict links between firms using their financial,
industrial, and geographical features. Additional efforts have been carried
out to adapt maximum-entropy models Hooijmaaijers and Buiten (2019); Mattsson
_et al._ (2021); Ialongo  _et al._ ; Squartini  _et al._ (2018); Almog  _et
al._ (2019); Squartini and Garlaschelli (2011); Squartini  _et al._ (2015),
already popular for models of international trade Squartini and Garlaschelli
(2014); Garlaschelli and Loffredo (2004, 2005); Garlaschelli  _et al._ (2007);
Almog  _et al._ (2019), to the reconstruction of firm-level networks. The
motivation of this research effort is that economic models conceived to
represent the economy at the firm level require a good knowledge of the
production network and should lead to a better understanding of economic
dynamics and forecast. But the converse should also be true: supply chains are
vital in a firm’s production, and they should leave a trace on the dynamics of
a firm, something that has been observed when considering natural disasters
Carvalho  _et al._ (2020) or the dynamics of companies’ market capitalisation
Abergel and Akar (2022). Is it possible to work backwards from this, and infer
the network topology from firm dynamics?

The study of firm dynamics, through the statistical analysis of their growth
rates, has a long history which dates back to the work of Gibrat Sutton
(1997). Gibrat’s model is a multiplicative growth model initially proposed to
explain the distribution of firm sizes (proxied, e.g., by sales or number of
employees). The model assumes that a firm grows by a random percentage of its
current size from one period to the next. This random variable is thought of
as being independent across firms, and was initially also modelled as having
the same distribution for all companies. Although this last hypothesis has
been weakened in past work, showing for example that the volatility of firm
growth decreases with their size in a non-trivial way Amaral  _et al._ (1997),
and even that it is necessary to think of the volatility of growth as being
firm-dependent Moran  _et al._ , the hypothesis of independence has not been
explicitly questioned thus far. We propose to go beyond this, making the
dependencies between firm growth explicit by studying the correlations between
them and leveraging this information to reconstruct the firm network.

This paper is organised as follows. Section I gives an overview of the data we
use for our paper, which we use in conjunction with the methods we outline in
Section II. Section III presents clear empirical evidence of the link between
the supply chain and firm growth. Section IV makes use of these observations
to reconstruct the production network from firm growth time series. We detail
both the optimisation algorithm used to carry out this reconstruction as well
as the results we obtain. Finally, Section V concludes.

##  I Data

The primary data sources used in this article are the FactSet Fundamentals and
FactSet Supply Chain Relationships datasets. Together, they provide a coherent
environment from which companies’ financial information (such as their
quarterly sales or market capitalisation), legal information (e.g., their
industrial classification or headquarters location) and supply chain
connections can be retrieved. Although it is very large, it should be noted
that this dataset has a strong bias in covering mainly US firms.

The first dataset contained in this environment, FactSet Fundamentals,
contains firms’ financial, balance sheet, and legal information. The dataset
spans a time range going from the early 1980s to the present day and covers
developed and emerging markets worldwide for a total of around $100,000$
active and inactive companies. From 1995 onwards, data on firms’ sales,
capitalisation and investments is available for each quarter.

The second dataset, FactSet Supply Chain Relationships, is assembled by
FactSet using multiple sources. The most prominent of these are filings
required by the US Federal Accounting Standards, whereby each firm must report
its most important suppliers and clients, and import-export declarations from
bills of lading. These sources are complemented with insight mined by FactSet
from news, press releases, company websites, and other sources of business
intelligence, which permit the inference of a link between two companies. Each
record of a link between two companies can be represented by a temporal
network, using directed links connecting a supplier to its customers. The
temporal dimension of this data is also provided by FactSet: each link is
assigned specific timestamps indicating the first time the connection was
reliably attested and when the connection is known to have ended, when this is
the case. 333Note that this procedure implies that persistent links appear
multiple times, as they are reported over many years.

To simplify our analysis, we have discarded the temporal dimension by
aggregating all the links into a single network that only considers whether a
link between two companies was ever present in the time period we consider.
Another simplification we perform is to aggregate firms that may be part of
large conglomerates at the ultimate parent level using ownership structure
data. Thus, the total sales, market capitalisation and any other balance sheet
data of these aggregated entities are the sum of these quantities for each of
the constituting entities. At the network level, this procedure has the effect
of deleting possible self-loops, as, for example, two branches of the same
conglomerate that are present in separate countries can trivially be reported
to have supply chain linkages between them. These aggregated entities
constitute what we understand by “firms” or “companies” in the remainder of
this paper.

Finally, we have only retained firms in the global supply chain’s weakly
largest connected component444 A weakly connected component is a set of nodes
such that for any two nodes $A$ and $B$, there exists a directed path starting
at $A$ and arriving at $B$ or from $B$ to $A$, but not necessarily the other
way around. When both a path $A\to\ldots\to B$ and $B\to\ldots\to A$ exist for
any two nodes $A$ and $B$ in the component, a much more restrictive condition,
then it is said to be strongly connected. , whose financial information was
available for at least eight years.555The reason for this is to remove time
series that are too short for our analysis, as the reader will appreciate
later. Our final sample is composed of $16,401$ firms connected by $178,911$
links. Appendix B details how to transform FactSet’s original tables into our
working dataset.

Number of firms | $16,401$
---|---
Number of links | $178,911$
Density | $6.7\times 10^{-4}$
Median degree | $7$
Max. degree | $1664$
Table 1: Network summary statistics

##  II Growth time series

We label firms with an index $i=1,\ldots,N$, calling $s_{i}(t)$ and $m_{i}(t)$
the sales and market capitalisation (the stock price multiplied by the number
of shares outstanding) of firm $i$ at time $t$ (counted in quarters). With
this, we define the annual growth rate of the sales of the firm as

$g_{i}(t):=\log\left(\frac{s_{i}(t+4)}{s_{i}(t)}\right).$ (1)

This quantity describes sales’ variations over the scale of a year, but
sampled with a quarterly frequency. We follow Ref. Moran  _et al._ in
describing sales growth rates with a random variable with a Gaussian central
region, although with fatter tails than a normal distribution, along with
firm-dependent mean and variance (volatility). This therefore leads us to
define the rescaled growth-rates,

$\begin{split}g^{\prime}_{i}(t):=\frac{g_{i}(t)-\mathbb{E}_{t^{\prime}}\left[g_{i}(t^{\prime})\right]}{\sqrt{\mathbb{V}_{t^{\prime}\neq
t}\left[g_{i}(t^{\prime})\right]}}\\\ \end{split}$ (2)

where the average is computed over all times $t^{\prime}$, but the variance is
computed from the time series where the observation corresponding to
$t^{\prime}=t$ has been removed. This corresponds to the leave-one-out
rescaling defined in Bouchaud and Potters (2003), where the denominator on the
right-hand side of Eqs. (2) allows one to rescale with respect to the
volatility when considering a variable with a fat-tailed
distribution.666Indeed, when the distribution is fat-tailed then the naive
estimator for the variance, related to $\sum_{t}g_{i}(t)^{2}$, may be
dominated by a single observation (the largest one in the sample) and
therefore introduce an artificial cut-off when dividing by the variance,
because in this case $\sum_{i}g_{i}(t)^{2}\approx\max_{t}g_{i}(t)^{2}$. When
rescaling the largest value in the sample, it is clear that it may be clipped
because of this. We drop the apostrophe below for clarity, as we will not use
the “bare” growth rates in the remainder of this article.

Our goal in the rest of this article is to infer the supply chain structure
from the correlation structure of the growth rates. Nonetheless, it is likely
that the growth rates of two companies are correlated because of reasons other
than their connection through the supply chain. This can be the case, for
instance, if two firms are in a given country that endures an exogenous
economic shock, as in the case of the Covid-19 pandemic. Our strategy
therefore will be to attempt to remove these common factors, assuming that
what remains in the correlations must be the more subtle effects due to the
supply chain. To illustrate the technique use for this, we shall resort to a
very simple model that is described below.

###  II.1 Removing common shocks

Let us propose first a very simple example, where one has $N$ time series
$x_{i}(t)$, with $1\leq i\leq N$ and $1\leq t\leq T$. Each time series
$x_{i}\left(t\right)$ is composed of an idiosyncratic term, driving time
series $i$ only and given by i.i.d. Gaussian terms, and a common term that
affects all the time series and that is also random. The model reads

$x_{i}(t)=\xi_{i}(t)+\sigma v(t),$ (3)

where $\xi_{i}(t)$ is a Gaussian random variable with
$\mathbb{E}[\xi_{i}(t)]=0$ and
$\mathbb{E}[\xi_{i}(t)\xi_{j}(t^{\prime})]=\delta_{ij}\delta_{tt^{\prime}}$,
with $\delta_{ij}$ the Kronecker delta (i.e., $\delta_{ij}=1$ if $i=j$ and $0$
otherwise). Similarly, $v(t)$ is a Gaussian random variable satisfying
$\mathbb{E}[v(t)v(t^{\prime})]=\delta_{tt^{\prime}}$ and
$\mathbb{E}[v(t)\xi_{i}(t^{\prime})]=0$.

In this case, where we know precisely the nature of the common shock, we can
estimate $v(t)$ when $N$ is large by writing:

$\frac{1}{N}\sum_{i=1}^{N}x_{i}(t)=\frac{1}{N}\sum_{i=1}^{N}\xi_{i}(t)+\sigma
v(t)\underset{N\gg 1}{\approx}\sigma v(t).$ (4)

The correlation matrix for the model’s time series reads

$\begin{split}C_{ij}&:=\mathbb{E}[\xi_{i}(t)\xi_{j}(t)]=\delta_{ij}+\sigma^{2},\end{split}$
(5)

which we can rewrite as
$\mathbf{C}=\mathbf{I}+N\sigma^{2}\mathbf{u}\mathbf{u}^{\intercal}$, with
$\mathbf{u}=\frac{1}{\sqrt{N}}\mathbf{1}$, and where $\mathbf{u}^{\intercal}$
indicates vector transposition.777This vector $\mathbf{u}$ is chosen to be
normalised. Because $\mathbf{C}$ is the sum of the identity matrix and a rank-
one matrix, it is easy to see that it has an eigenvalue $1+\sigma^{2}$,
corresponding to the eigenvector $\mathbf{u}$ as
$\mathbf{C}\mathbf{u}=(1+\sigma^{2})\mathbf{u}$, with all the other $N-1$
remaining eigenvalues equal to $1$, with eigenvectors corresponding to the
canonical basis of the vector space that is orthogonal to $\mathbf{u}$. We can
in fact go further in this geometric interpretation and bring meaning to the
vector $\mathbf{u}$ by focusing on the projection of the time series onto it.
What we mean by this is that for every time-step in the multi-dimensional time
series, we may consider the vector
$\mathbf{x}(t)=\left(x_{1}(t),\ldots,x_{2}(t)\right)$, and consider the
projected time series $\hat{v}(t)=\mathbf{u}\cdot\mathbf{x}(t)$.

In this case, we notice that for large $N$ we should have
$\hat{v}(t)=\frac{1}{N}\sum_{i=1}^{N}x_{i}(t)\approx\sigma v(t)$. We can
actually generalise this: if we replace Eq. (3) by

$x_{i}(t)=\xi_{i}(t)+\sigma u_{i}v(t),$ (6)

that is a model where each time series has a different exposure (or loading,
in factor-models’ jargon) to the common mode $v(t)$, then the correlation
matrix is the same and we still have an eigenvector
$\mathbf{u}=(u_{1},\ldots,u_{N})$.888This vector can be assumed to be
normalised, if not we can always replace $\sigma$ by
$\sqrt{\mathbf{u}^{2}}\sigma$ in the model. Doing the projection
$\mathbf{x}(t)\cdot\mathbf{u}(t)$ still leads to $\hat{v}(t)\approx v(t)$.

In fact, we can also consider the orthogonal projector to $\mathbf{u}$, given
by $\mathbf{P}=\mathbf{I}-\mathbf{u}\mathbf{u}^{\intercal}$, or equivalently
$P_{ij}=\delta_{ij}-u_{ij}$. We can now apply this projector to our time
series, as $\mathbf{y}(t)=\mathbf{P}\mathbf{x}$, or equivalently by defining
$\mathbf{Y}=\mathbf{P}\mathbf{X}$. It is straightforward to check that
$y_{i}(t)=x_{i}(t)-\hat{v}(t)\approx\xi_{i}(t)$.

To address our general problem of removing common fluctuations from time
series, we can adopt the following procedure to remove the common mode and be
left only with the idiosyncratic fluctuations. Assuming that the common mode
$v(t)$ is the primary driver of time series variations ($\sigma\gg 1$), we
can:

  1. 1.

Take the time series and compute the empirical correlation matrix,

  2. 2.

Diagonalise the correlation matrix and rank the eigenvalues and eigenvectors
according to the magnitude of the eigenvalue,

  3. 3.

Project the time series onto the eigenvector corresponding to the largest
eigenvalue to get the dynamics of the common mode,

  4. 4.

Remove the dynamics of the common mode from the time series by using the
orthogonal projector to the corresponding eigenvector.

Naturally, we can repeat this procedure and remove also the mode corresponding
to the second largest eigenvalue and so on, so that it is easily generalisable
to other, more complex situations than the one of Eq. (3) (see Fig. 1 for an
example where the common mode $v(t)$ is a sinusoidal wave).

Figure 1: (A) The time series $x_{i}\left(t\right)$ are created by adding a
sine wave and an idiosyncratic random noise. (B) The spectrum of the empirical
correlation matrix
$\widehat{C}_{ij}=\frac{1}{T}\sum_{t=1}^{T}x_{i}\left(t\right)x_{j}\left(t\right)$,
along with the random benchmark given by the Marčenko-Pastur distribution.
Note the presence of an eigenvector at $\lambda\approx 16$, beyond the random
benchmark (C). The eigenmode $\hat{v}(t)$, obtained by projecting the time
series onto the vector $\widehat{\mathbf{u}}$ corresponding to the largest
eigenvalue, tracks the collective oscillations of the system.

The issue, however, is that this relies on the assumption that the empirical
correlation matrix is a reliable estimator of the “true” underlying
correlation matrix from which the data is generated.999At least in this model.
In reality, when analysing time series with this point of view we are making
the more stringent assumption that the correlation structure of data is time-
invariant. Although there has been work to relax this assumption in e.g.
financial data Bongiorno  _et al._ (2021), these approaches are difficult, if
not impossible, to adopt to the time series we analyse because of their
relatively small length and sampling frequency. Naturally, this is not true,
and one expects some estimation error when the length of the time series $T$
is finite. In our toy model above, it is in fact possible to separate the
contribution of the idiosyncratic noise, as
$\widehat{\mathbf{C}_{0}}:=\frac{1}{T}\left(\boldsymbol{\xi}\boldsymbol{\xi}^{\intercal}\right)_{ij}=\frac{1}{T}\sum_{t=1}^{T}\xi_{i}(t)\xi_{j}(t)$.
Because the elements of $\boldsymbol{\xi}$ are i.i.d. Gaussian random
variables, this empirical correlation matrix is known as a Wishart matrix
Wishart (1928), and the statistical properties of its spectrum are known to be
determined by the Marčenko-Pastur distribution Marčenko and Pastur (1967). For
a more in-depth understanding of this and other links with random matrix
theory, we invite the reader to consult Potters and Bouchaud (2020), but we
will explain the main results we need below.

Because
$\widehat{\mathbf{C}_{0}}\underset{T\to\infty}{\longrightarrow}\mathbf{I}$, we
expect naturally that for large time series the spectrum of
$\widehat{\mathbf{C}_{0}}$ should be concentrated around $1$. In practice
however, because of measurement error, we don’t expect all of its eigenvalues
to be equal to $1$. Thus, we intuitively expect the full spectrum of
$\widehat{\mathbf{C}}$ to be constituted of $N-1$ eigenvalues close to $1$,
which constitute the contribution coming from $\mathbf{C}_{0}$, and a single
peaked eigenvalue close to $\sigma^{2}$, which is the contribution coming from
the dynamics of $v(t)$ that couples all of the $N$ time series. For the full
empirical correlation matrix $\widehat{\mathbf{C}}$, we also expect that the
eigenvector corresponding to its largest eigenvalue will satisfy,
$\widehat{\mathbf{u}}\approx\mathbf{u}$. However, the result of Marčenko-
Pastur is that in the limit where both $N,T\to\infty$, but with the ratio
$q=\frac{N}{T}$ fixed, the spectrum of $\mathbf{C}_{0}$ is concentrated in the
interval $(1-\sqrt{q},1+\sqrt{q})$, called the “bulk”, and may also have a
delta-peak at $0$ if $q<1$. For finite $N,T$ we also expect some eigenvalues
to be slightly out of this interval. This sheds light on why in practice
finding the common mode may be difficult: if, say, $\sigma$ is of the order of
$q$, then the eigenvalue “spike” at $1+\sigma^{2}$ will in fact be inside the
Marčenko-Pastur interval. This is linked to the so-called Baik-Ben Arous-Péché
(BBP) transition Baik  _et al._ (2005), and in this case it is not possible to
reconstruct the common mode.

We can indeed imagine that we run the model and execute the procedure
described above first for a value of $\sigma\gg q$, and then reduce $\sigma$
progressively until we reach $\sigma\approx q$. When diagonalising the
empirical correlation matrix $\widehat{\mathbf{C}}$ and considering the
eigenvector corresponding to its largest eigenvalue, $\widehat{\mathbf{u}}$,
this eigenvector will match the “true” eigenvector $\mathbf{u}$ when
$\sigma\gg q$, so that for example $\widehat{\mathbf{u}}\cdot\mathbf{u}\approx
1$. However, as $\sigma\to q$ this overlap will decrease, and the intuition
then is that when the outlier eigenvalue reaches the Marčenko-Pastur bulk,
then its associated eigenvector $\widehat{\mathbf{u}}$ cannot now reliably be
thought of as an estimator of $\mathbf{u}$, and will instead point in any
random direction. In this case $\mathbf{u}\cdot\widehat{\mathbf{u}}$ will be
of order $1/\sqrt{N}$ (see (Potters and Bouchaud, 2020, Section 14.2.2), and
also Allez  _et al._ (2014) for an intuition for this phenomenon using Dyson
Brownian motion). In this case, the usage of the projectors, or steps 3 and 4
of our procedure, will not lead to the identification of common modes.

The conclusion from this is that we are indeed capable of identifying common
factors in time series using this approach, but we must first make sure that
these modes correspond to eigenvalues of the correlation matrix that are not
compatible with a random benchmark. Indeed, the example above corresponds to
time series of equal length, where each entry of the time series is drawn at
random from a Gaussian distribution. In this case, the random benchmark for
the spectrum is determined by the Marčenko-Pastur distribution, as said above.
The case of our time series is, however, different since sales data is not
available for every company at any time. Growth time series can have different
starting points and lenghts, and the period over which one can compute their
correlation is different for any pair of firms. Our data therefore has a lot
of missing values, and two firms present in non-overlapping times for example
will be set two have a correlation of $0$. Another issue is that the growth-
rate distribution is not Gaussian, and has slightly heavier tails.
Understanding the correlation spectrum of heavy-tailed processes is feasible
(see for example Biroli  _et al._ (2007)), but very difficult to do for any
distribution.

We can nonetheless establish a random benchmark for the correlation spectrum
computationally and use it to identify eigenvalues indicating correlated
modes. We achieve this by creating a surrogate of the growth-rate time series
where the missing data structure is preserved and where the individual growth-
rates are drawn at random from their empirical distribution. This is similar
to the procedure used in Vodenska  _et al._ (2016), where the authors randomly
shuffle a time-series to benchmark the eigenvalues of correlation matrices
that can be distinguished from noise.

Figure 2 shows that the real correlation spectrum has several eigenvalues that
are beyond the bulk corresponding to the random benchmark, both on the left
and on the right side of the bulk. Note that the presence of negative
eigenvalues is a consequence of missing data, and is something that one does
not obtain for standard Wishart matrices. The largest eigenvalue corresponds
to the market mode, a collective trend shared by all the firms in the supply
chain. This collective mode concerns all firms, as shown by the fact that the
entries of the corresponding eigenvector have (roughly) all the same sign and
magnitude101010This is similar to the toy model presented in Section II.1..
Thus, this mode corresponds to a common factor in the economy, and all the
firms move coherently with it. Interpreting the modes corresponding to
eigenvalues outside the bulk is more challenging: contrary to what is observed
in the correlation structure of financial returns, we have not been able to
identify them with specific industrial sectors or geographies. Because we are
unable to give these eigenvector a clear interpretation, and since they could
potentially carry information about the production network, we have decided to
remove only the first eigenmode from the time series. In the rest of our
paper, we will refer to growth time series cleaned of the system’s first
eigenmode as “cleaned” time series $\tilde{g}_{i}\left(t\right)$, and to their
correlation as the “cleaned” correlation.111111We attract the reader’s
attention to the fact that we mean “cleaning” in a sense that is the opposite
of what is done for returns’ correlation matrices in finance: there, usually
one discards the modes corresponding to the smaller eigenvalues (see e.g. Bun
_et al._ (2017)). We, however, discard the largest mode because we want to
remove reasons for firm co-movement that are distinct from supply-chain
induced co-movement.

Figure 2: (A) The distribution $\rho\left(g\right)$ of the growth rates for
every firm $i$ and time $t$. A normal distribution is provided as a reference.
(B) Growth time series correlation spectrum. The two random benchmarks are
obtained by sampling random time series from the empirical distribution
$\rho\left(g\right)$ (Empirical benchmark) and the normal distribution
(Gaussian benchmark). The starting points and duration of the random time
series match those of the real ones. The spectrum shown is the average of $10$
sets of random time series.

##  III Network correlation and random benchmarks

We have introduced the main object of our analysis, firms’ growth time series
$\mathbf{g}_{i}(t)$. We will now show that the supply chain induces specific
correlations between firms, a necessary step to later justify our usage of
correlations in supply-chain reconstruction. We define the following
correlation matrices,121212Note that here we use the notation
$\mathbb{E}_{t}[\cdot]=\frac{1}{T}\sum_{t=1}^{T}\cdot(t)$ to indicate the
empirical average across the time variable. The notation $\mathbf{E}$ used in
the previous section corresponds instead to the “true” average value of our
stochastic model, computed over the distribution of the noise $\xi_{i}$ and
$v$. Similarly, $\mathbf{E}_{ij}$ indicates an empirical average taken by
summing over the variables $i$ and $j$.

$\begin{split}C_{ij}(\tau)&=\mathbb{E}_{t}\left[g_{i}(t)g_{j}(t+\tau)\right],\\\
\tilde{C}_{ij}(\tau)&=\mathbb{E}_{t}\left[\widetilde{g}_{i}(t)\tilde{g}_{j}(t+\tau)\right].\\\
\end{split}$ (7)

We can compute the average value of the elements of the matrix $\mathbf{C}$
and $\widetilde{\mathbf{C}}$ across the pairs of firms $\left(i,j\right)$
linked in the production network, defining averaged client/supplier
correlation functions. Given any (binary) adjacency matrix $\mathbf{A}$ we
define

$C_{\mathbf{A}}\left(\tau\right)=\mathbb{E}_{ij}\left[C_{ij}\left(\tau\right)|A_{ij}=1\right],$
(8)

and

$\widetilde{C}_{\mathbf{A}}\left(\tau\right)=\mathbb{E}_{ij}\left[\tilde{C}_{ij}\left(\tau\right)|A_{ij}=1\right],$
(9)

where the average runs over all pairs $1\leq i\leq j\leq N$. In other words,
$C_{\mathbf{A}}$ and $\widetilde{C}_{\mathbf{A}}$ are the average correlation
between two neighbours in a graph with an adjacency matrix $\mathbf{A}$. This
average can be computed using the true adjacency matrix of the production
network, $\mathbf{S}$, or over the adjacency matrix of any other network.

###  III.1 Random benchmarks

We first compute the correlations averaged over the adjacency matrix
$\mathbf{S}$ of FactSet’s production network, where $S_{ij}=1$ if $j$ either
supplies or is a client of $i$, and compare their value to those obtained with
several random network models: the Erdős-Rényi model Erdős and Rényi (1959),
the Stochastic Block Model Karrer and Newman (2011), and the Configuration
Model Newman (2003). We describe all three models and their parameters in
detail below.

We randomly sample $n=50$ networks of each model, with adjacency matrices
$\mathbf{R}_{1},\ldots,\mathbf{R}_{n}$ and compute the mean and standard
deviation of the sets
$\left\\{C_{\mathbf{R}_{1}},\ldots,C_{\mathbf{R}_{n}}\right\\}$ and
$\left\\{\tilde{C}_{\mathbf{R}_{1}},\ldots,\tilde{C}_{\mathbf{R}_{n}}\right\\}$.
All of the models are parametrised to match the empirical properties of the
supply-chain network.

For the Erdős-Rényi network, we fix its density $p$ to match that of the
production network, namely

$p=\frac{1}{N\left(N-1\right)}\sum_{i=1}^{N}\sum_{j>i}^{N}S_{ij}.$

The Erdős-Rényi network has no real structure, and in particular no clear
community structure is apparent in it. We therefore also used stochastic block
models, which we initialised with several different block schemes.
Specifically, we divided firms into blocks
$\left\\{B_{1},\ldots,B_{m}\right\\}$ depending on their industrial sector (at
their SIC code’s third-digit level of aggregation), their country, or their
network community as identified by the Louvain community-detection algorithm
Blondel  _et al._ (2008). The network densities within- and across- blocks are
chosen to be equal to their empirical counterparts,

$\rho_{ij}=\frac{1}{\left|B_{i}\right|\left(\left|B_{j}\right|-\delta_{ij}\right)}\sum_{k\in
B_{i},l\in B_{j}}A_{kl}.$ (10)

Finally, we use the configuration model to produce networks with a degree
distribution that matches exactly the empirical one.

Figure 3 compares the average correlation measured on the true production
network $\mathbf{S}$ and on the random network benchmarks. The value of
$C_{\mathbf{S}}\left(0\right)$ is twice as high as the average correlation
measured on the Erdős-Rényi graph, and $\approx 50\%$ higher than the
correlation measured for the configuration model. The result for
$\tilde{C}_{\mathbf{S}}\left(0\right)$ are even more striking, with the
residual correlation on the supply chain being still $\sim 0.1$ and most of
the random benchmarks dropping close to zero. This highlights the usefulness
of our cleaning procedure, as it significantly increases our signal-to-noise
ratio.

Figure 3: (A): Average correlation on the production network
$C\left(\tau\right)_{S}$ and several random network benchmarks. (B): Average
"cleaned" correlation on the production network
$\tilde{C}\left(\tau\right)_{S}$ and several random network benchmarks.
(C):Correlations along the supply chain decays with distance. At distance
$d=4$ ($d=3$ for the cleaned correlation), firms’ average correlation is the
same of the Erdos-Renyi benchmark. Results for the cleaned time series are
flagged with a (C)

###  III.2 Relationship with network distance

A second way to show that the supply chain induces correlations in the
dynamics of firm sales is to study how the correlation behaves with respect to
network distance. Intuitively, we expect that two firms that are close to each
other on the supply chain will be more correlated than two firms that are far
apart.

To see this, we start again from the binary adjacency matrix $\mathbf{S}$ of
the production network, and define recursively

$S_{ij}^{(k)}=\sum_{l_{1},\ldots,l_{k-1}}\mathbf{1}\left(S_{il_{1}}S_{l_{1}l_{2}}\ldots
S_{l_{k-1}j}>0\right)\prod_{m=1}^{k-1}\left(1-S_{ij}^{(m)}\right),$ (11)

where $S_{ij}^{(1)}=S_{ij}$. The first factor in the right hand side is equal
to $1$ if and only if there exists a path $i\to l_{1}\to\ldots\to j$ of length
$k$ linking $i$ to $j$. The second factor is $0$ if it exists a shorter path
from $i$ to $j$ in the network. Thus defined, $S^{(k)}_{ij}$ is equal to one
only if the shortest path between $i$ and $j$ is of length $k$.

Figure 4: An illustration of network distance. Nodes 2, 3 and 4 are at a
distance $k=1$ from node 10. Even though the path $1\to 3\to 4$ exists, we do
not consider $4$ to be at distance $k=2$ from $1$

We can see how these correlations decay with distance, by computing the values

$D_{S}(k)=\mathbb{E}_{ij}\left[C_{ij}(0)|S_{ij}^{(k)}=1\right],$ (12)

and

$\widetilde{D}_{S}(k)=\mathbb{E}_{ij}\left[\tilde{C}_{ij}(0)|S_{ij}^{(k)}=1\right],$
(13)

namely the average of the non-lagged growth correlation between any two firms
that are $k$-steps apart in the supply chain. We show this on Figure 3, C. The
correlation between firms decays as their distance in the production networks
increases, revealing again that the production networks mediates growth
correlations between firms.

##  IV Supply Chain Reconstruction

In the previous Sections, we have established that the supply-chain induces
correlations between firms, and we have also established that our cleaning
procedure increases the signal-to-noise ratio of these correlations with
respect to the real supply chain. We next propose a procedure to reconstruct
the supply chain using the cleaned correlation matrix.

Inferring networks from observations, or graph learning Dong  _et al._ (2019),
is a problem that encompasses several branches of natural and social sciences.
Following Dong  _et al._ (2019), we define the problem of graph learning as
follows: given $T$ observations on $N$ entities, represented by a data matrix
$\mathbf{X}\in\mathbb{R}^{N\times T}$, and taking some prior knowledge as
given, we seek to infer relationships between our $N$ entities and represent
these relationships as a graph $\mathcal{G}$.

A possible approach to solve this problem is to assume that $\mathcal{G}$
encodes some statistical relationship between the entities. Specifically,
probabilistic graphical models assume that the structure of $\mathcal{G}$
determines the joint probability distribution of the observations on the data
entities: the presence or absence of edges in the graphs encodes the
conditional independence among the random variables represented by the
vertices. In particular, Markov Random Fields consider a graph
$\mathcal{G}=\\{\mathcal{V},\mathcal{E}\\}$ and a set of random variables
$\mathbf{x}=\left\\{x_{i}:v_{i}\in\mathcal{V}\right\\}$ satisfying the
pairwise Markov property,

$\left(v_{i},v_{j}\right)\notin\mathcal{E}\Leftrightarrow
p\left(x_{i}|x_{j},\mathbf{x}\backslash\left\\{x_{i},x_{j}\right\\}\right)=p\left(x_{i},\mathbf{x}\backslash\left\\{x_{i},x_{j}\right\\}\right),$
(14)

which simply states that two variables $x_{i}$ and $x_{j}$ are conditionally
independent if there is no edge between the corresponding vertices $v_{i}$ and
$v_{j}$. In Markov Random Fields, the joint probability distribution of the
variables $x_{1},\ldots,x_{N}$ may also be represented as

$p\left(\mathbf{x}\right)=\frac{1}{Z}\prod_{i=1}^{K}\phi_{i}\left(\mathbf{D}_{i}\right),$
(15)

where $\mathbf{D}_{i},\ldots,\mathbf{D}_{K}$ are a set of graph’s cliques
(i.e., groups of nodes), $Z$ is a normalisation factor known as partition
function, and $\phi_{i}$s are generic functions known as factors. It is
straightforward to see that the exponential family of distributions with a
parameter matrix $\mathbf{\Theta}\in\mathbb{R}$,

$p\left(\mathbf{x}|\mathbf{\Theta}\right)=\frac{1}{Z\left(\mathbf{\Theta}\right)}\exp\left(\sum_{v_{i}\in\mathcal{V}}\theta_{ii}x_{i}^{2}+\sum_{(v_{i},v_{j})\in\mathcal{E}}\theta_{ij}x_{i}x_{j}\right),$
(16)

is compatible with this formalism; the multivariate Gaussian distribution with
precision matrix $\mathbf{\Theta}$,

$p\left(\mathbf{x}|\mathbf{\Theta}\right)=\frac{\left|\mathbf{\Theta}\right|^{1/2}}{\left(2\pi\right)^{N/2}}\exp\left(-\frac{1}{2}\mathbf{x}^{T}\mathbf{\Theta}\mathbf{x}\right),$
(17)

belongs to this family. The subclass of Markov random fields that adopt Eq. 17
as the parametrisation for the joint probability distribution $p$ are called
Gaussian Markov Random Fields or Gaussian Graphical Models. In Gaussian
Graphical models, the problem of finding the graph $\mathcal{G}$ is reduced to
that of estimating a precision matrix $\mathbf{\Theta}$ that encodes the
conditional relationship between the nodes. In the previous section, we saw
that the production network influences the correlation of firms’ growth
$g_{i}$. If we consider each vector $\mathbf{g}\left(t\right)$ as a drawn from
a joint probability distribution where the correlations are driven by the
supply chain, Gaussian graphical models seem well equipped to reconstruct the
production network if one ignore the fact that the growth rates do not have a
Gaussian distribution.131313Indeed, the marginal distribution of $x_{i}$ in
Eq. (17) is clearly a Gaussian distribution. We think nonetheless that,
because the growth-rates show a Gaussian-like central region, as shown by
Moran  _et al._ , it is reasonable to use this model to attempt a
reconstruction.

We therefore use the Graphical Lasso method to construct an estimator
$\widehat{\mathbf{\Theta}}$ of $\mathbf{\Theta}$ by solving the following
optimisation problem:141414 This is the result of applying Bayes theorem
assuming a constant prior for $\mathbf{\Theta}$.

$\widehat{\mathbf{\Theta}}=\text{argmax}_{\mathbf{\Theta}}\log\det\mathbf{\Theta}-\mathrm{tr}\left(\widehat{\mathbf{C}}\mathbf{\Theta}\right)-\alpha\|\mathbf{\Theta}\|_{1},$
(18)

with $\widehat{\mathbf{C}}=\frac{1}{T}\mathbf{G}\mathbf{G}^{T}$ the sample
covariance matrix, $\det\left(\cdot\right)$ the determinant and
$\mathrm{tr}\left(\cdot\right)$ the trace. The first two terms can be thought
of as the log-likelihood of $\mathbf{\Theta}$ in the Gaussian Graphical Model,
while $\alpha\left|\mathbf{\Theta}\right|$ is an $L^{1}$ regularisation term
with parameter $\alpha$. This approach will,in general, recover a matrix
$\mathbf{\Theta}$ with both positive and negative entries. In this setting, a
positive off-diagonal entry $\theta_{ij}$ of the precision matrix implies a
negative partial correlation between $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$,
whose interpretation is problematic since we would like $\mathbf{\Theta}$ to
proxy the adjacency matrix of the network.

References Lake and Tenenbaum (2010); Daitch  _et al._ (2009); Hu  _et al._
(2015) suggest instead searching for the precision matrix among the set
$\mathcal{S}_{\mathbf{\Theta}}$ of possible Graph Laplacian matrices,

$\mathcal{S}_{\mathbf{\Theta}}=\left\\{\mathbf{\Theta}|\theta_{ij}=\theta_{ji}<0\
\text{for}\ i\neq j,\theta_{ii}=-\sum_{j\neq i}\theta_{ij}\right\\}.$ (19)

Conditioning $\widehat{\mathbf{\Theta}}$ to be in the set of possible graph
Laplacians has two interesting consequences. First, the graph Laplacian
$\mathbf{L}$ uniquely determines the adjacency matrix $\mathbf{W}$ of the
graph; thus, the problem in 18 with the assumption
$\mathbf{\Theta}\in\mathcal{S}_{\mathbf{\Theta}}$ creates a direct connection
between the data and the topology of the network. Second, since the time
series $\mathbf{g}_{i}$ has zero mean, we can write the trace
$\left(\widehat{\mathbf{C}}\mathbf{\Theta}\right)$ as

$\mathrm{tr}\left(\widehat{\mathbf{C}}\mathbf{\Theta}\right)=\frac{1}{T}\mathrm{tr}\left(\mathbf{G}\mathbf{G}^{T}\mathbf{\Theta}\right)=\frac{1}{T}\sum_{i,j}\sum_{t=1}^{T}\theta_{ij}\left(g_{i}(t)-g_{j}(t)\right)^{2}.$
(20)

The term on the right hand of the equation measures the (squared) difference
between the observation on firms $i$ and $j$ ($\mathbf{g}_{i}$ and
$\mathbf{g}_{j}$), computed over couples of connected firms ($\theta_{ij}>0$);
it is generally known as the quadratic energy function and quantifies the
smoothness of $\mathbf{G}$ over the graph with Laplacian L. For an economic
interpretation, the second term in 18,
$\mathrm{tr}\left(\widehat{\mathbf{C}}\mathbf{\Theta}\right)$, can be
interpreted as a penalty term affecting networks over which $\mathbf{G}$ is
not smooth, i.e., a productio network that exhibits large differences between
the growth rates of connected firms.

In Kumar  _et al._ (2019) (see Appendix A), the authors propose an efficient
algorithm to solve the problem in Eq. 18 while also enforcing some (soft)
constraints on the spectrum $\text{Sp}(\mathbf{\Theta})$ of the Laplacian
matrix. The problem becomes

$\begin{split}\widehat{\mathbf{\Theta}}=&\text{argmax}_{\mathbf{\Theta}}\log\det\mathbf{\Theta}-\mathrm{tr}\left(\widehat{\mathbf{C}}\mathbf{\Theta}\right)-\alpha\|\mathbf{\Theta}\|_{1},\\\
&\textrm{subject to}\quad\mathbf{\Theta}\in\mathcal{S}_{\mathbf{\Theta}},\
\text{Sp}(\mathbf{\Theta})\subset\mathcal{S}_{\lambda}\\\ \end{split}$ (21)

where $\mathcal{S}_{\lambda}$ is the set of admissible spectra that we choose.
Because the spectrum of the Laplacian encodes information about the underlying
network’s topology, choosing $\mathcal{S}_{\lambda}$ appropriately allows us
to enforce high-level topological features on the reconstructed network.

We therefore attempt to use the algorithm provided in Kumar  _et al._ (2019)
to reconstruct the production network. In the following, we assume that we
know the network’s density in advance, and that we also have a reliable
estimate for the number of links within and across different sectors. This
information would not be available directly in a real-world situation, but the
literature on production networks and other available data sources as input-
output tables allow informed guesses (see, e.g., Bacilieri  _et al._ (2022)).
This means that our results should be placed halfway between a proof of
concept and a realistic use case.

Figure 5: (A) A stylised representation of an adjacency matrix with two
sectors. The density of links between the $n_{A}$ firms in sector A is
$\rho_{A}$, the density of links between the $n_{B}$ firms in sector B is
$\rho_{B}$, and the density of links across the two sectors is $\rho_{AB}$.
(B) Another adjacency matrix. There are two group of firms of size $n_{A}$
(right bottom corner of the matrix) and $n_{B}$ (top left corner of the
matrix). The density within firms in the first group is $\rho_{B}$, the
density between firms in the second group is $\rho_{A}$, and the density
across the groups is $\rho_{AB}$. The graph Laplacian of the matrix in (A) and
that of the matrix in (B) will have the same spectrum. However, the density
within and across sectors in (B) is different from that in (A).

We must however slightly modify this algorithm to apply it to our specific
situation. Indeed, a problem with the algorithm described in Kumar  _et al._
(2019) is that, while it is possible to encode a given community structure by
constraining the Laplacian, we are not be able to specify which firms should
go into which community (see Fig. 5).

To solve this, we have devised the following procedure. First, we split
$\widehat{\mathbf{C}}$ into diagonal and off-diagonal blocks based on firm
industries. Next we use the procedure defined in 21 to reconstruct each
diagonal block independently. Thirdly, we go through all the possible pairs of
diagonal blocks and – keeping the diagonal blocks equal to those that were
reconstructed in the previous step – we reconstruct the off-diagonal blocks.
Finally, we assemble all the blocks together to obtain the entire adjacency
matrix; this procedure is shown graphically in Fig. 6.

Figure 6: Reconstruction of the supply chain networks. The original
correlation matrix (A) is split in the different industry sectors. First, we
reconstruct the diagonal blocks (B). Then, we reconstruct the off-diagonal
blocks (C). Finally, we re-assemble the blocks together (D).

Every time we reconstruct a network, we choose the parameter $\alpha$ to match
the empirical network density. To reconstruct the diagonal blocks, we use the
spectrum obtained by averaging over the spectra 1000 Erdős-Rényi random
networks’ Laplacians, with probability $p$ equal to the desired density.
Similarly, to reconstruct the off-diagonal blocks, we use the spectrum
obtained by averaging over the spectra of 1000 block models’ Laplacians, where
the probabilities of links within and across each block are chosen to match
the desired density. We provide details on the reconstruction algorithm in
Appendix A.

We ran our procedure over several different subparts of the real production
network, each composed of a minimum of 300 to a maximum of 500 firms. We
compared our results to those of two random benchmarks: a Erdos-Renyi graph
and an industrial sector block model, built as in III.

Figure 7: True Positive rate (left), Accuracy (middle), and F1 Score (right)of
the reconstructed networks, plotted against the same metrics for the two
different random benchmarks.

Our approach beats the benchmark across several different metrics, often with
a considerable performance difference (Fig. 7).

##  V Conclusions

In this paper, we studied if the correlation between firms’ growth time series
could be useful in reconstructing production networks. Using FactSet’s supply
chain network as a use case and several random network models as benchmarks,
we have first shown that growths of firms connected in the production networks
are on average more correlated than those of randomly selected firms’ pairs.
We have shown that this effects fades gradually as one looks at the average
correlation between pair of firms at an increasing network distance along the
supply chain. Finally, we have framed the production network reconstruction in
the context of graph learning, and showed that some recent techniques
developed in the field can successfully be used to identify trade connections
between firms. Even if the results of our paper are closer to a proof of
concept than to a realistic use case, we believe that our approach is
promising. First, it relies on a mechanism that can be easily accepted as
universal: the growth of business partners is correlated. Second, it is a
fully "unsupervised" approach, which does not require the training of a model,
and is not prone to over-fitting. Third, it requires data that is easily
accessible (firms’ sales) and, to a certain extent, substitutable (e.g., we
obtained similar results when we looked at the correlation of firms’ stock
returns). Finally, it generates a network that matches a set of desired
topological features. This last point also highlights interesting avenues of
research: as more "universal" production networks features will be documented,
and better generative models for these networks will be developed, the more
effective our approach will be.

## Acknowledgements

We would like to thank Jean-Philippe Bouchaud, François Lafond, Doyne Farmer,
and Xiaowen Dong for their numerous suggestions for this work, and Andrea
Bacilieri for her help in handling the data. We would also like to thank the
participants of the 2022 CSH–INET Workshop on Firm-Level Production Networks
and the CCS 2022, in particular Christian Diem and Tobias Reisch, for the
useful feedback, and Stefan Thurner and the network economics group at CSH
Vienna for their hospitality and insight. This work was supported by Baillie
Gifford and the Institute for New Economic Thinking at the Oxford Martin
School.

## References

  * Quesnay (1766) François Quesnay, “Analyse de la formule arithmétique du tableau économique de la distribution des dépenses annuelles d’une nation agricole,” Journal d’agriculture, du commerce et des finances 5, 11–41 (1766).
  * Leontief (1936) Wassily W. Leontief, “Quantitative input and output relations in the economic systems of the united states,” The Review of Economics and Statistics 18, 105 (1936).
  * Bollard (2019) Alan Bollard, “The Peacenik who Helped Bombing Tactics: Wassily Leontief in the USA, 1943–4,” in  _Economists at War: How a Handful of Economists Helped Win and Lose the World Wars_ (Oxford University Press, 2019) https://academic.oup.com/book/0/chapter/321625127/chapter-ag-pdf/44488337/book_36631_section_321625127.ag.pdf .
  * Long and Plosser (1983) John B. Long and Charles I. Plosser, “Real business cycles,” Journal of Political Economy 91, 39–69 (1983).
  * Acemoglu  _et al._ (2012) Daron Acemoglu, Vasco Carvalho, Asu Ozdaglar, and Alireza Tahbaz-Salehi, “The network origins of aggregate fluctuations,” Econometrica 80, 1977–2016 (2012).
  * Carvalho  _et al._ (2020) Vasco M Carvalho, Makoto Nirei, Yukiko U Saito, and Alireza Tahbaz-Salehi, “Supply Chain Disruptions: Evidence from the Great East Japan Earthquake*,” The Quarterly Journal of Economics 136, 1255–1321 (2020), https://academic.oup.com/qje/article-pdf/136/2/1255/36725306/qjaa044.pdf .
  * Diem  _et al._ (2022) Christian Diem, András Borsos, Tobias Reisch, János Kertész, and Stefan Thurner, “Quantifying firm-level economic systemic risk from nation-wide supply networks,” Scientific reports 12, 1–13 (2022).
  * Dessertaine  _et al._ (2022) Théo Dessertaine, José Moran, Michael Benzaquen, and Jean-Philippe Bouchaud, “Out-of-equilibrium dynamics and excess volatility in firm networks,” Journal of Economic Dynamics and Control 138, 104362 (2022).
  * Poledna  _et al._ (2023) Sebastian Poledna, Michael Gregor Miess, Cars Hommes, and Katrin Rabitsch, “Economic forecasting with an agent-based model,” European Economic Review 151, 104306 (2023).
  * Hommes  _et al._ (2022) Cars Hommes, Mario He, Sebastian Poledna, Melissa Siqueira, and Yang Zhang, “Canvas: A canadian behavioral agent-based model,”  (2022), 10.34989/SWP-2022-51.
  * Pichler  _et al._ (2020) Anton Pichler, Marco Pangallo, R. Maria del Rio-Chanona, François Lafond, and J. Doyne Farmer,  _Production networks and epidemic spreading: How to restart the UK economy?_, INET Oxford Working Papers 2020-12 (Institute for New Economic Thinking at the Oxford Martin School, University of Oxford, 2020).
  * Bacilieri  _et al._ (2022) Andrea Bacilieri, András Borsos, Pablo Astudillo-Estevez, and François Lafond, “Firm-level production networks: what do we (really) know?” mimeo, University of Oxford (2022).
  * Reisch  _et al._ (2022) Tobias Reisch, Georg Heiler, Christian Diem, Peter Klimek, and Stefan Thurner, “Monitoring supply networks from mobile phone data for estimating the systemic risk of an economy,” Scientific Reports 12, 13347 (2022), number: 1, Publisher: Nature Publishing Group.
  * Brintrup  _et al._ (2018) A. Brintrup, P. Wichmann, P. Woodall, D. McFarlane, E. Nicks, and W. Krechel, “Predicting hidden links in supply networks,” Complexity 2018, 9104387 (2018).
  * Kosasih and Brintrup (2021) Edward Kosasih and Alexandra Brintrup, “A machine learning approach for predicting hidden links in supply chain with graph neural networks,” International Journal of Production Research (2021).
  * Mungo  _et al._ (2022) Luca Mungo, François Lafond, Pablo Astudillo-Estevez, and J. Doyne Farmer,  _Reconstructing production networks using machine learning_ , Tech. Rep. 2 (Institute for New Economic Thinking, 2022).
  * Hooijmaaijers and Buiten (2019) Sjoerd Hooijmaaijers and Gert Buiten,  _A methodology for estimating the Dutch interfirm trade network, including a breakdown by commodity_, Tech. Rep. (Technical report, Statistics Netherlands, 2019).
  * Mattsson  _et al._ (2021) Carolina E. S. Mattsson, Frank W. Takes, Eelke M. Heemskerk, Cees Diks, Gert Buiten, Albert Faber, and Peter M. A. Sloot, “Functional structure in production networks,” Frontiers in Big Data 4, 23 (2021).
  * (19) Leonardo Niccolò Ialongo, Camille de Valk, Emiliano Marchese, Fabian Jansen, Hicham Zmarrou, Tiziano Squartini, and Diego Garlaschelli, “Reconstructing firm-level interactions in the dutch input–output network from production constraints,”  12, 11847, number: 1 Publisher: Nature Publishing Group.
  * Squartini  _et al._ (2018) Tiziano Squartini, Guido Caldarelli, Giulio Cimini, Andrea Gabrielli, and Diego Garlaschelli, “Reconstruction methods for networks: The case of economic and financial systems,” Physics Reports 757, 1–47 (2018).
  * Almog  _et al._ (2019) Assaf Almog, Rhys Bird, and Diego Garlaschelli, “Enhanced gravity model of trade: Reconciling macroeconomic and network models,” Frontiers in Physics 7, 55 (2019).
  * Squartini and Garlaschelli (2011) Tiziano Squartini and Diego Garlaschelli, “Analytical maximum-likelihood method to detect patterns in real networks,” New Journal of Physics 13, 083001 (2011).
  * Squartini  _et al._ (2015) Tiziano Squartini, Rossana Mastrandrea, and Diego Garlaschelli, “Unbiased sampling of network ensembles,” New Journal of Physics 17, 023052 (2015).
  * Squartini and Garlaschelli (2014) Tiziano Squartini and Diego Garlaschelli, “Jan Tinbergen’s legacy for economic networks: From the gravity model to quantum statistics,” in  _Econophysics of Agent-Based Models_, edited by Frédéric Abergel, Hideaki Aoyama, Bikas K. Chakrabarti, Anirban Chakraborti, and Asim Ghosh (Springer International Publishing, Cham, 2014) pp. 161–186.
  * Garlaschelli and Loffredo (2004) Diego Garlaschelli and Maria I. Loffredo, “Fitness-dependent topological properties of the world trade web,” Phys. Rev. Lett. 93, 188701 (2004).
  * Garlaschelli and Loffredo (2005) Diego Garlaschelli and Maria I. Loffredo, “Structure and evolution of the world trade network,” Physica A: Statistical Mechanics and its Applications 355, 138–144 (2005), market Dynamics and Quantitative Economics.
  * Garlaschelli  _et al._ (2007) D. Garlaschelli, T. Di Matteo, T. Aste, G. Caldarelli, and M. I. Loffredo, “Interplay between topology and dynamics in the world trade web.” The European Physical Journal B 57, 159 – 164 (2007).
  * Abergel and Akar (2022) Frédéric Abergel and Adrien Akar, “Supply chain and correlations,” The Journal of Portfolio Management 49, 138–158 (2022).
  * Sutton (1997) John Sutton, “Gibrat’s legacy,” Journal of Economic Literature 35, 40–59 (1997).
  * Amaral  _et al._ (1997) Luís Amaral, Sergey Buldyrev, Shlomo Havlin, Heiko Leschhorn, Philipp Maass, Michael Salinger, H. Stanley, and Raymond Stanley, “Scaling behavior in economics: I. empirical results for company growth,” http://dx.doi.org/10.1051/jp1:1997180 7 (1997), 10.1051/jp1:1997180.
  * (31) José Moran, Angelo Secchi, and Jean-Philippe Bouchaud, In preparation.
  * Bouchaud and Potters (2003) Jean-Philippe Bouchaud and Marc Potters,  _Theory of Financial Risk and Derivative Pricing_ (Cambridge University Press, 2003).
  * Bongiorno  _et al._ (2021) Christian Bongiorno, Damien Challet, and Grégoire Loeper, “Cleaning the covariance matrix of strongly nonstationary systems with time-independent eigenvalues,”  (2021), 10.48550/ARXIV.2111.13109.
  * Wishart (1928) John Wishart, “The generalised product moment distribution in samples from a normal multivariate population,” Biometrika 20A, 32–52 (1928).
  * Marčenko and Pastur (1967) V A Marčenko and L A Pastur, “DISTRIBUTION OF EIGENVALUES FOR SOME SETS OF RANDOM MATRICES,” Mathematics of the USSR-Sbornik 1, 457–483 (1967).
  * Potters and Bouchaud (2020) Marc Potters and Jean-Philippe Bouchaud,  _A First Course in Random Matrix Theory_ (Cambridge University Press, 2020).
  * Baik  _et al._ (2005) Jinho Baik, Gérard Ben Arous, and Sandrine Péché, “Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices,” The Annals of Probability 33 (2005), 10.1214/009117905000000233.
  * Allez  _et al._ (2014) Romain Allez, Joël Bun, and Jean-Philippe Bouchaud, “The eigenvectors of gaussian matrices with an external source,”  (2014).
  * Biroli  _et al._ (2007) G Biroli, J.-P Bouchaud, and M Potters, “On the top eigenvalue of heavy-tailed random matrices,” Europhysics Letters (EPL) 78, 10001 (2007).
  * Vodenska  _et al._ (2016) Irena Vodenska, Hideaki Aoyama, Yoshi Fujiwara, Hiroshi Iyetomi, and Yuta Arai, “Interdependencies and causalities in coupled financial networks,” PLOS ONE 11, e0150994 (2016).
  * Bun  _et al._ (2017) Joël Bun, Jean-Philippe Bouchaud, and Marc Potters, “Cleaning large correlation matrices: Tools from random matrix theory,” Physics Reports 666, 1–109 (2017), cleaning large correlation matrices: tools from random matrix theory.
  * Erdős and Rényi (1959) P. Erdős and A. Rényi, “On random graphs I,” Publicationes Mathematicae Debrecen 6, 290–297 (1959).
  * Karrer and Newman (2011) Brian Karrer and M. E. J. Newman, “Stochastic blockmodels and community structure in networks,” Phys. Rev. E 83, 016107 (2011).
  * Newman (2003) M. E. J. Newman, “The structure and function of complex networks,” SIAM Review 45, 167–256 (2003), https://doi.org/10.1137/S003614450342480 .
  * Blondel  _et al._ (2008) Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre, “Fast unfolding of communities in large networks,” Journal of Statistical Mechanics: Theory and Experiment 2008, P10008 (2008).
  * Dong  _et al._ (2019) Xiaowen Dong, Dorina Thanou, Michael Rabbat, and Pascal Frossard, “Learning graphs from data: A signal representation perspective,” IEEE Signal Processing Magazine 36, 44–63 (2019).
  * Lake and Tenenbaum (2010) B. Lake and K. Tenenbaum, “Discovering structure by learning sparse graphs,” in  _Proceedings of the 32nd Annual Meeting of the Cognitive Science Society (CogSci)_ (Cognitive Science Society (CogSci), Portland, OR, 2010) pp. 778–784.
  * Daitch  _et al._ (2009) Samuel I. Daitch, Jonathan A. Kelner, and Daniel A. Spielman, “Fitting a graph to vector data,” in  _Proceedings of the 26th Annual International Conference on Machine Learning_, ICML ’09 (Association for Computing Machinery, New York, NY, USA, 2009) p. 201–208.
  * Hu  _et al._ (2015) Chenhui Hu, Lin Cheng, Jorge Sepulcre, Keith A. Johnson, Georges E. Fakhri, Yue M. Lu, and Quanzheng Li, “A spectral graph regression model for learning brain connectivity of alzheimer’s disease,” PLOS ONE 10, 1–24 (2015).
  * Kumar  _et al._ (2019) Sandeep Kumar, Jiaxi Ying, Jose Vinicius de Miranda Cardoso, and Daniel Palomar, “Structured graph learning via laplacian spectral constraints,” in  _Advances in Neural Information Processing Systems_, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Kumar  _et al._ (2020) Sandeep Kumar, Jiaxi Ying, José Vinícius de M. Cardoso, and Daniel P. Palomar, “A unified framework for structured graph learning via spectral constraints,”  21, 1–60 (2020).

##  Appendix A Network Reconstruction Algorithm

The algorithm used to solve the problem in Eq. 21 has first been proposed by
Kumar  _et al._ (2019, 2020) in the context of structured Graph Learning. The
authors formulate the problem as the following. Let
$\boldsymbol{x}=\left[x_{1},x_{2},\ldots,x_{p}\right]^{T}$ be a
$p$-dimensional, zero-mean, random vector (in the practical case, this would
be the collection of the "cleaned" time series
$\tilde{\mathbf{g}}_{1},\ldots,\tilde{\mathbf{g}}_{N}$) associated with an
undirected graph $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$, where
$\mathcal{V}=\left\\{1,2,\ldots,p\right\\}$ is a set of nodes corresponding to
the elements of $\boldsymbol{x}$, and
$\mathcal{E}\in\mathcal{V}\times\mathcal{V}$ is the set of edges connecting
nodes. In the Gaussian Graphical modelling framework, learning a graph
corresponds to solving the optimtime series problem

$\max_{\Theta\in\mathcal{S}^{p}_{++}}\log\det\left(\Theta\right)-\text{tr}\left(\Theta
S\right)-\alpha h\left(\Theta\right),$ (22)

where $\Theta\in\mathbb{R}^{p\times p}$ denotes the desired graph matrix,
$\mathcal{S}^{p}_{++}$ denotes the set of $p\times p$ positive definite
matrices, $S\in\mathbb{R}^{p\times p}$ is the covariance matrix obtained from
the data, $S=\frac{1}{n}\boldsymbol{x}^{T}\boldsymbol{x}$,
$h\left(\cdot\right)$ is a generic regularisation term, and $\alpha$ is a
coefficient tuning the strength of the regularisation. As we saw in IV, a
matrix $\Theta\in\mathbb{R}^{p\times p}$ is called a combinatorial graph
Laplacian matrix if it belongs to the set

$\mathcal{S}_{\Theta}=\left\\{\Theta|\theta_{ij}=\theta_{ji}<0\ \text{for}\
i\neq j,\theta_{ii}=-\sum_{j\neq i}\theta_{ij}\right\\}.$ (23)

The Laplacian Matrix $\Theta$ is a symmetric, positive semi definite matrix
with zero row sums. In the framework of network theory, a Laplacian matrix
$\Theta$ is computed from a graph’s adjacency matrix $A$ as $\Theta=D-A$,
where $D$ is a diagonal matrix and $D_{ii}$ is the degree of node $i$. It is
straightforward to see that the adjacency matrix of a graph can be recovered
from the Laplacian matrix simply as $A=\Theta\odot\left(I-\mathbb{1}\right)$,
where $I$ is the identity matrix, $\mathbb{1}_{ij}=1$, and $\odot$ is the
element-wise product. The structural properties of a graph are encoded in the
eigenvalues of its Laplacian, so that being able to constraint the spectrum of
the matrix $\Theta$ in the optimisation problem in Eq. 22 allows to enforce
some structural constraints on the reconstructed network. The goal hence
becomes that of solving the problem

$\displaystyle\max_{\Theta}$
$\displaystyle\log\text{gdet}\Theta-\mathrm{tr}\left(S\Theta\right)-\alpha
h\left(\Theta\right),$ (24) subject to
$\displaystyle\Theta\in\mathcal{S}_{\Theta},\
\mathbf{\lambda}\left(\Theta\right)\in\mathcal{S}_{\mathbf{\lambda}},$

where $\text{gdet}\left(\Theta\right)$ denotes the generalised
determinant151515 Note that in the main text we have not made explicit the
difference between $\text{gdet}\left(\Theta\right)$ and
$\det\left(\Theta\right)$ to improve readability.  of the matrix $\Theta$,
defined as the product of its non-zero eigenvalues,
$\mathbf{\lambda}\left(\Theta\right)$ denotes the set of eigenvalues of
$\Theta$, and $\mathcal{S}_{\lambda}$ is the set containing the spectral
constraints on the eigenvalues. As the authors in Kumar  _et al._ (2019) point
out, from the probabilistic perspective, if the data is generated from a
multivariate Gaussian distribution
$\mathcal{N}\left(0,\Theta^{\dagger}\right)$, then Eq. 24 can be viewed as a
penalised maximum likelihood estimation of the structured precision matrix of
an attractive Gaussian Markov Random Field model, while, if $\boldsymbol{x}$
is arbitrarily distributed, the problem in Eq. 24 corresponds to minimizing a
penalised log-determinant Bregman divergence (a common measure of distance for
probability distributions), and hence its solution should anyway result in a
meaningful graph. In the main body of the paper we saw how we assume to know
the spectrum $\bar{\boldsymbol{\lambda}}$ of the target matrix is known, so we
can define $\mathcal{S}_{\lambda}$ as

$\mathcal{S}_{\boldsymbol{\lambda}}=\left\\{\lambda_{i}=\bar{\lambda}_{i},\
\forall i\in\left[1,p\right]\right\\}.$ (25)

To solve the optimisation problem in Eq. 24, the authors in Kumar  _et al._
(2019) first introduce a Graph Laplacian linear operator $\mathcal{L}$ to
transform a generic, non-negative vector
$\boldsymbol{w}\in\mathbb{R}_{+}^{p\left(p-1\right)/2}$ to a Laplacian matrix
$\mathcal{L}\boldsymbol{w}\in\mathbb{R}^{p\times p}$. The linear operator
$\mathcal{L}:\boldsymbol{w}\in\mathbb{R}_{+}^{p\left(p-1\right)/2}\rightarrow\mathcal{L}\boldsymbol{w}\in\mathbb{R}^{p\times
p}$ is formally defined as

$\left(\mathcal{L}\boldsymbol{w}\right)_{ij}=\begin{cases}-w_{i+d_{j}}&i>j,\\\
\left(\mathcal{L}\boldsymbol{w}\right)_{ji}&i<j,\\\ \sum_{i\neq
j}\left(\mathcal{L}\boldsymbol{w}\right)_{ij}&i=j,\\\ \end{cases}$ (26)

where $d_{j}=-j+\frac{j-1}{2}\left(2p-1\right)$. The adjoint operator
$\mathcal{L}^{*}:Y\in\mathbb{R}^{\left(p\times
p\right)}\rightarrow\mathcal{L}^{*}T\in\mathbb{R}^{\frac{p\left(p-1\right)}{2}}$
is derived to satisfy
$\left<\mathcal{L}\boldsymbol{w},Y\right>=\left<\boldsymbol{w},\mathcal{L}^{*}Y\right>$.
While the definition of the two operators might seem cumbersome at the first
glance, their interpretation is fairly straightforward (see Fig, 8).

$w_{1}$, $w_{2}$, $w_{3}$, $w_{4}$, $w_{5}$, $w_{6}$ $\left(\vbox{\hrule
height=9.991pt,depth=9.991pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=9.991pt,depth=9.991pt,width=0.0pt}\right)$$\mathcal{L}$$\mathcal{L}^{-1}$$\sum_{i\in\left\\{1,2,3\right\\}}w_{i}$
$-w_{1}$ $-w_{2}$ $-w_{3}$$\ldots$ $\sum_{i\in\left\\{1,4,5\right\\}}w_{i}$
$-w_{4}$ $-w_{5}$$\ldots$ $\ldots$ $\sum_{i\in\left\\{2,4,6\right\\}}w_{i}$
$-w_{6}$$\ldots$ $\ldots$ $\ldots$
$\sum_{i\in\left\\{3,5,6\right\\}}w_{i}$$\left(\vbox{\hrule
height=49.51805pt,depth=49.51805pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=49.51805pt,depth=49.51805pt,width=0.0pt}\right)$ Figure 8: Given a
Laplacian matrix $Y$, the operator $\mathcal{L}^{-1}$ flattens the upper-
triangular part of $-Y$ into a vector $\boldsymbol{w}$. $\mathcal{L}$ inverts
the process. $y_{11}$ $y_{12}$ $y_{13}$ $y_{14}$ $y_{21}$ $y_{22}$ $y_{23}$
$y_{24}$ $y_{31}$ $y_{32}$ $y_{33}$ $y_{34}$ $y_{41}$ $y_{42}$ $y_{43}$
$y_{44}$ $\left(\vbox{\hrule
height=29.36497pt,depth=29.36497pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=29.36497pt,depth=29.36497pt,width=0.0pt}\right)$$\mathcal{L}^{*}$$y_{11}-y_{21}-y_{12}+y_{22}$
$y_{11}-y_{31}-y_{13}+y_{33}$ $y_{11}-y_{41}-y_{14}+y_{44}$
$y_{22}-y_{32}-y_{23}+y_{33}$ $y_{22}-y_{42}-y_{24}+y_{44}$
$y_{33}-y_{43}-y_{34}+y_{44}$ $\left(\vbox{\hrule
height=46.86436pt,depth=46.86436pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=46.86436pt,depth=46.86436pt,width=0.0pt}\right)$ Figure 9: The adjoint
operator $\mathcal{L}^{*}$ transforms a symmetric matrix in a vector. Above,
an example for a $4\times 4$ matrix.

The Laplacian operator $\mathcal{L}$ allows to reformulate the optimisation
problem in a simpler way. First, by the definition of $\mathcal{L}$, the set
of constraints in Eq. 23 can be expressed as
$\mathcal{S}_{\Theta}=\left\\{\Theta=\mathcal{L}\mathbf{w}|\mathbf{w}\geq
0\right\\}$. Second, if we choose $h\left(\Theta\right)$ to be the
$\mathcal{L}_{1}$-regularisation function, since
$\left(\mathcal{L}\boldsymbol{w}\right)_{ij}<0$ for $i\neq j$ and
$\left(\mathcal{L}\boldsymbol{w}\right)_{ij}>0$ for $i=j$, the regularisation
term $\alpha
h\left(\mathcal{L}\boldsymbol{w}\right)=\alpha\|\mathcal{L}\boldsymbol{w}\|_{1}$
can be written as $\text{tr}\left(\mathcal{L}\boldsymbol{w}H\right)$, where
$H=\alpha\left(2I-\mathbb{1}\right)$, which implies

$\text{tr}\left(\mathcal{L}\boldsymbol{w}S\right)+\alpha
h\left(\mathcal{L}\boldsymbol{w}\right)=\text{tr}\left(\mathcal{L}\boldsymbol{w}K\right),$
(27)

where $K=S+H$. We can now reformulate Eq. 24 as

$\displaystyle\min_{\boldsymbol{w},U}$
$\displaystyle-\log\text{gdet}\left(U\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)U^{T}\right)+\text{tr}\left(\mathcal{L}\boldsymbol{w}K\right)+\frac{\beta}{2}\|\mathcal{L}\boldsymbol{w}-U\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)U^{T}\|_{F}^{2},$
(28) subject to $\displaystyle\boldsymbol{w}>0,U^{T}U=I.$

where $\mathcal{L}\boldsymbol{w}$ is the Laplacian matrix that we would like
to decompose as
$\mathcal{L}\boldsymbol{w}=U\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)U^{T}$,
$\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)\in\mathbb{R}^{p\times p}$
is a diagonal matrix containing $\left\\{\bar{\lambda}_{i}\right\\}$ on its
diagonal , and $U\in\mathbb{R}^{p\times p}$ is an orthogonal matrix. The
constraints on the spectrum of the reconstructed matrix are enforced (softly)
thanks to the spectral penalty term
$\frac{\beta}{2}\|\mathcal{L}\boldsymbol{w}-U\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)U^{T}\|_{F}^{2}$.
It is well known that every Laplacian matrix $\Theta$ will have at least one
eigenvalue equal to zero, since $\Theta\cdot\mathbb{1}=0$ by definition.
Consequently, when solving 28, the first eigenvalue and the corresponding
eigenvector can be dropped from the optimisation formulation. Now
$\bar{\boldsymbol{\lambda}}$ only contains $q=p-1$ non zero eigenvalues in
increasing order, $\left\\{\lambda_{j}\right\\}_{j=2}^{p}$; we can replace the
generalised determinant in 28 with the standard determinant on
$\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)$, and redefine $U$ as
$U\in\mathcal{R}^{p\times q}$, containing the eigenvectors corresponding to
non-zero eigenvalues in the same order. The orthogonality constraint becomes
$U^{T}U=I_{q}$. In Kumar  _et al._ (2019), the authors show how the problem in
28 can be solved with an iterative approach. If we define the vector
$\boldsymbol{c}$,

$\boldsymbol{c}=\left[\mathcal{L}^{*}\left(U\text{Diag}\left(\bar{\boldsymbol{\lambda}}\right)U^{T}-\frac{1}{\beta}K\right)\right],$
(29)

and the function $f\left(\boldsymbol{w}\right)$,

$f\left(\boldsymbol{w}\right)=\frac{1}{2}\|\mathcal{L}\boldsymbol{w}\|^{2}_{F}-\boldsymbol{c}^{T}\boldsymbol{w},$
(30)

at each step $t$, we can update $\boldsymbol{w}$ and $U$, as

$\boldsymbol{w}^{t+1}=\left[\boldsymbol{w}^{t}-\frac{1}{2p}\nabla
f\left(\boldsymbol{w}^{t}\right)\right]^{+},$ (31)
$U^{t+1}=\Lambda\left(\mathcal{L}\boldsymbol{w}\right)[2:p],$ (32)

where $\Lambda\left(\mathcal{L}\boldsymbol{w}\right)$ is the matrix of the
eigenvectors of $\mathcal{L}\boldsymbol{w}$, sorted by the corresponding
eigenvalue. The algorithm can be run until convergence,
$\boldsymbol{w}^{t+1}=\boldsymbol{w}^{t}=\boldsymbol{w}^{*}$, and the vector
$\boldsymbol{w}^{*}$ can be used to reconstruct the Laplacian
$\Theta=\mathcal{L}^{*}\boldsymbol{w}^{*}$, and the corresponding adjacency
matrix. To reconstruct off-diagonal blocks of our Laplacian matrix, we have,
at each iteration steps, only updated the components of $\boldsymbol{w}$
corresponding to off-diagonal blocks, and again run the algorithm until
convergence. While there is no theoretical guarantee that the algorithm will
converge to the optimal solution of the optimisation problem, our results
suggest that this approach is still effective in reconstructing the network.

##  Appendix B Dataset construction

For the purposes of this paper, we accessed three different FactSet products:
Standard Datafeed - Fundamentals V3 - Advanced - Global, Standard Datafeed -
Supply Chain relationship, and APB - Standard Datafeed - Suppply Chain
Shipping Transaction. We parsed information on companies’ fundamentals (sales,
market capitalisation, capital expenditures, industrial sector, and geography)
from the first dataset and used the other two to identify supply-chain
relationship.

#### Fundamentals

The fundamentals dataset is built from the following FactSet files:

  1. 1.

Fundamentals

     * •

ff_basic_eu_v3_full_5315/ff_basic_af_eu.txt

     * •

ff_advanced_eu_v3_full_4524/ff_advanced_af_eu.txt

     * •

ff_basic_ap_v3_full_5276/ff_basic_af_ap.txt

     * •

ff_advanced_der_ap_v3_full_4460/ff_advanced_der_af_ap.txt

     * •

ff_basic_am_v3_full_5258/ff_basic_af_am.txt

     * •

ff_advanced_der_am_v3_full_4484/ff_advanced_der_af_am.txt

  2. 2.

FX Rates

     * •

fx_rates_usd.txt

  3. 3.

Symbology

     * •

sym_hub_v1_full_9915/sym_coverage.txt

     * •

sym_hub_v1_full_9915/sym_entity_sector.txt

     * •

f_sec_hub_v3_full_5299/ff_sec_entity_hist.txt

The Fundamentals files contain the (yearly) information regarding companies
sales, number of employees, and r&d expenses, and a currency column that
states the features’ currency. We can convert all of these features into USD
using the FX Rates table provided by FactSet. The original fundamentals files
are not at the security level, not at the company’s one. To create a dataset
at the company level, FactSet provided us with the following example query,

> Select a.factset_entity_id, c.fsym_id,c.date,c.ff_sales
> from [sym_v1].[sym_sec_entity] a
> join [sym_v1].[sym_coverage] b on a.fsym_id = b.fsym_id
> join [ff_v3].[ff_basic_qf] c on c.fsym_id = b.fsym_regional_id
> where a.factset_entity_id =’05HK0W-E’and a.fsym_id =
> b.fsym_primary_equity_id
> ,

that we "translated" to python. We used
sym_hub_v1_full_9915/sym_entity_sector.txt to assign the correct SIC code to
each of the firms.

#### Supply Chain edgelist

The Supply Chain’s edgelist is built from the following FactSet files:

  1. 1.

Supply Chain

     * •

ent_supply_chain_v1_full_2354/ent_scr_supply_chain.txt

  2. 2.

Shipments

     * •

sc_ship_trans_current_v1_full_1146/sc_ship_trans_curr_1.txt

     * •

sc_ship_trans_current_v1_full_1146/sc_ship_trans_curr_2.txt

     * •

sc_ship_trans_current_v1_full_1146/sc_ship_trans_curr_3.txt

     * •

sc_ship_trans_current_v1_full_1146/sc_ship_trans_curr_4.txt

  3. 3.

Mappings

     * •

ent_entity_advanced_v1_full_6896/factset_entity_structure.csv

     * •

sc_ship_trans_hub_v1_full_1120/sc_ship_parent.txt

The Supply Chain and Shipment files both contain an edgelist (supplier-to-
customer and shipper-to-consignee respectively). The mapping files have two
columns “FACTSET_ENTITY_ID” and “FACTSET_ULT_PARENT_ENTITY_ID”. We assume that
every FACTSET_ENTITY_ID that is not present in the mapping is a ultimate
parent company.

#### Coordinates

The firms’ geographical position was fetched from the following files:

  1. 1.

FactSet’s Addresses

     * •

ent_supply_chain_hub_v1_full_2355/ent_scr_address.txt

     * •

sc_ship_trans_hub_v1_full_1120/sc_ship_address_coord.txt

     * •

sym_hub_v1_full_9915/sym_address.txt’

##  Appendix C Other cleaning strategies

While working on the paper, we tested two other methods to process the
correlation matrix in a way to maximise the gap between the average
correlation along the supply chain and those of the random benchmarks (see
III). After cleaning the market mode, we tried to see whether we could remove
some sector-specific trends out of the time series. For each industrial sector
$\alpha$ we defined the quantity

$s_{\alpha}\left(t\right)=\sum_{i,i\in\alpha}x_{i}\left(t\right),$

where $x_{i}\left(t\right)$ is the growth time series of firm $i$, and the sum
runs on all the firms in sector $alpha$. We assumed that we could write the
time series $x_{i}\left(t\right)$as

$x_{i}\left(t\right)=\xi_{i}\left(t\right)+k_{i}s_{\alpha}\left(t\right),$

We estimated the coefficient $k_{i}$ as the correlation between $x_{i}$ and
$s_{\alpha}$, and cleaned the time series by computing the difference

$\xi_{i}(t)=x_{i}\left(t\right)-\hat{k}_{i}s_{\alpha}(t),$

where $\hat{k}_{i}$ is the estimated value for $k_{i}$.

We also investigated if more signal could be extracted by considering lags
between firms time series. We defined the lagged correlation matrix
$C\left(\tau\right)$, defined as

$C\left(\tau\right)=\mathbb{E}_{t}\left[x_{i}(t)x_{j}(t+\tau)\right],$

and its symmetrised version $C^{\prime}\left(\tau\right)$ as

$C^{\prime}(1)=\frac{1}{2}\left[C(1)+C(-1)\right].$

We then computed a linear combination
$\left[C^{\prime}(0)+C^{\prime}(1)\right]$, and computed the average value of
this matrix over the supply chain and the random benchmarks.

None of the two approaches improved significantly the outcomes we discussed.
However, we can’t exclude that a more thorough investigation of these
techniques, their combination, and the analysis of other time series (e.g.,
firms’ market returns) could improve the results of this paper.