# Learning production functions for supply chains
with graph neural networks

Serina Chang  Stanford University, Department of Computer Science  Zhiyin Lin
Stanford University, Department of Computer Science  Benjamin Yan  Stanford
University, Department of Computer Science  Swapnil Bembde  Hitachi America,
Ltd.  Qi Xiu  Hitachi America, Ltd.  Chi Heem Wong  Stanford University,
Department of Computer Science  Hitachi America, Ltd.
Yu Qin  Hitachi America, Ltd.  Tulane University  Frank Kloster  Hitachi
America, Ltd.  Alex Luo  Hitachi America, Ltd.  Raj Palleti  Stanford
University, Department of Computer Science  Hitachi America, Ltd.  Jure
Leskovec  Stanford University, Department of Computer Science

###### Abstract

The global economy relies on the flow of goods over supply chain networks,
with nodes as firms and edges as transactions between firms. While we may
observe these external transactions, they are governed by unseen production
functions, which determine how firms internally transform the input products
they receive into output products that they sell. In this setting, it can be
extremely valuable to infer these production functions, to better understand
and improve supply chains, and to forecast future transactions more
accurately. However, existing graph neural networks (GNNs) cannot capture
these hidden relationships between nodes’ inputs and outputs. Here, we
introduce a new class of models for this setting, by combining temporal GNNs
with a novel inventory module, which learns production functions via attention
weights and a special loss function. We evaluate our models extensively on
real supply chains data, along with data generated from our new open-source
simulator, SupplySim. Our models successfully infer production functions, with
a 6-50% improvement over baselines, and forecast future transactions on real
and synthetic data, outperforming baselines by 11-62%.

##  1 Introduction

Supply chains form the backbone of the global economy and disruptions can have
enormous consequences, costing trillions of dollars (Baumgartner et al.,,
2020) and risking national security (The White House,, 2021). Thus, modeling
supply chains and predicting how they will evolve, especially under shocks,
are essential. To model supply chains, we can represent them as temporal (or
dynamic) graphs, where nodes are firms and edges represent time-varying
transactions between firms. A unique feature of these graphs is that they are
governed by underlying production functions: firms receive input products
(e.g., wheels) from other firms, internally transform those inputs into
outputs (e.g., cars) via production functions, then sell those output products
to other firms or consumers (Carvalho and Tahbaz-Salehi,, 2019). These
production functions dictate which firms are connected to each other, as well
as the timing of transactions across the network. For example, if a firm
experiences a shortage in an input, its production of outputs that rely on
that input will be disrupted, but it can continue to produce other outputs for
some time. Conversely, if a firm needs to produce more of a particular output,
due to a change in demand, it will increase its orders of the inputs required
to produce that output, but orders for other inputs may remain unchanged.

Despite growing interest in developing graph neural networks (GNNs) for
temporal graphs (Huang et al.,, 2023; Longa et al.,, 2023), temporal graphs
governed by production functions—which we term temporal production graphs
(TPGs) in this work—have not yet been explored. Beyond supply chains, examples
also appear in biological domains, such as enzymes producing metabolites
within metabolic pathways, or in organizations, where teams rely on inputs
from other teams to produce their outputs (e.g., engineers interfacing with
product designers). TPGs introduce new challenges compared to other temporal
graphs, such as the need to infer the production functions. Inferring a TPG’s
production functions can help us better understand the graph’s mechanisms and
intervene more effectively, such as identifying bottlenecks in supply chains
and proposing plausible solutions to improve efficiency and sustainability
(Aigner and Chu,, 1968; Coelli et al.,, 2005). Furthermore, even standard
tasks, like link prediction, may be complicated under TPGs. For example, while
standard temporal GNNs might capture disruptions generally propagating across
connected firms in supply chains, they will miss the specific connections
between each firm’s inputs and outputs and, under a shortage of inputs, not be
able to precisely predict which outputs will be affected.

To fill this gap, here we introduce a new class of GNNs designed for TPGs,
focusing on two objectives: (1) learning the graph’s internal production
functions, (2) predicting its future edges. Prior temporal GNNs have focused
primarily on the second objective, but are not designed for the first. Our
models support both objectives, by combining temporal GNNs with a novel
inventory module, which learns production functions by explicitly representing
each firm’s time-varying inventory and updating it based on attention weights
that map external supply to internal consumption. Our module can be combined
with any GNN; to demonstrate this, here we combine it with two popular
temporal GNNs, Temporal Graph Network (TGN) (Rossi et al.,, 2020) and
GraphMixer (Cong et al.,, 2023), which we also extend in important ways. In
this work, we have rare access to real-world transaction-level supply chains
data, which we use to evalaute our models. We also build a new open-source
simulator, SupplySim, which generates realistic supply chain data that matches
the real data on key characteristics. Our simulator enables us to share data,
encourage further research development on TPGs, and test models under
controlled and varied settings (e.g., supply shocks, missing data). In
summary, our contributions are:

  1. 1.

Problem: a new graph ML problem setting, temporal production graphs (TPGs),
where each node’s edges are related via unobserved production functions,

  2. 2.

Models: a new class of models for TPGs, which combine temporal GNNs with a
novel inventory module to jointly learn production functions and forecast
future edges,

  3. 3.

Data: an open-source simulator, SupplySim, which generates realistic supply
chain data and enables model testing under a wide variety of settings,

  4. 4.

Results: experiments on two real-world supply chain datasets and synthetic
datasets from SupplySim, showing that our models effectively learn production
functions (outperforming baselines by 6-50%) and forecast future edges
(outperforming baselines by 11-62%).

While we focus on supply chains in this work, our models could be applied or
easily extended to other TPGs. Our work pushes forward the frontier of graph
ML and its application to real-world problems, by highlighting a new type of
real-world graph and introducing new methods to handle them.111Our code to run
all experiments, simulator SupplySim, and synthetic datasets are available at
https://github.com/snap-stanford/supply-chains/tree/tgb.

##  2 Learning from temporal production graphs

###  2.1 Problem definition

We define a new graph ML setting called temporal production graphs (TPGs),
illustrated in Figure 1a. TPGs are directed graphs with time-varying edges and
potentially a time-varying set of nodes.

Figure 1: (a) Illustration of our problem setting: we observe time-varying
transactions between firms and do not observe production functions within
firms. Our goals are to learn the production functions and predict future
transactions. (b) Example of our model architecture, combining our inventory
module with our extended version of TGN, SC-TGN.

What differentiates TPGs from other temporal graphs is that, in a TPG, each
node’s in-edges represent inputs to some internal production function, which
are transformed into outputs, represented by the node’s out-edges. In this
work, we focus on supply chain networks, as a canonical example of TPGs. We
are given a set of transactions $\mathcal{T}$ between firms where, for each
transaction, we know its timestamp, supplier firm, buyer firm, product sold,
and amount sold. We represent this data as a heterogeneous temporal graph,
$G_{\textrm{txns}}=\\{\mathcal{N},\mathcal{E}\\}$, where the nodes
$\mathcal{N}$ consist of $n$ firm nodes and $m$ product nodes and the edges
are $\mathcal{E}:=\\{e(s,b,p,t)\\}$, where $e(s,b,p,t)$ is a hyperedge between
supplier firm $s$, buyer firm $b$, and product $p$, representing a transaction
between them at time $t$.

In this setting, production functions define how firms internally transform
the products that they buy into products that they supply. Specifically, the
function $\mathcal{F}_{p}:\mathbb{R}_{+}\rightarrow\mathbb{R}^{m}_{+}$ for
product $p$ defines how much of each product is necessary to make $k$ amount
of $p$ (e.g., a car requires four wheels). The set of production functions
define a production graph, $G_{\textrm{prod}}$, which is a directed acyclic
graph where there is an edge from products $p_{1}$ to $p_{2}$ if $p_{1}$ is
required to make $p_{2}$.222For simplicity, we assume in this work that there
is one way to make each product, but future work may consider extensions, such
as firm-specific product graphs or products that are substitutes for each
other. Given the set of transactions $\mathcal{T}$, and the resulting temporal
graph $G_{\textrm{txns}}$, our goals are two-fold: (1) to infer
$G_{\textrm{prod}}$, which is entirely unobserved, (2) to predict future
transactions, i.e., future hyperedges in $G_{\textrm{txns}}$.

###  2.2 Model architecture

To learn from TPGs, we introduce a new class of models that combine temporal
GNNs with a novel inventory module to jointly learn production functions and
predict future edges. First, we describe our inventory module, which can
either operate as a stand-alone model or be attached to any GNN. Second, we
describe extended versions of two popular temporal GNNs, Temporal Graph
Network (Rossi et al.,, 2020) and GraphMixer (Cong et al.,, 2023), which we
refer to as SC-TGN (shown in Figure 1b) and SC-GraphMixer, respectively, with
SC standing for supply chain.

####  2.2.1 Inventory module

The basic idea of our inventory module is that it explicitly represents each
firm’s inventory of products, and it adds and subtracts from the inventory
based on products bought and products consumed, respectively. Let
$\mathbf{x}_{i}^{(t)}\in\mathbb{R}_{+}^{m}$ represent firm $i$’s inventory at
time $t$. We compute the total amount of product $p$ bought by firm $i$ at
time $t$ as

$\displaystyle\textrm{buy}(i,p,t)$
$\displaystyle=\sum_{e(s,i,p,t)\in\mathcal{E}}\textrm{amt}(s,i,p,t),$ (1)

where $\textrm{amt}(s,i,p,t)$ represents the amount of the product in the
transaction indexed by $s,i,p,t$. The amount bought per product is computed
directly from the data, since we observe the products that a firm buys through
its transactions. On the other hand, we cannot observe the products that a
firm consumes from their inventory; only the finished products that they
supply to others. So, we need to learn the function mapping from products
externally supplied to products internally consumed from the inventory. We
estimate the total amount of product $p$ consumed by firm $i$ at time $t$ as

$\displaystyle\textrm{consume}(i,p,t)$
$\displaystyle=\sum_{e(i,b,p_{s},t)\in\mathcal{E}}\alpha_{p_{s}p}\cdot\textrm{amt}(i,b,p_{s},t),$
(2)

where $\alpha_{p_{s}p}\in\mathbb{R}_{+}$ is a learned attention weight
representing how much of product $p$ is needed to make one unit of the product
supplied, $p_{s}$. In other words, each product attends to its parts. Finally,
let $\mathbf{b}_{i}^{(t)}$ represent the vector over all products of amount
bought by firm $i$ at time $t$, and let $\mathbf{c}_{i}^{(t)}$ be defined
analogously for consumption. Then the firm’s updated inventory is

$\displaystyle\mathbf{x}_{i}^{(t+1)}$
$\displaystyle=\max(0,\mathbf{x}_{i}^{(t)}+\mathbf{b}_{i}^{(t)}-\mathbf{c}_{i}^{(t)}).$
(3)

We take the elementwise max with 0 to ensure that the inventory stays non-
negative, but we also penalize whenever consumption exceeds the current
inventory amounts (see Eq. 6 below).

##### Attention weights.

If the inventory module is standalone, then we directly learn the pairwise
weights $\alpha_{p_{1}p_{2}}$ for all product pairs $p_{1},p_{2}$. If the
inventory module has access to product embeddings (e.g., from a GNN), we can
use product embedding $\mathbf{z}_{p}\in\mathbb{R}^{d}$ to inform the
attention weights as

$\displaystyle\alpha_{p_{1}p_{2}}$
$\displaystyle=\textrm{ReLU}(\mathbf{z}_{p_{1}}\mathbf{W}_{\textrm{att}}\mathbf{z}_{p_{2}}+\nu_{p_{1}p_{2}}),$
(4)

where $\mathbf{W}_{\textrm{att}}\in\mathbb{R}^{d\times d}$ and
$\nu_{p_{1}p_{2}}\in\mathbb{R}$ are learned parameters, and we apply ReLU to
ensure that the attention weights are non-negative. Compared to directly
learning the attention weights, here we treat
$\mathbf{z}_{p_{1}}\mathbf{W}_{\textrm{att}}\mathbf{z}_{p_{2}}$ as the base
rate and $\nu_{p_{1}p_{2}}$ as adjustments, which we encourage to be small in
magnitude with $L_{2}$ regularization (6). By using the embeddings to form the
base rate, instead of learning each pair independently, we can share
information across product pairs, which is especially useful given sparse
real-world data where we rarely observe most pairs.

##### Inventory loss.

To train the inventory module, we introduce a special loss function. For a
given firm $i$ at time $t$, its inventory loss is

$\displaystyle\ell_{\textrm{inv}}(i,t)$
$\displaystyle=\lambda_{\textrm{debt}}\sum_{p\in[m]}\max(0,\textrm{consume}(i,p,t)-\mathbf{x}_{i}^{(t)}[p])-\lambda_{\textrm{cons}}\sum_{p\in[m]}\textrm{consume}(i,p,t).$
(5)

That is, we penalize inventory debt, i.e., whenever consumption exceeds the
current inventory, but otherwise reward consumption. We penalize inventory
debt since firms should not be able to consume products that they never
received. Furthermore, penalizing inventory debt results in sparse attention
weights, since for most firms, we do not observe it buying most products, so
for any of those products that it does not buy, it would prefer to never
consume those products since it would immediately go into inventory debt if
so. On the other hand, we need the consumption reward in order to prevent
trivial solutions. Without it, the model could learn
$\alpha_{p_{1},p_{2}}\approx 0$ for all product pairs, and it would never
experience inventory debt. We use hyperparameters $\lambda_{\textrm{debt}}$
and $\lambda_{\textrm{cons}}$ to control the relative weight between
penalizing inventory debt and rewarding consumption, and in practice, we find
that choosing $\lambda_{\textrm{debt}}$ around $25\%$ larger than
$\lambda_{\textrm{cons}}$ works well (Table 4). All together, the inventory
loss is

$\displaystyle\ell_{\textrm{inv}}(t)$
$\displaystyle=\frac{1}{n}\sum_{i\in[n]}\ell_{\textrm{inv}}(i,t)+\lambda_{L_{2}}\sqrt{\sum_{p_{1},p_{2}\in[m]}\nu_{p_{1}p_{2}}^{2}}.$
(6)

####  2.2.2 SC-TGN

The first GNN we explore is Temporal Graph Network (TGN) (Rossi et al.,,
2020), which is one of the most established GNNs for dynamic link prediction,
outperforming other models in the Temporal Graph Benchmark (Huang et al.,,
2023). In our work, we have extended TGN to SC-TGN, by enabling it to (1)
perform message passing over hypergraphs, since we represent each transaction
as an edge between three nodes, (2) predict edge weights (i.e., transaction
amounts) in addition to edge existence. We also modified TGN in other ways
that improved performance; we document these changes in Appendix A.1. In SC-
TGN, each node $i$ has a time-varying memory $\mathbf{m}_{i}^{(t)}$. Each
transaction $e(s,b,p,t)$ sends three messages: one each to the supplier $s$,
buyer $b$, and product $p$. At the end of each timestep, each node aggregates
the messages it received and updates its memory, using a recurrent neural
network (RNN). Then, to produce node embedding $\mathbf{z}_{i}^{(t)}$, we
apply a GNN to the node memories, so that nodes can also learn from their
neighbors’ memories (this is useful for preventing staleness, if a node has
not had a transaction in a while).

####  2.2.3 SC-GraphMixer

We also explore GraphMixer (Cong et al.,, 2023), a recent model that showed
that temporal GNNs do not always need complicated architectures, such as RNNs
or self-attention (both of which are used by TGN), and strong performance can
sometimes be achieved with simpler models that only rely on multi-layer
perceptrons (MLPs). We similarly extend GraphMixer to SC-GraphMixer, so that
it can handle hypergraphs and predict edge weight in addition to edge
existence (Appendix A.2). In SC-GraphMixer, each node’s embedding
$\mathbf{z}_{i}^{(t)}$ is a concatenation of its node encoding and link
encoding. The node encoding is simply a sum of the node’s features and mean-
pooling over its 1-hop neighbors’ features. The link encoding summarizes the
recent edges that the node participated in, by applying an MLP to the
concatenated time encodings and features of the recent edges.

####  2.2.4 Decoder

Both models, SC-TGN and SC-GraphMixer, use the same decoder architecture, and
we use the same architecture (although separate decoders) for predicting edge
existence and edge weight. We model these as a two-step process: first,
predicting whether an edge exists; second, conditioned on the edge existing,
predicting its weight. The decoder is a two-layer MLP over the concatenation
of the supplier firm’s, buyer firm’s, and product’s embeddings, producing
$\hat{y}\in\mathbb{R}$:

$\displaystyle\hat{y}(s,b,p,t)$
$\displaystyle=\textrm{DEC}(\mathbf{z}_{s}^{(t)},\mathbf{z}_{b}^{(t)},\mathbf{z}_{p}^{(t)})=\textrm{MLP}([\mathbf{z}_{s}^{(t)}|\mathbf{z}_{b}^{(t)}|\mathbf{z}_{p}^{(t)}]).$
(7)

Real values are natural for both prediction tasks, since for edge existence,
we apply a softmax over alternatives and use cross-entropy loss to evaluate
the probability of the positive transaction compared to negative samples, and
for edge weight, we apply log-scaling and standard scaling to the transaction
amounts so negative predicted amounts are valid.

When the inventory module is attached, we may also allow it to inform edge
prediction. To help with predicting edge existence, the inventory module
penalizes impossible transactions, i.e., some $(s,b,p,t)$ where the supplier
$s$ does not have enough parts to make product $p$:

$\displaystyle pen(s,b,p,t)$
$\displaystyle=-\sum_{p^{\prime}\in[m]}\max(0,\alpha_{pp^{\prime}}-\mathbf{x}_{s}^{(t)}[p^{\prime}]).$
(8)

Then, our model’s new prediction for edge existence becomes the sum of the
original $\hat{y}$ in (7) and the penalty (8). To help with predicting edge
weight, i.e., transaction amount, the inventory module computes a cap on the
maximum amount of product $p$ that supplier $s$ could produce at this time:

$\displaystyle cap(s,b,p,t)$
$\displaystyle=\min_{p^{\prime}\in[m];\alpha_{pp^{\prime}}>0}\\{\frac{\mathbf{x}_{s}^{(t)}[p^{\prime}]}{\alpha_{pp^{\prime}}}\\}.$
(9)

Then, our new predicted amount is the minimum of the original predicted amount
(7) and the cap (after applying the same log-scaling and standard scaling to
the cap that we did to all amounts). Thus, as shown in Figure 1b, the
inventory module and GNN help each other, with the GNN’s embeddings informing
the inventory module’s attention weights, and the inventory module’s penalties
and caps affecting future edge prediction.

###  2.3 Model training and evaluation

##### Learning production functions.

Our model does not have access to any production functions during training,
but the goal is for the model to learn production functions via the inventory
module and its specialized loss. For each product $p$, we have the set of its
attention weights
$\\{\alpha_{pp_{1}},\alpha_{pp_{2}},\cdots,\alpha_{pp_{m}}\\}$, and we can
compute the ranking over all products from highest to lowest weight. We
compare this ranking to the product’s ground-truth parts, i.e., the parents of
$p$ in $G_{\textrm{prod}}$, and use average precision (17) to quantify
performance. Then, we compute the mean average precision (MAP) over all
products for which we have ground-truth parts. We have ground-truth parts for
all products in our synthetic data (Section 3.2) and for scanning electron
microscopes in real-world data (Section 3.1).

##### Edge existence and negative sampling.

We perform negative sampling such that each positive transaction $e(s,b,p,t)$
is paired with a set of negative transactions that did not actually occur at
time $t$. Following Temporal Graph Benchmark (Huang et al.,, 2023), we sample
two types of hard negatives: first, randomly perturbing one of the three
nodes; second, sampling a historical negative, meaning a transaction that
appeared in training but not at time $t$. Compared to a uniform random
negative, these negatives are far harder to distinguish from the positive
transaction. For the perturbations, they match the positive transaction on two
of the three nodes, requiring the model to make a distinction based on the
third node. The historical negatives are even harder: this transaction has
been observed before, but it does not occur at this timestep, requiring the
model to learn time-varying properties of the graph. For each positive
transaction, we sample 9 perturbation negatives and 9 historical negatives. To
evaluate performance, we use mean reciprocal rank (MRR), which evaluates the
rank of the positive transaction among the negatives (19). However, MRR is
non-differentiable, so during training, we use softmax cross-entropy as a
strong proxy loss for MRR (Bruch et al.,, 2019).

##### Edge weight.

Since we model edge prediction as a two-stage process, we only predict edge
weight (i.e., transaction amount) conditioned on the edge existing. Thus,
training and test edge weight prediction is simple: we only consider the
positive transactions and we compare the model’s predicted amount to the true
amount using root mean squared error (RMSE) (20). Unlike MRR, RMSE is
differentiable, so we also use it in the model loss during training.

##  3 Supply chains data

To evaluate our models, we collected two real-world datasets of supply chain
transactions, and built a simulator, SupplySim, which can generate realistic
supply chain data under varied settings. While we cannot release the real-
world datasets, due to their proprietary nature, we carefully document our
data sources and collection processes. Furthermore, we release SupplySim
(which matches the real data on key characteristics) and the three synthetic
datasets that we evaluate our models on in this work. We describe both data
types here, with details in Appendix B and data statistics in Table 3.

###  3.1 Real-world supply chains data

We acquired transactions data from TradeSparq, a third-party data provider
which aggregates data from authorized government sources (e.g., customs
departments) across 60+ countries. Their data sources include bills of lading,
receipts of reported transactions, and custom declarations. Each product being
transacted is described with a Harmonized System (HS) code, where HS6 is an
internationally recognized six digit number that describes products in terms
of their chapters (first two digits), headings (next two digits), and
subheadings (last two digits). For example, the HS code 840731 refers to the
chapter 84 (chapter: Machinery), 07 (heading: Spark-ignition Reciprocating
Piston Engines), and 31 (sub-heading: engines type not over 50CC cylinder
capacity). Along with HS codes, the transactions data also includes the
company name and ID of the supplier, the company name and ID of the buyer,
timestamp of the transaction, cost in USD, and so on. Using the TradeSparq
data (which we could only access via API calls, with costs per call), we
constructed two supply chains datasets. For both datasets, we began by
selecting specific HS codes that we were interested in, identifying their
prominent suppliers, then using breadth-first-search to collect their
neighbors upstream and downstream in the supply chain (see flowchart in Figure
3). This collection process allowed us to capture all interactions between
firms within multiple tiers of the supply chain.

Automotive dataset (Tesla). This dataset focuses on electric vehicles (EV) and
EV parts supplied by Tesla. We identified Tesla EV makers, their direct
suppliers and buyers, and their suppliers’ suppliers, and included all
transactions between these firms from January 1, 2019, to December 31, 2022.

Industrial equipment dataset (IED). This dataset focuses on microscopes along
with other specialized analytical and inspection equipment and their
manufacturers. We included makers of these products and their direct suppliers
and buyers, and included all transactions between these firms in 2023.
Compared to the Tesla dataset, this dataset is broader (i.e., more firms per
layer) but shorter in time and depth (i.e., fewer layers). For this dataset,
we have also estimated the ground-truth parts of microscopes in terms of their
HS codes, which we use to test our inventory module’s ability to infer
production functions from transactions.

###  3.2 Supply chains simulator, SupplySim

The TradeSparq data offers a rare opportunity to test models on real
transaction-level data, but it also has shortcomings: the transactions are
incomplete, most production functions are missing, and we cannot release the
data. Thus, we design a simulator, SupplySim, that addresses these
shortcomings and enables us to test the model under controlled settings (e.g.,
complete vs. incomplete transactions, stable supply vs. shocks). There are
three stages to the simulator: (1) constructing the production graph,
$G_{\textrm{prod}}$, (2) assigning products to suppliers and assigning
suppliers to buyers, (3) generating time-varying transactions between firms.
We ensure that our generated graphs match real-world supply chain networks on
known characteristics, such as power law degree distributions.

##### Constructing the production graph, $G_{\textrm{prod}}$.

Given the number of products, $m$, we partition the products into tiers and
sample a 2-dimensional position for each product from Uniform($0,1$). Then,
for each tier, we assign the products in the tier to parts from the previous
tier, with probability proportional to the inverse distance between their
positions. For each part-product pair, we also sample $u_{io}$, the number of
units of the part $p_{i}$ needed to make one unit of the product $p_{o}$. The
tier structure imitates tiers in real-world supply chains; for example,
products at the first tier are raw materials, e.g., minerals, and products at
the final tier are consumer products. The products’ positions capture the
product type, e.g., its industry, and assigning parts based on positions
reflects how parts should be similar to their products and naturally results
in commonly co-occurring parts.

##### Constructing supplier-buyer relations.

Given the number of firms, $n$, we also sample a 2-dimensional position from
Uniform($0,1$) for each firm, and we restrict the firm to a subset of two
consecutive tiers, meaning it can only produce products in those tiers. Then,
for each product, we select its suppliers from the set of firms that are
allowed to produce that product, again with probability proportional to
inverse distance. Now, each firm has a set of products that it is supplying,
which means, based on $G_{\textrm{prod}}$, we know which input parts it needs
to buy. For each pair $(b,p)$, where firm $b$ needs to buy product $p$, we
assign it to a supplier of $p$ with probability proportional to the number of
buyers that the supplier already has. This assignment mechanism, known as
preferential attachment (Newman,, 2001), yields power law degree
distributions, which are known to appear in real-world supply chain networks
(Fujiwara and Aoyama,, 2010; Zhao et al.,, 2019).

##### Generating transactions.

We generate transactions based on the ARIO model, an agent-based model widely
used in economics to simulate propagation over supply chains (Hallegatte,,
2008; Inoue and Todo,, 2019; Guan et al.,, 2020). At each timestep of the
simulation, each firm completes as many of its incomplete orders as it can
until it runs out of inventory. Each time it completes an order for $k$ amount
of output product $p_{o}$, it subtracts $k\cdot u_{io}$ of input product
$p_{i}$ from its inventory (i.e., we assume the production function
$\mathcal{F}_{p_{o}}$ is linear). At the end of the timestep, the firm places
orders to its suppliers, based on what it needs to complete its remaining
orders. Finally, the reported transactions are the completed orders in each
timestep. The simulator also keeps track of the exogenous supply for products
in the first tier, which do not require parts, and of exogenous demand for
products in the final tier, which are only bought by consumers and not by
other firms. By manipulating the exogenous supply, we can model shocks in the
supply chain.

##  4 Experiments

We run experiments on the two real-world supply chain datasets and three
synthetic datasets from SupplySim: a standard setting with high exogenous
supply (“SS-std”), a setting with shocks to exogenous supply (“SS-shocks”),
and a setting with missing transactions (“SS-missing”), where we sampled 20%
of firms uniformly at random and dropped all transactions where they were
buyer or supplier, reflecting firms that are missing from real-world supply
chains data. In all experiments, we order the transactions chronologically and
split them into train (the first 70%), validation (the following 15%), and
test (the last 15%). Here we describe our results on these datasets, with
additional details about the experimental set-up and supplemental results in
Appendix C.

###  4.1 Learning production functions

|  SS-std |  SS-shocks |  SS-missing |  IED
---|---|---|---|---
Random baseline |  0.124 (0.009) |  0.124 (0.009) |  0.124 (0.009) |  0.060 (0.002)
Temporal correlations |  0.745 |  0.653 |  0.706 |  0.128
PMI |  0.602 |  0.602 |  0.606 |  0.175
node2vec |  0.280 |  0.280 |  0.287 |  0.127
Inventory module (direct) |  0.771 (0.005) |  0.770 (0.006) |  0.744 (0.006) |  0.143 (0.004)
Inventory module (emb) |  0.790 (0.005) |  0.778 (0.011) |  0.755 (0.007) |  0.262 (0.005)
Table 1: Results for production learning, evaluated with mean average
precision (MAP $\uparrow$). For the models with randomness, we report mean and
standard deviation (in parentheses) over 10 seeds.

We try three baselines, which each compute scores for output product $p_{o}$
and potential input $p_{i}$:

1) Temporal correlations: we expect that inputs and outputs are temporally
correlated, so this method computes the maximum correlation, with possible
lags, between the buying timeseries of $p_{i}$ and supplying timeseries of
$p_{o}$, averaged over all firms that buy $p_{i}$ and supply $p_{o}$.

2) Pointwise Mutual Information (PMI): we expect that input-output pairs
appear with greater frequency in the firm-product graph,333The firm-product
graph is a static bipartite graph of firm and product nodes, where an edge
between a firm and product indicates that the firm buys or supplies the
product. so this method computes the probability that a firm buys $p_{i}$ and
supplies $p_{o}$, divided by the product of their individual probabilities
(26).

3) node2vec: we expect that inputs are close to outputs, so this method
computes the cosine similarity of $p_{i}$ and $p_{o}$’s node2vec embeddings,
learned over the firm-product graph.

The three baselines capture the usefulness of temporal information, 1-hop
neighbors in the graph, and the entire graph, respectively. In contrast, our
inventory module captures both temporal and structural information. We try two
versions of the inventory module, one that learns the attention weights
directly and one that uses product embeddings, as described in (4). We also
report a random baseline, which produces uniform random rankings of parts for
each product.

We summarize the production learning results in Table 1. We find that the
inventory module significantly outperforms the baselines, with especially
large margins on the real-world data (IED) and in the synthetic data, when
there are shocks in supply. In the standard synthetic data, when supply is
plentiful, temporal correlations are a strong predictor since firms place
orders for inputs, receive them promptly, then supply their own outputs
shortly after. However, once there are shocks, the firm will receive inputs at
different times, due to delays, and since they cannot produce their outputs
until all inputs have arrived, the correlation in time of buying inputs and
supplying outputs is seriously worsened. On the other hand, the inventory
module is robust to such delays, since it does not rely on similarity in
timeseries; simply that an input must go into the inventory before the output
comes out. The inventory module is also remarkably robust to missing data: the
MAP only drops by 3.5 points (4.4%) when we drop 20% of firms in the synthetic
data. We also see that using product embeddings, instead of learning the
attention weights directly, consistently helps the inventory module.

Figure 2 visualizes our results, showing that the inventory module effectively
learns the true production functions. The block structure visible in all three
subfigures reflects how our products were partitioned into tiers (products 0-4
in tier 0, 5-14 in tier 1, 15-24 in tier 2, and so on) and products in one
tier served as viable parts for the next tier. The inventory module recovers
the tier structure much more precisely, while temporal correlations show
positive weights in larger blocks. The inventory module also learns attention
weights of similar magnitude as the true production functions, while the
temporal correlations are constrained between 0 and 1, and the other baselines
are similarly not comparable to the true production functions on magnitude,
only ranking. While only ranking matters for MAP, magnitude is essential if
using the inventory module to inform edge prediction with penalties (8)-(9).
Finally, in Figure 6, we show that our inventory module’s loss function (6) is
well-correlated with MAP, which is why it can effectively learn production
functions without observing any of them.

Figure 2: True production functions (left), predictions from inventory module
(middle), predictions from temporal correlations (right), trained on SS-std.

###  4.2 Predicting future edges

For edge existence, where we seek to predict whether future transaction
$(s,b,p,t)$ exists, we compare our models to the following baselines:

1) Edgebank: the “binary” version predicts 1 if $(s,b,p,*)$ appeared before in
the train set; 0 otherwise. The “count” version predicts the number of times
that $(s,b,p,*)$ appeared in the train set.

2) Static: learns a static vector to represent each node (no GNN).

3) Graph transformer: learns a static embedding to represent each node, using
the UniMP model (Shi et al.,, 2021) which utilizes a graph transformer
architecture and multi-head attention.

While Edgebank simply memorizes the train set, it serves as a strong baseline,
as the Temporal Graph Benchmark (Huang et al.,, 2023) found that even the
binary version outperformed some GNNs on dynamic link prediction. The latter
two baselines test two types of static node embeddings, which allow us to
isolate the benefit of dynamic node embeddings. For each of our models (SC-TGN
and SC-GraphMixer), we try them alone and with the inventory module (+inv).
For edge weight, we test the same set of methods, except for Edgebank, where
we now predict the average amount for $(s,b,p,*)$ in the train set, or 0
(which is the mean amount, after standard scaling) if the triplet never
appeared in train.

|  SS-std |  SS-shocks |  SS-missing |  Tesla |  IED
---|---|---|---|---|---
Edgebank (binary) |  0.174 |  0.173 |  0.175 |  0.131 |  0.164
Edgebank (count) |  0.441 |  0.415 |  0.445 |  0.189 |  0.335
Static |  0.439 (0.001) |  0.392 (0.002) |  0.442 (0.001) |  0.321 (0.001) |  0.358 (0.001)
Graph transformer |  0.431 (0.003) |  0.396 (0.024) |  0.428 (0.003) |  0.507 (0.020) |  0.613 (0.045)
SC-TGN |  0.522 (0.003) |  0.449 (0.004) |  0.494 (0.004) |  0.820 (0.007) |  0.842 (0.004)
SC-TGN+inv |  0.540 (0.003) |  0.461 (0.009) |  0.494 (0.004) |  0.818 (0.004) |  0.841 (0.008)
SC-GraphMixer |  0.453 (0.005) |  0.426 (0.004) |  0.446 (0.003) |  0.690 (0.027) |  0.791 (0.009)
SC-GraphMixer+inv |  0.497 (0.004) |  0.448 (0.004) |  0.446 (0.002) |  0.681 (0.014) |  0.791 (0.008)
Edgebank (avg) |  0.341 |  0.387 |  0.349 |  1.148 |  0.489
Static |  0.343 (0.008) |  0.425 (0.019) |  0.374 (0.027) |  1.011 (0.007) |  0.504 (0.018)
Graph transformer |  0.340 (0.005) |  0.398 (0.025) |  0.361 (0.016) |  0.885 (0.024) |  0.425 (0.008)
SC-TGN |  0.303 (0.003) |  0.359 (0.007) |  0.313 (0.002) |  0.796 (0.012) |  0.428 (0.011)
SC-TGN+inv |  0.312 (0.003) |  0.370 (0.009) |  0.312 (0.002) |  0.801 (0.015) |  0.422 (0.011)
SC-GraphMixer |  0.318 (0.003) |  0.384 (0.005) |  0.330 (0.005) |  0.774 (0.077) |  0.457 (0.008)
SC-GraphMixer+inv |  0.320 (0.004) |  0.378 (0.005) |  0.328 (0.003) |  0.767 (0.054) |  0.454 (0.012)
Table 2: Results for edge prediction. Top 8 rows are edge existence, evaluated
with mean reciprocal rank (MRR $\uparrow$). Bottom 7 rows are edge weight
(i.e., transaction amount), evaluated with root mean squared error (RMSE
$\downarrow$). We report mean and standard deviation (in parentheses) over 10
seeds.

We summarize the edge prediction results in Table 2. First, we find that our
models, SC-TGN and SC-GraphMixer, outperform the baselines on both tasks and
all five datasets, by 11-62% on edge existence and 1-13% on edge weight.
Second, we find that SC-TGN consistently outperforms SC-GraphMixer; we
hypothesize this is because of SC-TGN’s sophisticated updating of node
memories, while SC-GraphMixer only captures changes over time through its link
encoder, which encodes the features and timestamps of the node’s recent edges,
and its node encoder relies on the node’s static features (Appendix A.2).
Supporting this hypothesis is the fact that the static methods perform poorly,
demonstrating the need for TPGs to be represented with dynamic node
embeddings. Since we find that SC-TGN outperforms SC-GraphMixer, we also
compared SC-TGN to ablated versions: the original TGN, to test the value of
our extensions (which we document in Appendix A.1), and to SC-TGN where the
node memory is directly used as the embedding, instead of applying a GNN
(UniMP in our case) to the memories. We find that the full SC-TGN outperforms
both ablations substantially, by 34-45% and 53-100%, respectively (Table 5).

We also find in the synthetic data experiments, the setting with shocks is the
hardest, resulting in lower MRRs and higher RMSEs across models. This
demonstrates how shocks complicate prediction on TPGs, since shocks both delay
firms’ abilities to complete their orders and limit them to producing smaller
amounts, thus affecting both edge existence and edge weight. However, adding
the inventory module significantly improves MRRs for both SC-TGN and SC-
GraphMixer, on standard synthetic data and with shocks. In general, even
though adding the inventory module introduces an additional loss so that we
can learn production functions (16), we find that it does not hurt edge
prediction performance, and in several cases, improves performance.

##  5 Related work

Many real-world systems can be represented as temporal graphs, such as
transportation systems (Jiang and Luo,, 2022), human mobility patterns (Chang
et al.,, 2023), and biological networks (Prill et al.,, 2005). While most GNNs
are designed for static graphs, in recent years, there has been growing
interest in developing GNNs for temporal graphs, as noted by the Temporal
Graph Benchmark (Huang et al.,, 2023) and surveys on temporal GNNs (Skarding
et al.,, 2021; Longa et al.,, 2023). However, to the best of our knowledge,
GNNs have not yet been designed to handle temporal production graphs (TPGs),
as described in this work. Furthermore, only a handful of works have explored
GNNs for supply chains: Aziz et al., (2021) and Kosasih and Brintrup, (2022)
use GNNs to predict hidden links between firms in the static supply chain
network, and Wasi et al., (2024) provide a benchmark dataset for GNNs,
consisting of a graph where nodes are products and edges represent
relationships between products (e.g., same product group, same storage
location). In contrast, we construct a dynamic supply chain network consisting
of firm and product nodes, and our goals are to learn production functions and
to forecast future transactions on this network.

Our work also has connections to temporal causal discovery (Nauta et al.,,
2019; Löwe et al.,, 2022; Assaad et al.,, 2022), where the goal is to learn a
causal graph between variables based on their timeseries. Here, the variables
are the products, but our setting is different since the timeseries appear as
in-edges and out-edges embedded within the supply chain network, with each
product having multiple timeseries (for different firms that buy and supply
it). Our work is also related to inferring networks from node marginals (Kumar
et al.,, 2015; Maystre and Grossglauser,, 2017; Chang et al.,, 2024) and
combining deep learning with domain-specific principles, a theme which appears
in physics simulation (Wang et al.,, 2021) and epidemiological forecasting
(Liu et al.,, 2024).

Within the supply chains domain, our work builds on prior literature that
represents supply chains as networks (Fujiwara and Aoyama,, 2010) and models
the propagation of shocks over these networks (Acemoglu et al.,, 2012; Zhao et
al.,, 2019; Li et al.,, 2020; Carvalho et al.,, 2021). Our work is unique in
two key ways: first, much of the prior work relies on theoretical models and
synthetic networks, and, even among empirical studies, the data used is
typically industry-level or, at best, firm-level with static relations between
firms (Carvalho and Tahbaz-Salehi,, 2019). In contrast, we have access to
transaction-level data, which reveals essential time-varying information.
Second, integration of ML into supply chains has been limited, and only
recently have researchers begun exploring the use of ML in supply chains, with
calls for further exploration (Baryannis et al.,, 2019; Brintrup et al.,,
2020; Younis et al.,, 2022). Most prior models for modeling supply chain
propagation are mechanistic (Hallegatte,, 2008; Guan et al.,, 2020; Inoue and
Todo,, 2020; Li et al.,, 2021) and limited in their ability to even fit
aggregate counts, such as country-level GDP (Inoue and Todo,, 2019). In
contrast, by combining mechanistic models with GNNs, we can forecast
individual transactions with remarkable accuracy, while maintaining the
interpretability of the inventory module. Our model can also be used for
critical supply chain applications, such as demand forecasting (Seyedan and
Mafakheri,, 2020), early detection of risks (Sheffi,, 2015), and inventory
optimization (Vandeput,, 2020).

##  6 Conclusion

In this work, we have formalized TPGs, with supply chain networks as a
canonical example, and developed a new class of GNNs to handle them. Our
models successfully achieve two essential objectives in this setting—inferring
production functions and predicting future edges—while preexisting GNNs
focused on the latter. Furthermore, we build a new simulator, SupplySim, to
evaluate models under varied and controlled TPG settings, and release the
simulator and synthetic datasets to encourage further research development. In
future work, we hope to apply our model to TPGs in other domains, such as
biological, and to develop theoretical results for our inventory module, such
as establishing identifiability conditions and connecting it to causal
inference.

## Acknowledgements

S.C. was supported in part by an NSF Graduate Research Fellowship, the Meta
PhD Fellowship, and NSF award CCF-1918940. J.L. was supported in part by NSF
awards OAC-1835598, CCF-1918940, DMS-2327709, and Stanford Data Applications
Initiative. From Hitachi, we are thankful to Arnab Chakrabarti for his
continuous support.

## References

  * Acemoglu et al., (2012) Acemoglu, D., Carvalho, V. M., Ozdaglar, A., and Tahbaz-Salehi, A. (2012).  The network origins of aggregate fluctuations.  Econometrica, 80(5):1977–2016.
  * Aigner and Chu, (1968) Aigner, D. J. and Chu, S. F. (1968).  On estimating the industry production function.  The American Economic Review, 58(4):826–839.
  * Assaad et al., (2022) Assaad, C. K., Devijver, E., and Gaussier, E. (2022).  Survey and evaluation of causal discovery methods for time series.  Journal of Artificial Intelligence Research, 73.
  * Aziz et al., (2021) Aziz, A., Kosasih, E. E., Griffiths, R.-R., and Brintrup, A. (2021).  Data considerations in graph representation learning for supply chain networks.  In ICML 2021 Workshop on Machine Learning for Data.
  * Baryannis et al., (2019) Baryannis, G., Dani, S., and Antoniou, G. (2019).  Predicting supply chain risks using machine learning: The trade-off between performance and interpretability.  Future Generation Computer Systems, 101:993–1004.
  * Baumgartner et al., (2020) Baumgartner, T., Malik, Y., and Padhi, A. (2020).  Reimagining industrial supply chains.  McKinsey & Company.
  * Brintrup et al., (2020) Brintrup, A., Pak, J., Ratiney, D., Pearce, T., Wichmann, P., Woodall, P., and McFarlane, D. (2020).  Supply chain data analytics for predicting supplier disruptions: a case study in complex asset manufacturing.  International Journal of Production Research, 58(11):3330–3341.
  * Bruch et al., (2019) Bruch, S., Wang, X., Bendersky, M., and Najork, M. (2019).  An analysis of the softmax cross entropy loss for learning-to-rank with binary relevance.  In Proceedings of the 2019 ACM SIGIR International Conference on the Theory of Information Retrieval (ICTIR’19).
  * Carvalho et al., (2021) Carvalho, V. M., Nirei, M., Saito, Y. U., and Tahbaz-Salehi, A. (2021).  Supply chain disruptions: Evidence from the great east japan earthquake.  The Quarterly Journal of Economics, 136(2):1255–1321.
  * Carvalho and Tahbaz-Salehi, (2019) Carvalho, V. M. and Tahbaz-Salehi, A. (2019).  Production networks: A primer.  Annual Review of Economics, 11:635–663.
  * Chang et al., (2024) Chang, S., Koehler, F., Qu, Z., Leskovec, J., and Ugander, J. (2024).  Inferring dynamic networks from marginals with iterative proportional fitting.  In Proceedings of the 41st International Conference on Machine Learning (ICML’24).
  * Chang et al., (2023) Chang, S., Vrabac, D., Leskovec, J., and Ugander, J. (2023).  Estimating geographic spillover effects of covid-19 policies from large-scale mobility networks.  In Proceedings of the 37th AAAI Conference on Artificial Intelligence (AAAI’23).
  * Coelli et al., (2005) Coelli, T. J., Rao, D. S. P., O’Donnell, C. J., and Battese, G. E. (2005).  An Introduction to Efficiency and Productivity Analysis.  Springer Science & Business Media.
  * Cong et al., (2023) Cong, W., Zhang, S., Kang, J., Yuan, B., Wu, H., Zhou, X., Tong, H., and Mahdavi, M. (2023).  Do we really need complicated model architectures for temporal networks?  In Proceedings of the 11th International Conference on Learning Representations (ICLR’23).
  * Fujiwara and Aoyama, (2010) Fujiwara, Y. and Aoyama, H. (2010).  Large-scale structure of a nation-wide production network.  The European Physical Journal B, 77:565––580.
  * Grover and Leskovec, (2016) Grover, A. and Leskovec, J. (2016).  node2vec: Scalable feature learning for networks.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD’16).
  * Guan et al., (2020) Guan, D., Wang, D., Hallegatte, S., Davis, S. J., Huo, J., Li, S., Bai, Y., Lei, T., Xue, Q., Coffman, D., Cheng, D., Chen, P., Liang, X., Xu, B., Lu, X., Wang, S., Hubacek, K., and Gong, P. (2020).  Global supply-chain effects of covid-19 control measures.  Nature Human Behaviour, 4:577–587.
  * Hallegatte, (2008) Hallegatte, S. (2008).  An adaptive regional input-output model and its application to the assessment of the economic cost of katrina.  Risk Analysis, 28(3):779–799.
  * Huang et al., (2023) Huang, S., Poursafaei, F., Danovitch, J., Fey, M., Hu, W., Rossi, E., Leskovec, J., Bronstein, M., Rabusseau, G., and Rabbany, R. (2023).  Temporal graph benchmark for machine learning on temporal graphs.  In Proceedings of the 37th Annual Conference on Neural Information Processing Systems (NeurIPS’23).
  * Inoue and Todo, (2019) Inoue, H. and Todo, Y. (2019).  Firm-level propagation of shocks through supply-chain networks.  Nature Sustainability, 2:841–847.
  * Inoue and Todo, (2020) Inoue, H. and Todo, Y. (2020).  The propagation of economic impacts through supply chains: The case of a mega-city lockdown to prevent the spread of covid-19.  PLoS ONE, 15(9).
  * Jiang and Luo, (2022) Jiang, W. and Luo, J. (2022).  Graph neural network for traffic forecasting: A survey.  Expert Systems with Applications: An International Journal, 207(C).
  * Kazemi et al., (2020) Kazemi, S. M., Goel, R., Jain, K., Kobyzev, I., Sethi, A., Forsyth, P., and Poupart, P. (2020).  Representation learning for dynamic graphs: A survey.  Journal of Machine Learning Research, 21(70):1–73.
  * Kosasih and Brintrup, (2022) Kosasih, E. E. and Brintrup, A. (2022).  A machine learning approach for predicting hidden links in supply chain with graph neural networks.  International Journal of Production Research, 60(17):5380–5393.
  * Kumar et al., (2015) Kumar, R., Tomkins, A., Vassilvitskii, S., and Vee, E. (2015).  Inverting a steady-state.  In Proceedings of the 8th ACM International Conference on Web Search and Data Mining (WSDM’15).
  * Li et al., (2021) Li, Y., Chen, K., Collignon, S., and Ivanov, D. (2021).  Ripple effect in the supply chain network: Forward and backward disruption propagation, network health and firm vulnerability.  European Journal of Operational Research, 291:1117–1131.
  * Li et al., (2020) Li, Y., Zobel, C. W., Seref, O., and Chatfield, D. (2020).  Network characteristics and supply chain resilience under conditions of risk propagation.  International Journal of Production Economics, 223.
  * Liu et al., (2024) Liu, Z., Wan, G., Prakash, B. A., Lau, M. S. Y., and Jin, W. (2024).  A review of graph neural networks in epidemic modeling.  In arXiv.
  * Longa et al., (2023) Longa, A., Lachi, V., Santin, G., Bianchini, M., Lepri, B., Liò, P., Scarselli, F., and Passerin, A. (2023).  Graph neural networks for temporal graphs: State of the art, open challenges, and opportunities.  In arXiv.
  * Löwe et al., (2022) Löwe, S., Madras, D., Zemel, R., and Welling, M. (2022).  Amortized causal discovery: Learning to infer causal graphs from time-series data.  In Proceedings of the 1st Conference on Causal Learning and Reasoning (CLeaR’22).
  * Maystre and Grossglauser, (2017) Maystre, L. and Grossglauser, M. (2017).  Choicerank: Identifying preferences from node traffic in networks.  In Proceedings of the 34th International Conference on Machine Learning (ICML’17).
  * Nauta et al., (2019) Nauta, M., Bucur, D., and Seifert, C. (2019).  Causal discovery with attention-based convolutional neural networks.  Machine Learning and Knowledge Extraction, 1(1):312–340.
  * Newman, (2001) Newman, M. E. J. (2001).  Clustering and preferential attachment in growing networks.  Physical Review E, 64.
  * Prill et al., (2005) Prill, R. J., Iglesias, P. A., and Levchenko, A. (2005).  Dynamic properties of network motifs contribute to biological network organization.  PLoS Biology, 3(11).
  * Rossi et al., (2020) Rossi, E., Chamberlain, B., Frasca, F., Eynard, D., Monti, F., and Bronstein, M. (2020).  Temporal graph networks for deep learning on dynamic graphs.  In ICML 2020 Workshop on Graph Representation Learning.
  * Seyedan and Mafakheri, (2020) Seyedan, M. and Mafakheri, F. (2020).  Predictive big data analytics for supply chain demand forecasting: methods, applications, and research opportunities.  Journal of Big Data, 7(53).
  * Sheffi, (2015) Sheffi, Y. (2015).  Preparing for disruptions through early detection.  MIT Sloan Management Review, 57(1).
  * Shi et al., (2021) Shi, Y., Huang, Z., Feng, S., Zhong, H., Wang, W., and Sun, Y. (2021).  Masked label prediction: Unified message passing model for semi-supervised classification.  In Proceedings of the 13th International Joint Conference on Artificial Intelligence (IJCAI’21).
  * Skarding et al., (2021) Skarding, J., Gabrys, B., and Musial, K. (2021).  Foundations and modeling of dynamic networks using dynamic graph neural networks: A survey.  IEEE Access, 9:79143–79168.
  * The White House, (2021) The White House (2021).  Executive order on america’s supply chains.  Available at https://www.whitehouse.gov/briefing-room/presidential-actions/2021/02/24/executive-order-on-americas-supply-chains/.
  * Tolstikhin et al., (2021) Tolstikhin, I., Houlsby, N., Kolesnikov, A., Beyer, L., Zhai, X., Unterthiner, T., Yung, J., Steiner, A., Keysers, D., Uszkoreit, J., Lucic, M., and Dosovitskiy, A. (2021).  Mlp-mixer: An all-mlp architecture for vision.  In Proceedings of the 35th Conference on Neural Information Processing Systems (NeurIPS’21).
  * Vandeput, (2020) Vandeput, N. (2020).  Inventory Optimization: Models and Simulations.  Walter de Gruyter GmbH, Berlin, Germany.
  * Wang et al., (2021) Wang, R., Walters, R., and Yu, R. (2021).  Incorporating symmetry into deep dynamics models for improved generalization.  In Proceedings of the 9th International Conference on Learning Representations (ICLR’21).
  * Wasi et al., (2024) Wasi, A. T., Islam, M. S., and Akib, A. R. (2024).  Supplygraph: A benchmark dataset for supply chain planning using graph neural networks.  arXiv.
  * Younis et al., (2022) Younis, H., Sundarakani, B., and Alsharairi, M. (2022).  Applications of artificial intelligence and machine learning within supply chains:systematic review and future research directions.  Journal of Modelling in Management.
  * Zhao et al., (2019) Zhao, K., Zuo, Z., and Blackhurst, J. V. (2019).  Modelling supply chain adaptation for disruptions: An empirically grounded complex adaptive systems approach.  Journal of Operations Management, 65(2):190–212.

## Appendix

##  Appendix A Details about model

We provide additional details about our model architecture and evaluation
framework here.

###  A.1 SC-TGN

Here we describe SC-TGN, our extended version of TGN (Rossi et al.,, 2020), in
greater detail. Like TGN, our model consists of the following modules:

Memory. In SC-TGN, each node has a time-varying memory $\mathbf{m}_{i}^{(t)}$.
The memory captures the node’s state, influenced by all of the transactions we
have seen for this node up to this point. We initialize each node’s memory to
a learnable vector $\mathbf{v}_{i}^{(0)}$. The node’s memory is updated each
time the node participates in a transaction.

Message function. Each transaction, $e(s,b,p,t)$, sends a message to the three
participating nodes: the supplier firm $s$, buyer firm $b$, and product $p$.
Following the original work, we represent the raw message $\textrm{msg}_{e}$
as the concatenation of the nodes’ memories and a time encoding:

$\displaystyle\textrm{msg}_{e}$
$\displaystyle=[\mathbf{m}_{s}^{(t)}|\mathbf{m}_{b}^{(t)}|\mathbf{m}_{p}^{(t)}|\textrm{enc}(t)],$
(10)

where $\textrm{enc}(t)$ is a temporal encoding of time $t$, implemented as a
simple linear layer. Then, the message to the supplier is
$M_{s}(\textrm{msg}_{e})$, where $M_{s}$ is a supplier-specific message
encoder, implemented as a linear layer followed by a ReLU activation. $M_{b}$
and $M_{p}$ are buyer-specific and product-specific message encoders,
respectively, with the same architecture as $M_{s}$.

Memory updater. First, for each firm $i$, we aggregate over the messages that
it received at time $t$ to produce $\bar{msg}_{i}^{(t)}$. Following the
original work, we simply take the mean over messages. Then, we update node
$i$’s memory with its aggregated message, using a recurrent neural network:

$\displaystyle\mathbf{m}_{i}^{(t+1)}=RNN(\mathbf{m}_{i}^{(t)},\bar{msg}_{i}^{(t)}).$
(11)

Embedding. Finally, the embedding module produces an embedding
$\mathbf{z}_{i}^{(t)}$ of node $i$ at time $t$. Following the original work,
we explore both an “ID” embedding, where the embedding is simply the node’s
memory, and embeddings that are constructed by applying a GNN to the memories.
For the GNN, we use the UniMP model (Shi et al.,, 2021), which is the model
used in the TGN implementation by Temporal Graph Benchmark (Huang et al.,,
2023). UniMP is a graph transformer architecture with multi-head attention; we
refer the reader to the original text for details. The advantage of the GNN-
based embedding over the ID embedding is that the GNN allows the node’s
embedding to adapt over time, even if the node itself has not been involved in
a transaction in a while, thus avoiding the memory staleness problem (Kazemi
et al.,, 2020). In Table 5, we show that using a GNN to construct node
embeddings from the node memories, instead of directly using the node
memories, greatly improves performance of SC-TGN.

##### New elements in SC-TGN.

SC-TGN extends TGN in the following ways:

  * •

Hyperedges. To adapt TGN to handle hyperedges, we modify the message function
so that it sends a message to three nodes, instead of two nodes. To capture
the nodes’ varying roles within each hyperedge (supplier, buyer, product), we
also use separate message encoders $M_{s}$, $M_{b}$, and $M_{p}$,
respectively. We also modified the decoder to predict relationships between
three nodes instead of two nodes.

  * •

Predicting edge weight. The original TGN only predicts edge existence, with a
decoder that produces the probability of an edge existing between two nodes.
Our model has two decoders: one to predict the probability of a transaction
existing and, conditioned on a transaction existing, one to predict the amount
of the transaction (i.e., edge weight).

  * •

Update penalty. We found it useful to regularize the memories over time, so
that the memory update (11) was encouraged to be small.

$\displaystyle\ell_{\textrm{update}}=\frac{\lambda_{\textrm{update}}}{m+n}\sum_{i\in[m+n]}||\mathbf{m}_{i}^{(t+1)}-\mathbf{m}_{i}^{(t)}||_{2}.$
(12)

  * •

Learnable initial memory. We initialize the node’s memory to a learnable
vector $\mathbf{v}_{i}^{(0)}$, instead of initializing all new nodes to the
zero vector, which is what was done in the original TGN. Making the initial
memory learnable allows the model to learn representations for nodes even
before they have participated in any transactions.

  * •

Training following the negative sampling scheme. In the original paper, the
authors train TGN with only perturbation negatives and not historical
negatives, but they test TGN on both kinds of negatives. We found that test
performance increased substantially (+10 or more points in MRR) if the model
could be also be trained on historical negatives.

###  A.2 SC-GraphMixer

Here we describe SC-GraphMixer, our extended version of GraphMixer (Cong et
al.,, 2023), in greater detail. Like GraphMixer, our model consists of a link
encoder and a node encoder:

##### Link encoder.

The purpose of the link encoder is to capture information from recent links
associated with the node.

  * •

Time encoding. Given a timestep $t$, GraphMixer encodes the timestep by
projecting it to a $d$-dimensional vector. The time-encoding function utilizes
features
$\mathbf{\omega}=\\{\alpha^{-(d^{\prime}-1)/\beta}\\}_{d^{\prime}=1}^{d}$,
where $\alpha$ and $\beta$ are predefined hyperparameters, and projects $t$ to
$\cos(t\times\mathbf{\omega})\in[-1,+1]$.

  * •

To represent a given node $i$ at time $t$, GraphMixer constructs the matrix
$\mathbf{T}_{i}(t)$. Each row of this matrix corresponds to a recent link of
node $i$’s, where the link is represented as the concatenation of
$\cos((t-t_{e})\times\mathbf{\omega})$, where $t_{e}$ is the timestep of the
link, and the link’s features (in our case, the transaction amount). Encoding
$t-t_{e}$ instead of $t_{e}$ itself captures how recent the link was, and thus
how much influence it should have over the current timestep.

  * •

Then, GraphMixer applies a 1-layer MLP-mixer (Tolstikhin et al.,, 2021) to
$\mathbf{T}_{i}(t)$:

$\displaystyle\mathbf{H}_{\textrm{token}}$
$\displaystyle=\mathbf{T}_{i}(t)+\mathbf{W}^{(2)}_{\textrm{token}}GeLU(\mathbf{W}^{(2)}_{\textrm{token}}LayerNorm(\mathbf{T}_{i}(t)))$
(13) $\displaystyle\mathbf{H}_{\textrm{channel}}$
$\displaystyle=\mathbf{H}_{\textrm{token}}+GeLU(LayerNorm(\mathbf{H}_{\textrm{token}})W^{(1)}_{\textrm{channel}})W^{(2)}_{\textrm{channel}}.$
(14)

Finally, node $i$’s link encoding at $t$, $\mathbf{t}_{i}(t)$, is the mean of
$\mathbf{H}_{\textrm{channel}}$.

##### Node encoder.

The purpose of the node encoder is to capture the node’s identity and its
features, along with its neighbors’ features. To represent a given node $i$ at
time $t$, its node encoding, $\mathbf{s}_{i}(t)$, is

$\displaystyle\mathbf{s}_{i}(t)$
$\displaystyle=\mathbf{x}_{i}^{\textrm{node}}+\frac{1}{|\mathcal{N}(i;t-T,t)|}\sum_{j\in\mathcal{N}(i;t-T,t)}\mathbf{x}_{j}^{\textrm{node}},$
(15)

i.e., the sum of the node’s own features and the mean of the node’s neighbors’
features, where neighbors from time $t-T$ to $t$ are included and $T$ is a
predefined hyperparameter.

Finally, GraphMixer’s node embedding of node $i$ at time $t$ is the
concatenation of its node encoding and link encoding:
$\mathbf{z}_{i}^{(t)}=[\mathbf{s}_{i}(t)|\mathbf{t}_{i}(t)]$.

##### New elements in SC-GraphMixer.

SC-GraphMixer extends GraphMixer in the following ways:

  * •

Hyperedges. To adapt GraphMixer to handle hyperedges, we compute both link and
node embeddings for all three nodes instead of two nodes, and the input to the
decoder is the concatenation of their three embeddings (7). For a product
node, its recent links involve all recent transactions that it was involved
in. For a firm node, its recent links involve all recent transactions where it
was the supplier or the buyer. As in GraphMixer, the time encoding encodes the
difference in time between when the transaction occurred and the current
timestep.

  * •

Learnable node features. The original GraphMixer treats node features as per-
node fixed input to models. Inspired by the memory module of TGN, we
introduced learnable node feature $\mathbf{x}_{i}^{\textrm{node}}$ for each
node, when the dataset does not provide raw node features. The node features
are randomly initialized and optimized as model parameters. Introducing
learnable node features also means that the static node embedding is a
strictly nested version of SC-GraphMixer, making it a proper ablation of both
of our models.

  * •

Predicting edge weight in addition to edge existence, as described above for
SC-TGN.

  * •

Training following the negative sampling scheme, as described above for SC-
TGN.

###  A.3 Penalties and caps from inventory module

As described in the main text, when our inventory module is combined with a
GNN, such as SC-TGN or SC-GraphMixer, we can optionally use the inventory
module to help with predicting future edges, in addition to its primary
purpose, which is to learn production functions. To help with predicting edge
existence, the inventory module penalizes impossible transactions, as
described in (8), and adds its penalty to the original prediction $\hat{y}$.
To help with predicting edge weight, the inventory module computes a cap,
i.e., the maximum amount possible in the transaction, as described in (9), and
takes the minimum of its cap and $\hat{y}$. Note that, if the inventory
module’s attention weights match the ground-truth production functions and the
inventories reflect what the firm truly has at the time of the transaction,
these penalties and caps can only help, not hurt, the prediction. This is
because the penalty (8) will always be 0 for a true transaction while it can
be negative for a false transaction, so adding the penalties can only help.
Similarly, the amount cap (9) is guaranteed to be greater than or equal to the
true transacted amount, so taking the minimum of the cap and the original
predicted amount can also only help.

However, even if the learned attention weights match the ground-truth
production functions closely, it is not guaranteed that the inventory reflects
what each firm truly has at the time of the transaction. One minor case that
we see of this is that the inventory is only updated up to time $t$, but the
current batch could include transactions at time $t+1$ or later; thus, we only
apply penalties or caps to transactions in time $t$, since at time $t+1$ and
later, the inventory will look different. A more serious case is when we are
systematically missing transactions in the dataset, so as a result, we are
underestimating the number of products that each firm has in its inventory, so
penalties are too large and caps are too low. We investigate this issue in
Appendix C.2 and show that, when some firms are missing from the synthetic
data, the penalties and caps indeed hurt edge prediction performance, even if
the attention weights are perfectly learned (we set them to the ground-truth
production functions). So, it is preferable not to use the penalties and caps
in those settings. In future work, it would be useful to explore how to make
the inventory module more robust to missing data, so that we do not
underestimate the number of products in the firms’ inventories even when
transactions are missing. For example, some smoothing could be applied such
that when we observe a firm buying product $p$, we also add a bit of
$p^{\prime}$ to its inventory, where $p^{\prime}$ is a “similar” product,
where similarity could be defined based on co-buying patterns.

###  A.4 Model training and evaluation

In summary, our model loss consists of four components,

$\displaystyle\mathcal{L}$
$\displaystyle=\ell_{\textrm{exist}}+\ell_{\textrm{weight}}+\ell_{\textrm{inv}}+\ell_{\textrm{update}},$
(16)

which correspond to the softmax cross-entropy loss from predicting edge
existence; RMSE loss from predicting edge weight (20); the inventory loss (6);
and an update loss that regularizes how much the node memories change between
consecutive timesteps (12). Aside from $\ell_{\textrm{exist}}$, each of the
other losses has hyperparameters that control their relative weight. In Table
4, we document the hyperparameters we used in our experiments.

##### Learning production functions.

We use mean average precision (MAP) to evaluate how well the inventory module
learned the ground-truth production functions. Let $\mathcal{P}_{p_{o}}$
represent the ground-truth set of parts for product $p_{o}$, let
$\alpha_{p_{o}}=\\{\alpha_{p_{o}p_{1}},\alpha_{p_{o}p_{2}},\cdots,\alpha_{p_{o}p_{m}}\\}$
represent the inventory module’s learned attention weights for $p_{o}$, and
let $\textrm{pos}(\alpha_{p_{o}},k)$ represent the product at rank $k$ in
$\alpha_{p_{o}}$. We define average precision as

$\displaystyle\textrm{AvePrec}(p_{o})$
$\displaystyle=\frac{1}{|\mathcal{P}_{p_{o}}|}\sum_{k=1}^{m}\textrm{Prec@K}(\alpha_{p_{o}},\mathcal{P}_{p_{o}},k)\cdot\mathbbm{1}[\textrm{pos}(\alpha_{p_{o}},k)\in\mathcal{P}_{p_{o}}]$
(17) $\displaystyle\textrm{Prec@K}(\alpha,\mathcal{P},k)$
$\displaystyle=\frac{1}{k}\sum_{k^{\prime}=1}^{k}\mathbbm{1}[\textrm{pos}(\alpha,k^{\prime})\in\mathcal{P}].$
(18)

In other words, average precision takes the average of Precision@K over all
ranks $k$ where there is a relevant item. Average precision and Precision@K
are both common metrics in information retrieval, but we prefer average
precision here since it can account for varying numbers of relevant items over
queries, which applies to our setting since the size of $\mathcal{P}_{p_{o}}$
varies across products.

##### Predicting edge existence.

As described in the main text, we sample two types of negatives, perturbation
negatives and historical negatives, to accompany each positive transaction.
Temporal Graph Benchmark (Huang et al.,, 2023) also samples both types of
negatives, although we modify the procedure slightly for hyperedges. For a
positive transaction $e(s,b,p,t)$, we have three types of perturbation
negatives: $e(s^{\prime},b,p,t)$, where a random supplier is sampled;
$e(s,b^{\prime},p,t)$, where a random buyer is sampled; and
$e(s,b,p^{\prime},t)$, where a random product is sampled. A historical
negative is a transaction $e(s^{\prime},b^{\prime},p^{\prime},t)$ such that
$s^{\prime},b^{\prime},p^{\prime}$ appeared in the train set, but not at time
$t$. To make historical negatives even harder, we sample historical negatives
proportionally to their count in the train set. For both types of negatives,
we ensure that the sampled negative did not actually occur at time $t$. We
sample 9 negatives following the perturbation procedure (with 3 of each type)
and 9 negatives following the historical procedure. It is possible, though,
that more than 9 negatives are historical (if the perturbation procedure
results in a historical negative) or more than 9 negatives are perturbations
(if a historical negative happens to be a one-node perturbation).

We use mean reciprocal rank (MRR) to evaluate performance for predicting edge
existence. For each positive transaction $e$ and its corresponding set of
negatives $\mathcal{N}_{e}$ in batch $B$, MRR quantifies the position of the
positive transaction among the negative samples:

$\displaystyle MRR=\frac{1}{|B|}\sum_{e\in
B}\left(\frac{\sum_{n\in\mathcal{N}_{e}}\mathbbm{1}[\hat{y}_{n}<\hat{y}_{e}]+\sum_{n\in\mathcal{N}_{e}}\mathbbm{1}[\hat{y}_{n}\leq\hat{y}_{e}]}{2}+1\right)^{-1}.$
(19)

Following the Temporal Graph Benchmark (Huang et al.,, 2023), we compute an
“optimistic” MRR where we tie-break in favor of the positive transaction and a
“pessimistic” MRR where we tie-break in favor of the negative transactions,
and take the average of the two. We refer to 0.1 as the random baseline, since
we sample 18 negatives in total and, if all transactions received the same
score, then the MRR would be $((0+18)/2+1)^{-1}=0.1$. Similarly, if we
assigned uniform random scores to all 19 alternatives, the expected rank of
the positive is 10 (9 negatives before and 9 negatives after), so the
reciprocal is 0.1.

##### Predicting edge weight.

We preprocess transaction amounts by applying log-scaling then standard
scaling, where we subtract the mean and divide by the standard deviation, so
that amounts now have zero mean and variance of 1. We use the scaled amount as
an edge feature, to stabilize training, and as a prediction target, to make
RMSE comparable across datasets. However, note that we use the raw amount to
update the inventory module, as required in (1)-(2). Then, we compare the
predicted amount $\hat{y}$ to the true amount, $y$, using root mean squared
error:

$\displaystyle RMSE=\sqrt{\frac{1}{|B|}\sum_{e\in
B}(amt(e)-\hat{y}_{e})^{2}}.$ (20)

Since RMSE is differentiable, we use it both as the training loss and
evaluation metric.

##  Appendix B Supply chains data

Attribute Counts |  SS-std |  Tesla |  IED
---|---|---|---
# Product Nodes |  50 |  2,690 |  3,029
# Firms Nodes |  119 |  11,628 |  2,583
# Transactions |  71646 |  581,002 |  279,712
Timespan (Days) |  198 |  1683 |  359
Table 3: Dataset statistics.

###  B.1 Real-world supply chains data

For our experiments, we constructed two large-scale, real-world supply chain
datasets by compiling data from TradeSparq, a third-party data provider. We
rely on a wide variety of data sources from authorized government sources
across over 60 countries, including authorities like the Automated Manifest
System, operated by US Customs and Border Protection, and other government
customs departments. The data appear in the form of bills of lading (BoL),
receipts of reported transactions, and custom declarations. Each transaction
includes information such as data source, buyer title, supplier title, cost of
reported transaction, product (in terms of HS code and description), import
and export country, amount (weight and quantity unit), among other columns.

We have access to the data through an API, which allows us to make queries
with associated costs. Due to these costs, we cannot reconstruct the entire
supply chain network, but we focus on certain products and grow the network
outwards via a breadth-first-search-like process. The construction begins by
selecting product-specific HS codes and identifying prominent Tier-0 supplier
firms involved in high-volume trading. The process captures all transactions
from Tier-0 suppliers, then aggregates their group companies to identify
Tier-1 suppliers, recording all corresponding transactions. This method is
repeated to determine Tier-2 suppliers by analyzing the group companies of
Tier-1 firms. All transactions between Tier-1 and Tier-2 firms are also
collected. This systematic approach ensures a thorough mapping of supplier
interactions up to Tier-2, providing detailed insights into the supply chain
network. We illustrate our process as a flowchart in Figure 3.

Start Select product HS codes Select high amount trading prominent Tier-0
supplier firms of the product Get all transactions between Tier-0 firms
setting as suppliers and their buyers Aggregate group companies of the Tier-0
firms & setting them as buyers get all Tier-1 supplier firms Get all
transactions between Tier-0 & Tier-1 firms Aggregate group companies of the
Tier-1 firms & setting them as buyers get all Tier-2 supplier firms Get all
transactions between Tier-1 & Tier-2 firms End

Figure 3: Overview of supply chain network data construction process.

##### Automotive dataset.

Following this process, we constructed real-world automotive dataset within
the Tesla, Inc. group companies involved in the electric vehicle (EV) and EV
parts supply chain. The analysis covered transactions from January 1, 2019, to
December 31, 2022. This process started by identifying Tesla EV maker firms
and capturing all their transactions with Tier-1 suppliers, i.e. Tesla’s
direct suppliers. Subsequently, Tier-1 suppliers were identified by
aggregating Tier-1 firm groups and analyzing their transactions with Tier-2
supplier firms. This structured approach provided a comprehensive mapping of
Tesla’s supply chain interactions, offering detailed insights into the
dynamics and dependencies within their EV supply chain network.

##### Industrial equipment dataset.

This dataset focuses on microscopes along with other specialized analytical
and inspection equipment and their manufacturers. Spanning the years 2022 to
2023, the dataset covers approximately 630 makers. The process began with
identifying Tier-0 suppliers associated with these specialized products and
capturing all relevant transactions with their buyers. Further deep-tiers were
constructed by aggregating transactions at Tier-1 and Tier-2 supplier levels,
providing a detailed insight into the supply chain dynamics and the extensive
network of suppliers involved in the production of these critical industrial
tools. To construct the bill of material (BoM) for microscopes,we make use of
an extensive survey of direct suppliers, inquiring about their products and
their own suppliers. This survey, conducted in Japanese, required translating
each supplied part name from Japanese to English. We then mapped these part
names to HS codes using the HTS standard and manually verified the mappings.
We successfully mapped and verified 87% of parts to HS codes, due to
challenges posed by difficult and miscellaneous part names. Of the mapped
parts, we could find 75% parts within our constructed supply chain.
Consequently, we were able to construct 65% of the BoM for microscopes, which
we then used to evaluate our inventory module.

##### Limitations of real-world supply chains data.

TradeSparq data provides a rare opportunity to access transaction-level
information, while most supply chain studies are done on much more aggregated
data (e.g., imports and exports aggregated over country) or, at best, static
firm-level networks, without time-varying transactions (Inoue and Todo,, 2019;
Carvalho and Tahbaz-Salehi,, 2019; Guan et al.,, 2020). However, there are
still limitations of the real-world data, which motivate the need for
synthetic data. We list some of these limitations below:

  * •

Missing transactions. While TradeSparq covers 60+ countries, many transactions
are also missing, especially domestic transactions that would not be reported
in customs. Some companies may also deliberately obscure or omit details of
their sensitive transactions in publicly available data to protect their
business interests.

  * •

Missing production functions. For most firms and products, we are missing
their true production functions, which limits our ability to evaluate our
inventory module.

  * •

Proprietary data. Since the TradeSparq data is proprietary, we are not allowed
to release it.

  * •

Incomplete records. Essential fields such as transaction suppliers,
quantities, and amounts may be missing, making it difficult to construct a
complete picture of the transaction depending on the reporting country
authority.

  * •

Duplicate entries can arise from multiple sources reporting the same
transaction.

  * •

Transactions are reported at different times than actual production dates.

  * •

Due to logistic constraints, multiple parts or products may be reported
together.

###  B.2 Details about SupplySim

To address the limitations of real-world supply chains data, we build a
simulator, SupplySim, which offers the following benefits:

  * •

We can test the model under complete data settings, where all transactions are
observed and all necessary data per transaction is provided. This allows us to
disentangle model performance from data quality.

  * •

We can also test how model performance changes as data quality changes,
controlling how and to what extent it changes, such as dropping 20% of firms
at random.

  * •

We can generate any number of shocks, of varying size and source, and test the
model’s ability to predict how transactions change under the shocks. In real
data, large shocks are rare and confounded by many other contemporaneous
factors so it is difficult to isolate their effect (e.g., shocks during the
COVID-19 pandemic).

  * •

We have production functions for all products, which enables us to evaluate
the inventory module, and we can, again, vary the data settings to test how
robust the inventory module’s performance is to missing data.

  * •

Unlike the real-world data, we can release SupplySim and our generated
synthetic data, improving the reproducibility and transparency of our
research, along with encouraging future research in this domain.

Thus, SupplySim addresses all the limitations of real-world supply chains
data, at the cost of being synthetic, instead of real-world, data. To mitigate
this cost, we develop SupplySim to be as realistic as possible: for example,
we ensure that the generated networks match real-world networks on key
characteristics, such as degree distribution and clustering, and we use a
realistic agent-based model from economics to generate time-varying
transactions over the network.

As described in the main text, our simulator progresses in three steps: (1)
constructing the production graph, $G_{\textrm{prod}}$, which describes static
product-part relations, (2) assigning products to supplier firms and assigning
supplier firms to buyer firms, which are also static relations, (3) generating
time-varying transactions between firms.

####  B.2.1 Static graphs

##### Constructing the production graph, $G_{\textrm{prod}}$.

First, we define the number of exogenous products $n_{\textrm{exog}}$, number
of consumer products $n_{\textrm{consumer}}$, number of products per inner
tier $n_{\textrm{tier}}$, and number of inner tiers $T$. Then, the products
are assigned to tiers as follows:

  * •

Tier 0: product 0 to product $n_{\textrm{exog}}-1$,

  * •

Tier 1: product $n_{\textrm{exog}}$ to product
$n_{\textrm{exog}}+n_{\textrm{tier}}-1$,

  * •

Tier 2: product $n_{\textrm{exog}}+n_{\textrm{tier}}$ to product
$n_{\textrm{exog}}+2\cdot n_{\textrm{tier}}-1$,

  * •

$\cdots$

  * •

Tier $T$: product $n_{\textrm{exog}}+(T-1)\cdot n_{\textrm{tier}}$ to product
$n_{\textrm{exog}}+T\cdot n_{\textrm{tier}}-1$,

  * •

Tier $T+1$: product $n_{\textrm{exog}}+T\cdot n_{\textrm{tier}}$ to product
$n_{\textrm{exog}}+T\cdot n_{\textrm{tier}}+n_{\textrm{consumer}}-1$.

Thus, there are $n_{\textrm{exog}}+T\cdot
n_{\textrm{tier}}+n_{\textrm{consumer}}$ products and $T+2$ tiers overall. In
our experiments, we use $n_{\textrm{exog}}=5$, $n_{\textrm{consumer}}=5$,
$n_{\textrm{tier}}=10$, and $T=4$, resulting in 50 products. We also sample a
2-dimensional position from Uniform(0, 1) for each product.

Figure 4: Visualizing products in our synthetic datasets. Each point
represents the position of one of the 50 products, and points are color-coded
by the product’s tier. We also denote part-product relations between Tier 1
and Tier 2 products, where an arrow from product $p_{1}$ to product $p_{2}$
means that $p_{1}$ is required to make $p_{2}$. For each product, we sample
its number of parts from {1, 2, 3, 4} uniformly, then assign its parts to the
closest products in the previous tier, resulting in commonly co-occurring
parts.

For each product, we sample its number of parts uniformly at random from a
predefined range (we use 2 to 4, inclusive). Then, we select parts for the
product, where the viable parts come from the previous tier. We assign the
products to parts based on the distance between their positions, so that parts
closer to the products are likelier to be chosen, with the probability
proportional to $d^{\gamma}$, where $d$ is their distance. In our experiments,
we deterministically assign products to their closest parts for simplicity,
which is essentially equivalent to choosing an extreme $\gamma$, but it would
be straightforward to introduce more randomness in the selection. In Figure 4,
we visualize our products’ positions and the part-product relations between
Tier 1 and Tier 2, demonstrating how this assignment procedure results in
commonly co-occurring parts, while every product still has a unique set of
parts. Commonly co-occurring parts are realistic, such as nails and bolts, and
they create possible “similar” parts that can be learned by GNNs. For each
pair of input part $p_{i}$ and output product $p_{o}$, we also sample $u_{io}$
uniformly at random (we use 1 to 4, inclusive), where $u_{io}$ is the number
of units of the part $p_{i}$ needed to make one unit of the product $p_{o}$.

##### Constructing supplier-buyer relations.

Now, we define the number of firms per “group”, $n_{\textrm{group}}$, where
each group is assigned a set of $n_{\textrm{consec}}$ consecutive tiers, such
that the firms in the group are only allowed to supply products in this set of
tiers. We define groups and tier sets as follows:

  * •

Group 0: Tier 0 to Tier $n_{\textrm{consec}}-1$,

  * •

Group 1: Tier 1 to Tier $n_{\textrm{consec}}$,

  * •

Group 2: Tier 2 to Tier $n_{\textrm{consec}}+1$,

  * •

$\cdots$

  * •

Group $T-n_{\textrm{consec}}+2$: Tier $T-n_{\textrm{consec}}+2$ to Tier $T+1$.

In our experiments, we use $n_{\textrm{group}}=30$ and
$n_{\textrm{consec}}=2$, resulting in 120 firms overall. For each firm, we
also sample a 2-dimensional position from Uniform(0, 1).

Then, for each product, we sample its number of suppliers uniformly at random
from a predefined range (we use 4 to 8, inclusive). Its set of viable supplier
firms are the firms within the groups that supply products in this product’s
tier. For example, in our experiments, a product at Tier 2 can be supplied by
firms in Groups 1 and 2 (so 60 firms). Similar to how we assign parts to
products, we assign suppliers to products based on the closest distances
between the product and firm’s positions. Our procedure for assigning
suppliers to products, where we restrict firms to consecutive tiers, results
in more realistic supply chain graphs, since firms tend to specialize and we
would not expect the same firm to produce products at very different tiers.
Furthermore, by restricting firms to tiers and choosing suppliers based on
closest distances, we also result in firms that supply similar sets of
products, which is another form of similarity that can be captured in GNN
representations.

Based on the assignment of suppliers to products, and the existing production
graph $G_{\textrm{prod}}$, we can now determine which input products each firm
needs to buy in order to make the products it supplies. We iterate through
pairs $(b,p)$, where firm $b$ needs to buy product $p$, and select a supplier
firm $s$ from all suppliers of $p$, with probability proportional to $s$’s
current number of buyers. This assignment mechanism, known as preferential
attachment (Newman,, 2001), reflects a rich-get-richer dynamic that results in
power law degree distributions, which are known to appear in real-world supply
chain networks (Fujiwara and Aoyama,, 2010; Inoue and Todo,, 2019; Zhao et
al.,, 2019). Real-world firm-firm networks are also reported to have
disassortativity (negative correlation between the degrees of connected nodes)
and low clustering, which we also observe in our generated firm-firm networks:
the degree assortativity is -0.253 and average clustering coefficient is
0.246, which is similar to the average clustering coefficient of 0.179
reported in Zhao et al., (2019).

####  B.2.2 Generating time-varying transactions

Now, given the static graphs, we want to generate time-varying transactions
between firms. We base our simulation model on the ARIO model, an agent-based
model widely used in economics to simulate propagation over supply chains
(Hallegatte,, 2008; Inoue and Todo,, 2019; Guan et al.,, 2020). At any point
in the simulation, we keep track of each firm’s inventory
$\mathbf{x}_{i}^{(t)}$; all incomplete orders $\mathcal{I}^{(t)}$, where an
order is defined as $o(b,s,p,k)$, from buyer firm $b$ to supplier firm $s$ for
$k$ amount of product $p$; and all orders that were newly completed at time
$t$, $\mathcal{C}^{(t)}$. Then, at each timestep $t$, we iterate through
firms, with each firm $f$ taking the following actions.

  * •

New supply: first, the firm adds its newly received products to its inventory,
based on completed orders in $\mathcal{C}^{(t-1)}$ where it was the buyer. So,
the firm is only allowed to use products from orders completed at time $t-1$
at time $t$ or later.

  * •

Production: then, the firm goes through its incomplete orders and tries to
complete as many as it can before it runs out of inventory. Each time the firm
completes an order $o(b,f,p,k)$, it adds that order to $\mathcal{C}^{(t)}$,
removes it from $\mathcal{I}^{(t)}$, and subtracts $u_{p_{i}p}\cdot k$ of
input product $p_{i}$ from its inventory, for all inputs $p_{i}$ of product
$p$. The firm goes through its incomplete orders from oldest to latest, and
stops when $u_{p_{i}p}\cdot k$ exceeds what it has in its inventory. Thus,
orders are completed first-in-first-out. This is a reasonable principle, since
earlier orders should be completed earlier, but extensions could include
prioritization or trying smaller, later orders after larger, earlier orders
are no longer feasible.

  * •

New demand: finally, the firm places new orders to its suppliers. First, for
each input product $p_{i}$, it calculates the necessary amount
$k_{f,p_{i}}^{(t)}$ as the amount it would need to complete all of its own
incomplete orders minus the amount that is already pending in its incomplete
orders to its suppliers:

$\displaystyle k_{f,p_{i}}^{(t)}$
$\displaystyle=\sum_{(b,f,p,k)\in\mathcal{I}^{(t)}}u_{p_{i}p}\cdot
k-\sum_{(f,s,p_{i},k)\in\mathcal{I}^{(t)}}k.$ (21)

Then, the firm adds a new order, $o(f,s,p,k_{f,p_{i}}^{(t)})$ to
$\mathcal{I}^{(t)}$, where $s$ is the firm’s default supplier for product
$p_{i}$ with probability $p_{\textrm{default}}$ and, with probability
$1-p_{\textrm{default}}$, we select a supplier of $p_{i}$ uniformly at random.
We use $p_{\textrm{default}}=0.8$, reflecting a strong level of “stickiness”
in supplier-buyer relations, while allowing some randomness or one-off deals.

All completed orders at time $t$, $\mathcal{C}^{(t)}$, are reported as
transactions for that timestep. The simulation also requires special logic,
which we describe below, for products in Tier 0, which have no parts, and
products in Tier $T+1$ (the final tier), which are not used as parts for any
other products and thus only bought by consumers.

##### Tier 0 products and exogenous supply.

Figure 5: Visualizing the number of transactions per timestep, over the three
synthetic datasets: SS-std, SS-shocks, and SS-missing.

For each Tier 0 product $p_{0}$, we define time-varying exogenous supply
$supply(p_{0},t)$, such that any firm $f$ that supplies $p_{0}$ can only
complete orders $o(b,f,p_{0},k)$ where $k\leq supply(p_{0},t)$. In the
standard synthetic data (SS-std), exogenous supply is plentiful, so orders for
Tier 0 products are completed almost immediately. However, we can also
introduce supply shocks, which we model as follows: each Tier 0 product at
each timestep $t$ has an independent probability $p_{\textrm{shock}}$ of
experiencing a shock, at which point
$supply(p_{0},t)=supply_{\textrm{shock}}$. Then, we allow the supply to
recover over the subsequent timesteps, by multiplying it in each timestep by
some recovery rate $r>1$, unless it experiences another shock.

In our synthetic data with shocks, SS-shocks, we use
$p_{\textrm{shock}}=0.01$, resulting in 10 shocks over 200 timesteps: shocks
to product 0 at times 32 and 144, product 2 at times 42 and 155, product 3 at
times 51, 90, and 144, and product 4 at times 16, 58, and 131. We allow a
recovery rate of $r=1.25$, resulting in a recovery time of approximately 30
timesteps from the initial shock to full recovery (since we set the stable
supply to 1000x of $supply_{\textrm{shock}}$). In Figure 5, we visualize the
number of transactions over time in each synthetic data (SS-std, SS-shocks,
and SS-missing). While the missing firms reduce the number of transactions by
a consistent fraction, the shocks dramatically disrupt the pattern of
transactions, slowing them down whenever there is a shortage then trying to
catch up whenever the Tier 0 products have recovered.

##### Final tier products and exogenous demand.

For each product $p$ in the final tier (i.e., Tier $T+1$), we define time-
varying exogenous demand $d(p,t)$, as well as firm-specific demand $d(f,p,t)$.
To simulate demand from consumers, with temporal patterns, we assign these
products to one of three types: uniform, weekday, or weekend. We define demand
as follows:

$\displaystyle\nu$ $\displaystyle\sim\mathcal{N}(0,0.1)$ (22)
$\displaystyle\lambda$ $\displaystyle=\begin{cases}1\textrm{ if type =
uniform},\\\ 2\textrm{ if type = weekday and $t$ mod $7<5$},\\\ 0.5\textrm{ if
type = weekday and $t$ mod $7\geq 5$},\\\ 0.5\textrm{ if type = weekend and
$t$ mod $7<5$},\\\ 2\textrm{ if type = weekend and $t$ mod $7\geq 5$}.\\\
\end{cases}$ (23) $\displaystyle d(p,t)$
$\displaystyle=\lambda\cdot(d(p,t-1)+\nu)$ (24) $\displaystyle d(f,p,t)$
$\displaystyle=\textrm{Poisson}(d(p,t))$ (25)

In other words, $\nu$ captures random temporal drift and $\lambda$ captures
fixed weekly patterns, where uniform is the same over the week, weekday favors
the first five days of the week, and weekend favors the last two days. Then,
$d(f,p,t)$ is added as an order with buyer “consumer” to the incomplete orders
$\mathcal{I}^{(t)}$ at the beginning of each timestep $t$. In our experiments,
we allow demand to vary steadily like this, without any sharp changes, but it
would be possible to simulate changes in consumer demand (e.g., due to new
trends) by changing the demand schedule for final tier products.

##  Appendix C Details about experiments

Here we provide details about our experimental set-up, including our baselines
and model training / testing procedures.

###  C.1 Production learning experiments

##### Temporal correlations.

If product $p_{i}$ is a part of product $p_{o}$, then we expect that some
firms buy $p_{i}$ and supply $p_{o}$, and it would be natural for those
timeseries to be correlated, with buying slightly preceding supplying. So, for
each product pair $p_{1}$ and $p_{2}$, first we find all firms that buy
$p_{1}$ and supply $p_{2}$. If there is no such firm, then we give the pair a
score of 0. If there are such firms, then for each firm $f$, we compute the
buying timeseries of $p_{1}$, i.e., the total amount bought by $f$ of $p_{1}$
on $t=1,2,\cdots$, along with the supply timeseries of $p_{2}$, i.e., the
total amount supplied by $f$ of $p_{2}$ on $t=1,2,\cdots$. Then, we compute
the maximum correlation between the buying and supply timeseries, allowing the
supplying timeseries to lag the buying timeseries by 0 to 7 timesteps.
Finally, the score given to the pair is the average maximum correlation over
firms that buy $p_{1}$ and supply $p_{2}$.

##### Pointwise Mutual Information (PMI).

If product $p_{i}$ is a part of product $p_{o}$, then we would also expect
that firms that supply $p_{o}$ are unusually likely to buy $p_{i}$, compared
to the base rate of firms buying $p_{i}$. This idea is naturally quantified by
PMI, where we compare the joint probability to the product of the individual
probabilities. Let $buy(p_{1})$ represent whether a firm buys $p_{1}$ and
$supply(p_{2})$ represent whether a firm supplies $p_{2}$. Then, we compute
PMI as

$\displaystyle PMI(p_{1},p_{2})$
$\displaystyle=\log\left(\frac{\Pr(buy(p_{1})\wedge
supply(p_{2}))}{\Pr(buy(p_{1}))\cdot\Pr(supply(p_{2}))}\right),$ (26)

where the probabilities are computed as fractions of all firms.

##### node2vec.

If product $p_{i}$ is a part of product $p_{o}$, then we would also expect
that the products are close to each other in embedding space. To test this, we
train node2vec embeddings (Grover and Leskovec,, 2016) on the firm-product
graph, which is a static bipartite graph where an edge between a firm and a
product indicates that the firm supplies or buys the product. We use embedding
dimension 64 for the node2vec embeddings. Then, the score for product pair
$p_{1}$, $p_{2}$ is the cosine similarity of their node2vec embeddings.

##### Inventory module.

To learn production functions, we train our inventory module on the
transactions data (only the train set, for consistency with the edge
prediction experiments) and update it according to the inventory loss (6). In
practice, we find it effective to choose a debt penalty
$\lambda_{\textrm{debt}}$ that is 25% larger than the consumption reward
$\lambda_{\textrm{cons}}$, as reported in Table 4. We experiment with learning
the inventory module’s attention weights directly and with product embeddings,
following (4). We experiment with product embeddings from node2vec, from a
simple GNN trained on link prediction over the collapsed product-product graph
(collapsed from the firm-product graph, so an edge from product $p_{1}$ to
$p_{2}$ indicates how often a firm buys $p_{1}$ and supplies $p_{2}$), and
from SC-TGN. We find that product embeddings from node2vec help the inventory
module the most on synthetic data and embeddings from the simple GNN help it
the most on the microscope dataset. In Figure 6, we also show that our
inventory loss (6) is well-correlated with MAP, despite the inventory loss not
having access to any ground-truth parts.

Figure 6: Comparing inventory module’s loss (6) vs. MAP on ground-truth
production functions, trained on SS-std.

###  C.2 Edge prediction experiments

##### Model training.

| Synthetic data | Tesla | IED
---|---|---|---
SC-TGN
Memory dimension | 500 | 1000 | 1000
Embedding dimension | 500 | 1000 | 1000
Time dimension | 100 | 100 | 100
# neighbors for node embedding | 20 | 20 | 100
Update penalty $\lambda_{\textrm{update}}$ | 1 | 1 | 1
SC-GraphMixer
# MLPMixer layers | 2 | 2 | 2
Node encoding dimension | 500 | 50 | 300
Link encoding dimension | 100 | 10 | 10
# neighbors for node encoding | 20 | 100 | 10
# neighbors for link encoding | 20 | 10 | 2
Inventory module
Debt penalty $\lambda_{\textrm{debt}}$ | 5 | 5 | 5
Consumption reward $\lambda_{\textrm{cons}}$ | 4 | 4 | 4
Adjustment penalty $\lambda_{\textrm{adjust}}$ | 4 | 4 | 4
Training parameters
Batch size | 30 | 30 | 100 (SC-TGN), 30 (SC-GraphMixer)
Learning rate | 0.001 | 0.001 | 0.001
Max # epochs | 100 | 100 | 100
Patience | 10 | 10 | 10
Table 4: Hyperparameters that we used in our experiments. Figure 7:
Performance over 10 random seeds of SC-TGN on SS-std. Figure 8: Performance
over 10 random seeds of SC-GraphMixer on SS-std.

As discussed in the main text, for all datasets, we order transactions
chronologically and split them into train (70%), validation (15%), and test
(15%). In these experiments, we train the models on the train set, choose
hyperparameters based on the validation set, and report final results on the
test set. When training the models, we train for a maximum of 100 epochs, but
we allow for early stopping, if the validation MRR does not improve for over
10 epochs. In our main experiments (Table 2), we tested the following models:
Edgebank, Static, Graph Transformer, SC-TGN, SC-TGN+inv, SC-GraphMixer, and
SC-Graphmixer+inv. In Table 4, we report the hyperparameters we used for SC-
TGN and SC-GraphMixer, which we found by sweeping over the hyperparameters.
For all experiments, we ran 10 trials, each with a different random seed, and
reported the mean and standard deviation in test MRR and test RMSE. As shown
by the small standard deviations, performance is quite stable over random
seeds, which we also show in Figures 7-8. These figures also show that SC-TGN
and SC-Graphmixer both make steady progress in train and validation metrics
over epochs: train loss decreases, validation MRR increases, and validation
RMSE decreases. It is reassuring to see that the improvements in train loss
correspond well to improvements in validation MRR and RMSE, since the train
loss (16) combines multiple objectives, and, furthermore, we cannot directly
train on MRR, which is non-differentiable, so instead we use softmax cross-
entropy as a proxy loss.

|  Tesla |  IED
---|---|---
TGN |  0.612 (0.009) |  0.582 (0.016)
SC-TGN (id) |  0.537 (0.021) |  0.422 (0.015)
SC-TGN |  0.820 (0.007) |  0.842 (0.004)
Table 5: Ablations of SC-TGN: comparing to original TGN (Rossi et al.,, 2020)
and SC-TGN with ID embeddings, i.e., use memory directly as embedding, instead
of applying GNN to memories. We only evaluate edge existence here, with mean
reciprocal rank (MRR $\uparrow$), and leave out edge weight, since the
original TGN did not predict edge weight. We report mean and standard
deviation (in parentheses) over 10 seeds.

##### Ablations.

We can view the Static and Graph Transformer baselines as ablations of our
models, so we used the SC-TGN hyperparameters wherever applicable for fair
comparison (e.g., embedding dimension, batch size). Since SC-TGN learns the
initial memory, the Static baseline is equivalent to SC-TGN if no memory
updates are applied and the embedding module is ID. The Static baseline is
also a nested version of SC-GraphMixer, when we allow SC-GraphMixer to learn
its initial node features. The Graph Transformer baseline is equivalent to SC-
TGN if no memory updates are applied and the embedding module is UniMP (Shi et
al.,, 2021). Since SC-TGN was our best-performing model, we also ran further
ablations of SC-TGN to test the value of our extensions. We compared it to the
original TGN, where we stripped away most of the extensions discussed in
Appendix A.1: it no longer predicted edge weight, we removed the update
penalty and learnable initial memory, and we trained it only on perturbation
negatives. We also compared it to a version of SC-TGN that applied memory
updates, but used ID for the embedding module, i.e., uses the node memory
directly instead of applying a GNN to the memories. As shown in Table 5, we
found that SC-TGN greatly outperformed both of these ablations as well.

##### Analysis of +inv experiments.

For the +inv experiments, we found that validation performance was best if we
first pretrained the inventory module on only inventory loss (6), then co-
trained the inventory module and SC-TGN / SC-GraphMixer, following the joint
loss (16). Furthermore, in the +inv experiments, we experimented with whether
to include the inventory module’s edge prediction penalties or not (8)-(9).
Based on validation performance, we found that the penalties were useful on
synthetic data, in the standard setting and under supply shocks, but not when
there were missing transactions, or for the real-world datasets, so for those,
we turned off the penalties. We believe this is because, while the inventory
module is robust to missing transactions when learning production functions
(as shown by our production learning results in Table 1), the penalties are
not robust to missing transactions, since they rely on having an approximately
correct and complete inventory. For example, if product $p_{i}$ is an input
part to product $p_{o}$, and we are missing transactions where firm $f$
receives $p_{i}$, then the inventory module will incorrectly believe that $f$
has zero inventory of $p_{i}$ and penalize transactions where $f$ is supposed
to supply $p_{o}$.

|  SS-std |  SS-shocks |  SS-missing
---|---|---|---
SC-TGN |  0.522 (0.003) |  0.449 (0.004) |  0.494 (0.004)
SC-TGN+inv∗ |  0.548 (0.003) |  0.474 (0.003) |  0.476 (0.003)
SC-GraphMixer |  0.453 (0.005) |  0.426 (0.004) |  0.446 (0.003)
SC-GraphMixer+inv∗ |  0.477 (0.005) |  0.450 (0.005) |  0.430 (0.003)
Table 6: Testing the impact of inventory module on edge existence prediction,
when the inventory module is provided the ground-truth production functions.
We report mean and standard deviation (in parentheses) over 10 seeds.

To test this theory, we experimented with providing the inventory module with
the ground-truth production functions, and tested the usefulness of the edge
prediction penalties on the three synthetic datasets. We found that, on SS-std
and SS-shocks, the penalties help SC-TGN and SC-GraphMixer achieve stronger
performance on both edge existence and edge weight prediction, but for SS-
missing, performance worsens (Table 6). These results show that the penalties
hurting edge prediction performance are not indicative of how well the
inventory module has learned the production functions, since here we provided
the ground-truth production functions. Rather, missingness in the synthetic
data causes the penalties to hurt, not help, performance. This analysis
motivates the need to explore penalties in the future that are more robust to
missingness, as discussed in Appendix A.3, so that even with missing data, the
inventory module can not only learn production functions, but also aid edge
prediction so that the two objectives become more synergistic.