Retail Demand Forecasting with Graph Neural Networks: Mathematical Principles, Controlled Experiments, and Production Boundaries

At 10 a.m. on a Friday, a neighbourhood store of a supermarket chain is suddenly packed with shoppers. Watching a long line form in front of the chilled display cases, the store manager is baffled: today is not a public holiday, the store is running no discounts or promotions, and the backend single-store forecasting system recommended only normal everyday stock yesterday—where did this sudden rush come from?

The answer does not lie in this store’s own historical sales ledger. It lies two kilometers down the street in the same commercial district: a sister store down the block launched a major promotion this morning, but the eastern regional warehouse supplying that neighbour suffered a sorting breakdown the night before, leaving its shelves empty shortly after opening. Shoppers drawn into the district by the promotion came away empty-handed and walked down the street to this store, which receives normal replenishment from the western regional warehouse.

In our earlier posts—Retail Demand Forecasting with Prophet: From Model Mechanics to Engineering Practice and From a Single Tree to a Forest: Mathematical Principles, Algorithmic Evolution, and Engineering Practices of Tree Models—whether we decomposed trend and seasonality or engineered lag and calendar features, we relied on the same implicit premise: a store’s future is governed mainly by its own past. In a real retail chain, however, neighbouring stores constantly pull customer traffic back and forth in two opposite ways:

• Promotional cannibalization (also called traffic diversion or internal substitution): when a neighbouring store runs a deep discount with ample stock on hand, it draws away price-sensitive shoppers who would otherwise have visited our store;
• Stockout spillover (also called stockout substitution or demand overflow): when that neighbouring store sells out because an upstream warehouse failed to replenish it, disappointed shoppers turn around and flood into our store instead.

Once we step into a retail network like this, any forecasting method that looks only at a single store’s own history suffers from a structural information blind spot. Forecasting is no longer a matter of extrapolating an isolated time series; it must answer two intertwined questions at once: How do warehouse and store states evolve step by step over time, and how do supply and sales signals (known across retail and logistics as fulfilment) propagate across the warehouse–store network?

This is precisely where stochastic processes and graph neural networks (GNNs) meet: the former uses probability to describe how states evolve along the time axis, while the latter learns how signals relay across supply chains and neighbourhood street networks. But does upgrading our mathematical lens from an isolated series to a relational network automatically solve inventory forecasting in production? Rather than writing a triumphalist story in which graph models effortlessly crush traditional models, this article walks through a complete engineering arc—mechanism breakdown \(\to\) controlled verification \(\to\) production reality check:

1. Seeing the mechanism (Chapters 1–3): Starting from random state transitions over time, we progressively introduce entity-relationship graphs and map out the applicability boundaries of single-store time series, spatial graph propagation, and spatio-temporal graph models;
2. Controlled verification (Chapter 4): In a purpose-built synthetic network of “2 warehouses + 10 stores + staggered promotions,” we strip away model-complexity confounders to test under what exact conditions graph structure recovers signals that single-store models miss;
3. Confronting the production gap (Chapter 5): Stepping beyond paper benchmarks, we examine why a graph algorithm that beats every traditional reference model in theory can prove unexecutable in live retail systems once it collides with “on-book stock, empty-shelf” inventory distortion, data latency, single-store glitches cascading across the network, and strict operational red lines.

The companion synthetic data generator, multi-toolkit benchmarks, and plotting scripts are open-sourced at GitHub: retail-gnn-demand-forecast.

1. How Do Stochastic Processes Characterize Temporal State Evolution? #

Retail demand fluctuates across calendar time, and upstream warehouses switch among different operating states. Before introducing graph structure, Chapter 1 uses stochastic processes to establish the mathematical language for how states evolve over time; Chapter 2 then expands our view from a single entity to a relational network of entities.

1.1 Stochastic Processes, the Markov Property, and Transition Matrices #

A stochastic process is a family of random variables indexed by time or another parameter:

\[ \{X_t:t\in T\} \]

Here, \(T\) is the index set and \(X_t\) is the random variable associated with index \(t\). The index usually represents time, though it can also represent spatial location or another ordered dimension; \(X_t\) may be a discrete state, continuous value, or vector. The definition itself requires neither the Markov property nor a stationary distribution.

Adopting a discrete-time Markov chain further assumes that a well-designed current state already contains the historical information needed to predict the next step:

\[ \Pr(X_{t+1}\mid X_t,X_{t-1},\ldots,X_0)=\Pr(X_{t+1}\mid X_t) \]

This is the Markov property (often called memorylessness). It does not claim that earlier history has no business value; rather, it requires the model’s defined “current state” to already summarize the relevant history. If the next step still depends on how long a state has lasted or on its recent trajectory, that information should be incorporated into the state definition or modelled with a higher-order process.

For a finite state space \(\mathcal S=\{s_1,s_2,\ldots,s_m\}\), a row-stochastic matrix \(\mathbf P\) (also called a state transition matrix or Markov matrix) represents one-step state transitions. Because a system in state \(s_a\) today must transition to some state in \(\mathcal S\) tomorrow, every entry in \(\mathbf P\) is non-negative and each row sums strictly to \(1\) (\(100\%\)):

\[ P_{ab}=\Pr(X_{t+1}=s_b\mid X_t=s_a),\qquad \boldsymbol\pi_{t+1}=\boldsymbol\pi_t\mathbf P,\qquad \boldsymbol\pi_{t+k}=\boldsymbol\pi_t\mathbf P^k \]

Here, \(\boldsymbol\pi_t\) is the state probability distribution at time \(t\) (a row vector recording the probability of being in each state on day \(t\)). If the transition probabilities in \(\mathbf P\) do not change with calendar date \(t\), the Markov chain satisfies time homogeneity (or time invariance), and the matrix power \(\mathbf P^k\) extends one-step transitions into \(k\)-step distributional evolution. Real applications must estimate these transition probabilities from historical data and test whether time homogeneity holds.

1.2 Two Kinds of “Stationarity”: From Store Sales and Warehouse States to Catchment Foot Traffic #

Across time-series analysis and stochastic processes, we encounter two similarly named concepts that characterize retail system stability from different angles:

1. Weak Stationarity of a Time Series—Illustrated by Daily Store Sales

Weak stationarity (also called covariance stationarity or second-order stationarity) requires a series to maintain a constant mean and variance over calendar time, with autocovariance (how strongly sales on two days co-vary) depending only on the number of days between them rather than on the specific calendar date. In a real retail store, raw daily sales are almost never stationary: weekend traffic exceeds weekday traffic, summer and winter transitions create seasonal waves, and major promotions or public holidays produce sharp demand spikes. That is why single-store forecasting models (such as Prophet) explicitly introduce trend terms, seasonal components, and promotional covariates (commonly called external features in machine learning feature engineering, or exogenous variables / extra regressors in econometrics and time-series tools like Prophet)—stripping away calendar-driven non-stationary components first so that the remaining random fluctuations return to an approximately stationary, predictable track.

2. Stationary Distribution of a Markov Chain—From Warehouse Regimes to Cross-Store Foot Traffic

Rather than asking whether a series’ mean drifts across the calendar, a Markov chain’s stationary distribution (also called its equilibrium or steady-state distribution) asks: when the same transition operator \(\mathbf P\) is applied repeatedly, what long-run proportions does the system settle into? It is defined as a probability vector \(\boldsymbol\pi^*\) satisfying:

\[ \boldsymbol\pi^*=\boldsymbol\pi^*\mathbf P \]

When a finite Markov chain is irreducible and aperiodic (meaning any state can eventually reach any other state without getting trapped in an isolated pocket or locked into a rigid alternating cycle), repeatedly applying \(\mathbf P\) converges to the same \(\boldsymbol\pi^*\) regardless of the initial state. In our retail network, this single equation illuminates equilibrium behaviour along both the temporal and spatial axes:

• Temporal equilibrium (warehouse operating states): Suppose a regional warehouse occupies one of three daily states—normal supply, sorting congestion, or stockout disruption—switching from day to day according to a fixed transition matrix \(\mathbf P\). Even if the warehouse suffers a disruption today and recovers tomorrow, as long as the underlying transition mechanism remains unchanged, its long-run state proportions converge to \(\boldsymbol\pi^*\)—telling us what fraction of operating days the warehouse supplies normally, experiences congestion, or suffers stockouts over the long run. Those steady-state probabilities provide the theoretical basis for sizing regular safety stock in supply-chain planning.
• Spatial equilibrium and shock disruption (customer flow across stores): If we replace the states of a single warehouse with the multiple stores inside a shared commercial district, the transition matrix \(\mathbf P\) describes how shoppers circulate and substitute among neighbouring stores from day to day. During normal periods without external shocks, inter-store customer exchange stays balanced and store traffic shares hover near \(\boldsymbol\pi^*=\boldsymbol\pi^*\mathbf P\). The moment one store launches a deep promotion, however—or runs out of stock because of an upstream warehouse disruption—the district’s spatial equilibrium is broken on that day, forcing foot traffic to surge along the inter-store network as cannibalization or spillover. Conversely, if we multiply by the same transition matrix repeatedly without bound (\(\mathbf P^k,\,k\to\infty\)), no matter which store initially held the traffic, every store’s share is eventually flattened into the same steady-state proportion—keep this mathematical property in mind (that repeated multiplication washes out individual differences), as we will meet it again in Chapter 3 when discussing why graph models cannot be stacked arbitrarily deep.

As this comparison shows, a single Markov chain from Chapter 1 can track how one warehouse or one store evolves over time; capturing how an upstream warehouse disruption breaks the customer-flow equilibrium across neighbouring stores requires expanding our lens from isolated entities to a relational network—the starting point for Chapter 2.

2. How Graph Neural Networks Learn from Entity Relationships #

Chapter 1 described how a single entity’s state evolves over time. We now expand our view to a relational network of warehouses and stores. The same product may sit across dozens or hundreds of stores: a regional warehouse supplies multiple districts, and a promotion or stockout at one store can directly alter customer traffic at neighbouring stores.

Suppose we want to forecast tomorrow’s product demand at a specific store. A traditional single-store time-series model predicts from that store’s own sales history, prices, holidays, and promotions; yet if a nearby store runs a major promotion, or stocks out after an upstream warehouse delay, shoppers may switch stores. A single store’s own historical series cannot fully capture these cross-store shifts. The forecasting problem therefore changes: a system’s state depends not only on its own history, but also on connected entities.

2.1 Understanding Graph Structure Through a Minimal Retail Network #

To avoid switching scenarios across formulas, this chapter tracks one local network: warehouse \(W_1\) supplies store \(S_1\), warehouse \(W_2\) supplies store \(S_2\), and the two stores can influence each other’s customer traffic. Our task is to jointly predict next-day demand at \(S_1\) and \(S_2\).

A graph represents business entities as nodes and relationships between entities as edges. Write a graph as \(G=(\mathcal V,\mathcal E)\), where \(\mathcal V\) is the node set and \(\mathcal E\) is the edge set. This four-node teaching subgraph is:

\[ \begin{aligned} \mathcal V&=\{W_1,W_2,S_1,S_2\},\\ \mathcal E&=\{(W_1,S_1),(W_2,S_2),(S_2,S_1),(S_1,S_2)\}. \end{aligned} \]

Here, \(W_1,W_2\) are warehouse nodes and \(S_1,S_2\) are store nodes. Edges \((W_1,S_1)\) and \((W_2,S_2)\) represent supply relationships, while the spatial relationship between the two stores expands into two directed edges, \((S_1,S_2)\) and \((S_2,S_1)\). In compact arrow notation:

\[ W_1\overset{\mathrm{sup}}{\longrightarrow}S_1 \overset{\mathrm{spa}}{\longleftrightarrow}S_2 \overset{\mathrm{sup}}{\longleftarrow}W_2 \]

Above the arrows, \(\mathrm{sup}\) denotes a supply relationship and \(\mathrm{spa}\) denotes a spatial relationship. Arrow direction indicates that information can flow from the source node to the target node: supply information flows from warehouse to store, whereas the spatial relationship uses both directions so that either store’s information can inform the other’s forecast.

Because these edges distinguish source from target, this is a directed graph. For example, \(W_1\to S_1\) allows supply information to flow only from warehouse \(W_1\) to store \(S_1\), while \(S_2\to S_1\) allows the state of \(S_2\) to inform the forecast for \(S_1\). Direction specifies which entity’s information may update another entity; it does not represent a fixed direction of physical customer movement or predetermine whether demand increases or decreases.

Because this example includes two node types (warehouses and stores) and two edge types (supply and spatial), it forms a heterogeneous graph, also called a heterogeneous information network. “Heterogeneous” emphasizes multiple types of nodes or relationships; a graph with a single node type and a single edge type is homogeneous.

Distances and travel times in the diagram are edge features describing the relationships themselves. Assigning each edge a scalar \(w_e\) to represent connection strength yields a weighted graph. The two concepts should be kept distinct: an edge can carry multiple features such as distance and travel time, whereas “weight” usually refers specifically to the scalar strength used in aggregation or computation.

Finally, an edge in a graph only records a relationship that the model is allowed to use; it does not automatically prove causality. When \(S_2\) runs a promotion with ample stock, it may draw traffic away from \(S_1\); when \(S_2\) stocks out during a promotion, the same spatial connection may spill unmet demand over to \(S_1\). Graph topology answers which entities are linked; the direction and magnitude of a given day’s impact depend on the current states of the nodes and edges.

2.2 How Do We Represent Graph Structure and Business State as Data? #

The graph \(G\) in Section 2.1 only specifies which entities are connected by which relationships; it does not record demand, promotions, or supply availability on a given day. To pass this graph to a model, we must attach the business state of its nodes and edges on date \(t\). Because the graph is heterogeneous, warehouses and stores can use different node features, and supply and spatial edges can use different edge features, without padding everything with zeros into a single large matrix.

Using the same four-node network, we now construct a snapshot for date \(t\). Each row of a matrix corresponds to a node or an edge, and each column corresponds to a feature; the row order must remain fixed so that each value maps to the right business entity. All numbers below are author-designed teaching values and do not represent industry averages.

First split the node set into warehouses and stores:

\[ \mathcal V^{\mathrm W}=\{W_1,W_2\},\qquad \mathcal V^{\mathrm S}=\{S_1,S_2\}. \]

Warehouses retain only supply availability \(a_{W_i,t}^{\mathrm W}\in[0,1]\): 1 indicates normal outbound supply, 0 indicates no supply, and intermediate values indicate partial supply capacity. Stores use a demand signal \(b_{S_i,t}\) (the most recent demand already available when forecasting date \(t+1\), measured in units, excluding the realized outcome from date \(t+1\)) and a same-day promotion indicator \(p_{S_i,t}\in\{0,1\}\). With fixed row order, the two node-feature matrices can be written compactly as:

\[ \begin{aligned} \mathbf X_t^{\mathrm W} &= \begin{bmatrix} a_{W_1,t}^{\mathrm W}\\ a_{W_2,t}^{\mathrm W} \end{bmatrix} = \begin{bmatrix} 1\\ 0 \end{bmatrix} \quad \begin{array}{l} \leftarrow W_1\text{ (supply availability }a^{\mathrm W}\text{)}\\ \leftarrow W_2 \end{array} \\[6pt] \mathbf X_t^{\mathrm S} &= \begin{bmatrix} b_{S_1,t}&p_{S_1,t}\\ b_{S_2,t}&p_{S_2,t} \end{bmatrix} = \begin{bmatrix} 100&0\\ 80&1 \end{bmatrix} \quad \begin{array}{l} \leftarrow S_1\text{ (cols: demand }b\text{, promo }p\text{)}\\ \leftarrow S_2 \end{array} \end{aligned} \]

Thus, \(\mathbf X_t=\{\mathbf X_t^{\mathrm W},\mathbf X_t^{\mathrm S}\}\) represents the set of node features on date \(t\). Rather than one large zero-padded matrix, it is a typed collection consisting of a warehouse matrix and a store matrix.

Beyond node states, we also represent attributes of the relationships themselves. First partition the edges into supply and spatial relationships and fix the row order within each set:

\[ \mathcal E^{\mathrm{sup}}=\{W_1\to S_1,\,W_2\to S_2\},\qquad \mathcal E^{\mathrm{spa}}=\{S_1\to S_2,\,S_2\to S_1\}. \]

Which set an edge belongs to already identifies whether it is a supply or spatial relationship, so no extra one-hot type column (a binary 0/1 indicator column for each category) is needed. Letting \(d\) denote distance (km) and \(\tau\) denote travel time (min), the two edge-feature matrices are:

\[ \begin{aligned} \mathbf R_t^{\mathrm{sup}} &= \begin{bmatrix} d_{W_1\to S_1}&\tau_{W_1\to S_1,t}\\ d_{W_2\to S_2}&\tau_{W_2\to S_2,t} \end{bmatrix} = \begin{bmatrix} 10&20\\ 12&25 \end{bmatrix} \quad \begin{array}{l} \leftarrow W_1\to S_1\text{ (cols: distance km, time min)}\\ \leftarrow W_2\to S_2 \end{array} \\[6pt] \mathbf R_t^{\mathrm{spa}} &= \begin{bmatrix} d_{S_1\to S_2}&\tau_{S_1\to S_2,t}\\ d_{S_2\to S_1}&\tau_{S_2\to S_1,t} \end{bmatrix} = \begin{bmatrix} 2&6\\ 2&6 \end{bmatrix} \quad \begin{array}{l} \leftarrow S_1\to S_2\text{ (cols: distance km, time min)}\\ \leftarrow S_2\to S_1 \end{array} \end{aligned} \]

Distances typically remain constant while travel times may vary by date, so the edge-feature matrices retain the subscript \(t\), written as \(\mathbf R_t=\{\mathbf R_t^{\mathrm{sup}},\mathbf R_t^{\mathrm{spa}}\}\). To keep the hand calculation in Section 2.4 simple, we leave distance and travel time unused there; a production model can incorporate these edge features directly into its computation.

At this point, a graph snapshot on date \(t\) can be written as \((G,\mathbf X_t,\mathbf R_t)\): \(G\) records the topology, \(\mathbf X_t\) records warehouse and store states, and \(\mathbf R_t\) records supply and spatial edge attributes. The next section explains how a model reads this input to predict next-day demand for both stores.

2.3 How Do We Define and Train a GNN Forecasting Task? #

If an ordinary neural network reads only single-store features, it remains blind to warehouse states and neighbouring store shifts; unless engineers manually pre-build those relational signals into the input table, the model will not exchange information along the connections in \(G\). A graph neural network (GNN) differs precisely here: it incorporates the graph topology into its computation, passes signals along edges, and updates connected node states using shared, learnable rules.

First define the forecasting task cleanly. Denote the entire GNN by \(f_\Theta\), where \(\Theta\) represents parameters learned from training data. Given the date-\(t\) network snapshot \((G,\mathbf X_t,\mathbf R_t)\), the model jointly outputs next-day demand forecasts for both stores on date \(t+1\):

\[ \widehat{\mathbf y}_{t+1}= \begin{bmatrix} \hat y_{S_1,t+1}\\ \hat y_{S_2,t+1} \end{bmatrix}=f_\Theta\!\left(G,\{\mathbf X_t^{\mathrm W},\mathbf X_t^{\mathrm S}\},\{\mathbf R_t^{\mathrm{sup}},\mathbf R_t^{\mathrm{spa}}\}\right). \]

During training, each historical date forms a “same-day input to next-day label” pair. For compact notation, write the complete input on date \(t\) as \(\mathbf z_t\) and the two-store label vector constructed from realized observations on date \(t+1\) as \(\mathbf y_{t+1}^{\mathrm{target}}\):

\[ \mathbf z_t=(G,\mathbf X_t,\mathbf R_t),\qquad \mathbf y_{t+1}^{\mathrm{target}}= \begin{bmatrix} y_{S_1,t+1}^{\mathrm{target}}\\ y_{S_2,t+1}^{\mathrm{target}} \end{bmatrix}. \]

Let \(\mathcal T_{\mathrm{train}}\) be the set of all training input dates; the training set is then the collection of these input–label pairs:

\[ \mathcal D_{\mathrm{train}}=\left\{\left(\mathbf z_t,\mathbf y_{t+1}^{\mathrm{target}}\right):t\in\mathcal T_{\mathrm{train}}\right\}. \]

If we flatten every entry of the four matrices into an individual column, a single date’s training sample can be stored as one row of a wide table (a relational table that spreads all warehouse and store attributes for a given day horizontally across a single row with many columns). To fit the page width, the table below transposes that wide table: database fields appear as rows, and records for three dates appear as columns. Column headers indicate the input date \(t\), while the bottom two rows come from realized demand observations on the next day, \(t+1\).

Field	Business meaning	2025-01-01	2025-01-02	2025-01-03
\(a_{W_1}^{\mathrm W}\)	\(W_1\) supply availability	1	1	0
\(a_{W_2}^{\mathrm W}\)	\(W_2\) supply availability	0	1	1
\(p_{S_1}\)	\(S_1\) same-day promotion indicator	0	0	1
\(p_{S_2}\)	\(S_2\) same-day promotion indicator	1	0	0
\(d_{W_1\to S_1}\)	\(W_1\to S_1\) distance (km)	10	10	10
\(\tau_{W_1\to S_1}\)	\(W_1\to S_1\) travel time (min)	20	20	30
\(d_{W_2\to S_2}\)	\(W_2\to S_2\) distance (km)	12	12	12
\(\tau_{W_2\to S_2}\)	\(W_2\to S_2\) travel time (min)	25	20	22
\(d_{S_1\to S_2}\)	\(S_1\to S_2\) distance (km)	2	2	2
\(\tau_{S_1\to S_2}\)	\(S_1\to S_2\) travel time (min)	6	6	8
\(d_{S_2\to S_1}\)	\(S_2\to S_1\) distance (km)	2	2	2
\(\tau_{S_2\to S_1}\)	\(S_2\to S_1\) travel time (min)	6	6	8
\(b_{S_1}\)	\(S_1\) latest demand signal (units)	100	115	104
\(b_{S_2}\)	\(S_2\) latest demand signal (units)	80	52	83
\(y^{\mathrm{target}}_{S_1,t+1}\)	\(S_1\) next-day demand label (units)	115	104	72
\(y^{\mathrm{target}}_{S_2,t+1}\)	\(S_2\) next-day demand label (units)	52	83	86

Except for the bottom two rows, the symbols omit the shared date subscript \(t\); for instance, \(a_{W_1}^{\mathrm W}\) in the table corresponds to \(a_{W_1,t}^{\mathrm W}\) above, and \(\tau_{W_1\to S_1}\) corresponds to \(\tau_{W_1\to S_1,t}\). The first fourteen values in each column come from that day’s \(\mathbf X_t^{\mathrm W}\), \(\mathbf X_t^{\mathrm S}\), \(\mathbf R_t^{\mathrm{sup}}\), and \(\mathbf R_t^{\mathrm{spa}}\); the last two rows are the supervision targets for date \(t+1\). We place the two demand signals directly above the target labels to make it easy to compare inputs and next-day outcomes vertically. An actual database still stores one row per date; rows and columns are swapped here purely for readability.

The labels are not set arbitrarily from the 2025-01-01 inputs, but constructed from the next day’s realized data. Specifically, the 2025-01-01 column takes that day’s graph snapshot as input, while \((y^{\mathrm{target}}_{S_1,t+1},y^{\mathrm{target}}_{S_2,t+1})=(115,52)\) comes from realized demand at the two stores on 2025-01-02. Moving to the 2025-01-02 column, those same observations become that day’s latest demand signals \((b_{S_1},b_{S_2})=(115,52)\), while its labels \((104,83)\) come from 2025-01-03. Stepping one column to the right thus shows how the previous column’s \(t+1\) labels align with the next column’s date-\(t\) demand signals.

These numbers only illustrate temporal alignment and storage for training data; they do not mean any single field in the table directly determines the target. Because the training set stores one row for each input date in \(\mathcal T_{\mathrm{train}}=\{t_1,t_2,\ldots,t_T\}\), we have \(|\mathcal D_{\mathrm{train}}|=T\).

The fixed graph topology \(G\) does not need to be repeated in the wide table. It can be stored separately as a node table and an edge table, using stable node and edge IDs to identify which entity or connection each field belongs to. This wide table works well for explaining our fixed four-node toy graph; in real retail networks with many stores, warehouses, or dynamic edges, engineers typically split node features, edge features, and labels into long tables (also called narrow tables, stacking “date–entity ID–feature value” records vertically with a fixed number of columns and many rows) and assemble \(\mathbf X_t\), \(\mathbf R_t\), and \(\mathbf y_{t+1}^{\mathrm{target}}\) at training time.

Given these supervised examples, the GNN learns \(\Theta\) by minimizing prediction error across all training dates and target stores. Using mean squared error as an example:

\[ \Theta^*=\underset{\Theta}{\operatorname{argmin}}\; \frac{1}{|\mathcal D_{\mathrm{train}}|\,|\mathcal V^{\mathrm S}|} \sum_{t\in\mathcal T_{\mathrm{train}}} \sum_{S_i\in\mathcal V^{\mathrm S}} \left(y_{S_i,t+1}^{\mathrm{target}}-\hat y_{S_i,t+1}\right)^2 \]

Here, \(|\mathcal V^{\mathrm S}|=2\) is the number of target stores, so the denominator is the total number of “date–store” samples entering the loss. Training iteratively adjusts \(\Theta\) so that the prediction \(\hat y_{S_i,t+1}\) approaches the label \(y_{S_i,t+1}^{\mathrm{target}}\). Here we must distinguish between the recorded label on the books and the underlying business quantity we actually want to forecast: retail operations frequently suffer from sales censoring (known in statistics as censored data, and in retail as stockout clipping or constrained demand)—for instance, if a store has only 50 units left in stock that day, the checkout register can record at most 50 units sold even if 80 shoppers wanted to buy; the extra 30 units of true demand are sliced off and become completely invisible in sales reports. If we train directly on stockout-capped sales as our label, the model learns only to predict capped observed sales; if our target is unconstrained latent demand, we must first filter out or reconstruct the clipped stockout samples (a process called demand uncensoring) or explicitly model the inventory cap. A GNN does not automatically remove label bias on its own.

2.4 How Does Message Passing Turn Graph Data into a Forecast? #

With the inputs, outputs, and training objective defined, let us examine how \(f_\Theta\) executes a single forecast. The most common family of GNN architectures is the Message Passing Neural Network (MPNN): it generates messages based on each edge’s type and endpoint states, aggregates messages arriving at the same node, updates the node representations, and finally uses a prediction head to turn the updated representations into demand forecasts.

The Complete Process: From Graph Data to Demand Forecast #

Section 2.3 wrote the entire model as \(f_\Theta\) without unpacking its internal calculations. We now start from the model’s final output and trace backward: what the prediction head reads, and how each store’s final representation is formed.

A GNN maintains a node representation (also called a node embedding, hidden state, or node feature vector in deep learning) for every node, written in bold as \(\mathbf h\). It serves as a summary vector of what the model currently knows about that node. The node features \(\mathbf X_t\) on date \(t\) are first converted into initial representations \(\mathbf h^{(0)}\); after \(K\) rounds of message passing, target store \(S_i\) obtains its final representation \(\mathbf h_{S_i,t}^{(K)}\). The next-day forecasts for the two stores can first be written as[1]:

\[ \begin{aligned} \widehat{\mathbf y}_{t+1} &= \begin{bmatrix} \hat y_{S_1,t+1}\\ \hat y_{S_2,t+1} \end{bmatrix} = \begin{bmatrix} g_\omega\!\left(\mathbf h_{S_1,t}^{(K)}\right)\\ g_\omega\!\left(\mathbf h_{S_2,t}^{(K)}\right) \end{bmatrix}. \end{aligned} \]

Here, \(g_\omega\) is the prediction head (also called the output layer or readout layer, typically a simple linear regression or fully connected layer). It no longer passes messages across the graph; instead, it maps a single store’s final representation \(\mathbf h_{S_i,t}^{(K)}\) to a concrete next-day demand forecast \(\hat y_{S_i,t+1}\). The symbol \(g\) denotes the mapping from final representation to predicted value, and \(\omega\) denotes the parameters of that mapping learned during training. Both stores share the same \(g_\omega\), yet they can produce different outputs because their final representations differ.

Next, consider where \(\mathbf h_{S_i,t}^{(K)}\) comes from. The equations below use \(S_i\) for the target store receiving information, and \(W_j\) and \(S_j\) for the warehouse and neighbouring store sending information. Lowercase bold \(\mathbf r\) denotes a single edge’s feature vector; for instance, \(\mathbf r_{W_j\to S_i,t}^{\mathrm{sup}}\) is the row in \(\mathbf R_t^{\mathrm{sup}}\) corresponding to that supply edge.

Each propagation round consists of three consecutive steps:

1. Generate messages: A message function \(\phi\) can read the source node, target node, and edge features to generate a message along an edge. Because supply edges and spatial edges carry different business meanings, they use separate functions \(\phi^{\mathrm{sup}}\) and \(\phi^{\mathrm{spa}}\). A specific message function does not have to use every available argument.

2. Aggregate messages: \(\operatorname{AGG}\) combines all messages of the same relation type pointing to the same target store. Lowercase \(m\) denotes a single edge’s message, whereas bold \(\mathbf m\) denotes the aggregated message.

3. Update representations: An update function \(\psi\) combines the store’s current representation with the two aggregated messages to produce the next round’s representation.

In equation form, one round of computation is:

\[ \begin{aligned} \mathbf m_{S_i,t}^{(k),\mathrm{sup}} &= \operatorname{AGG}\left\{ \phi^{\mathrm{sup}}\!\left(\mathbf h_{W_j,t}^{(k)}, \mathbf r_{W_j\to S_i,t}^{\mathrm{sup}}\right): W_j\to S_i\in\mathcal E^{\mathrm{sup}}\right\},\\ \mathbf m_{S_i,t}^{(k),\mathrm{spa}} &= \operatorname{AGG}\left\{ \phi^{\mathrm{spa}}\!\left(\mathbf h_{S_j,t}^{(k)}, \mathbf h_{S_i,t}^{(k)}, \mathbf r_{S_j\to S_i,t}^{\mathrm{spa}}\right): S_j\to S_i\in\mathcal E^{\mathrm{spa}}\right\},\\ \mathbf h_{S_i,t}^{(k+1)} &= \psi\!\left(\mathbf h_{S_i,t}^{(k)}, \mathbf m_{S_i,t}^{(k),\mathrm{sup}}, \mathbf m_{S_i,t}^{(k),\mathrm{spa}}\right), \qquad k=0,1,\ldots,K-1. \end{aligned} \]

Setting \(k=0,1,\ldots,K-1\), the model repeats this calculation \(K\) times:

\[ \mathbf h_{S_i,t}^{(0)} \longrightarrow \mathbf h_{S_i,t}^{(1)} \longrightarrow \cdots \longrightarrow \mathbf h_{S_i,t}^{(K)}. \]

The superscript \((0)\) denotes the initial state before any messages are passed; when \(k=K-1\), the update equation yields \(\mathbf h_{S_i,t}^{(k+1)}=\mathbf h_{S_i,t}^{(K)}\), which the prediction head then turns into the final forecast. Here, \(K\) only counts message-passing rounds within the same date-\(t\) snapshot; it does not advance business time by \(K\) days.

A practical model learns the parameters of the message functions, update function, and prediction head from training data. Grouping the parameters of the two message functions as \(\theta_\phi\) and the update function parameters as \(\theta_\psi\), the complete parameter set is \(\Theta=\{\theta_\phi,\theta_\psi,\omega\}\). This links the internal mechanics of this section directly back to the overall GNN \(f_\Theta\) in Section 2.3.

Numerical Example: How Does the Four-Node Network Complete a Forecast? #

Let us return to the four-node retail network and walk through a complete forward pass after training. Training in Section 2.3 uses input–label pairs across many dates to learn \(\Theta\); here we treat the parameters as fixed and convert a new snapshot on date \(t\) into a forecast for date \(t+1\). To keep the arithmetic easy to follow by hand, the constants \(-30\), \(2/3\), \(0.9\), and \(5\) are given directly as teaching placeholders; they are not estimated from the three-day table above and should not be interpreted as empirical retail coefficients.

In a general model, the node representation \(\mathbf h_{S_i,t}^{(k)}\) is a vector. This example simplifies it to a scalar \(h_{S_i}^{(k)}\) and drops the shared date subscript \(t\). Both stores initialize their representations with their latest demand signals: \(h_{S_1}^{(0)}=100\) and \(h_{S_2}^{(0)}=80\). The same snapshot also provides \(a_{W_1,t}^{\mathrm W}=1\), \(a_{W_2,t}^{\mathrm W}=0\), \(p_{S_1,t}=0\), and \(p_{S_2,t}=1\), all taken directly from \(\mathbf X_t^{\mathrm W}\) and \(\mathbf X_t^{\mathrm S}\) in Section 2.2. Throughout the calculation, \(G\), \(\mathbf X_t\), and \(\mathbf R_t\) for date \(t\) remain fixed; the superscripts \((0)\), \((1)\), and \((2)\) track the initial state and two propagation rounds within this single snapshot, not a calendar progression to date \(t+2\).

In graph theory, traversing one edge is called one hop (connecting a node’s first-order neighbourhood), and traversing two consecutive edges forms a two-hop path (reaching its second-order neighbourhood). Two propagation rounds here correspond directly to a two-hop path: Round 1 lets each store receive supply status from its directly connected warehouse (one hop), and Round 2 lets the two stores exchange their updated states (second hop). Thus, \(W_2\) reaches \(S_2\) in one hop along \(W_2\to S_2\), and then influences \(S_1\) on the second hop along \(S_2\to S_1\). In standard GNNs, one propagation round typically corresponds to one network layer, though self-loops (edges from a node back to itself to retain its own state) and residual (or skip) connections (adding the previous layer’s representation directly to the updated output so local features are not washed out) alter what information is retained, so layer count should not be equated mechanically with effective propagation distance.

The table below simplifies the general functions into four hand-calculable rules. A general model would first encode warehouse features into a warehouse representation \(\mathbf h_W\); this example skips that step and lets the supply message read warehouse availability \(a^{\mathrm W}\) directly. The spatial message reads source-store demand \(b_{S_j,t}\), promotion indicator \(p_{S_j,t}\), and Round-1 representation \(h_{S_j}^{(1)}\); distance and travel time remain stored in \(\mathbf R_t\) but do not enter this hand calculation.

Component	Inputs	Simplified rule
Supply message \(\phi^{\mathrm{sup}}\)	\(a_{W_j,t}^{\mathrm W}\)	\(m_{W_j\to S_i,t}^{(0),\mathrm{sup}}=-30(1-a_{W_j,t}^{\mathrm W})\)
Spatial message \(\phi^{\mathrm{spa}}\)	\(b_{S_j,t},\,p_{S_j,t},\,h_{S_j}^{(1)}\)	\(m_{S_j\to S_i,t}^{(1),\mathrm{spa}}=p_{S_j,t}\frac{2}{3}\max(0,b_{S_j,t}-h_{S_j}^{(1)})\)
Update function \(\psi\)	\(h_{\mathrm{old}},\,m\)	\(h_{\mathrm{new}}=h_{\mathrm{old}}+m\)
Prediction head \(g_\omega\)	\(h_{S_i}^{(2)}\)	\(\hat y=g_\omega(h)=0.9h+5\)

Message subscripts indicate which directed edge carries the message; superscripts \((0)\) and \((1)\) indicate whether the message is computed from the initial state or the Round-1 state; and \(\mathrm{sup}\) versus \(\mathrm{spa}\) distinguishes supply from spatial messages. In this prediction head, \(\omega\) consists of the slope \(0.9\) and intercept \(5\). Production models use richer parameterized functions; these four rules simply make one forward pass transparent.

We can write the complete inference pass upfront as the two-store forecast vector:

\[ \begin{aligned} \widehat{\mathbf y}_{t+1} &= \begin{bmatrix} \hat y_{S_1,t+1}\\ \hat y_{S_2,t+1} \end{bmatrix} = \begin{bmatrix} 0.9h_{S_1}^{(2)}+5\\ 0.9h_{S_2}^{(2)}+5 \end{bmatrix} = \begin{bmatrix} 0.9\times120+5\\ 0.9\times50+5 \end{bmatrix} = \begin{bmatrix} 113\\ 50 \end{bmatrix}. \end{aligned} \]

The values \(120\) and \(50\) come from the two propagation rounds below.

Round 1: Supply messages update store representations. When \(a^{\mathrm W}=1\), the warehouse supplies normally and sends a message of 0; when \(a^{\mathrm W}=0\), the message is \(-30\). Consequently, \(W_1\) leaves the representation of \(S_1\) unchanged, whereas the disrupted warehouse \(W_2\) reduces \(S_2\) from 80 to 50:

\[ \begin{aligned} m_{W_1\to S_1,t}^{(0),\mathrm{sup}} &= -30(1-1)=0,\\ m_{W_2\to S_2,t}^{(0),\mathrm{sup}} &= -30(1-0)=-30,\\[4pt] h_{S_1}^{(1)} &= 100+m_{W_1\to S_1,t}^{(0),\mathrm{sup}}=100,\\ h_{S_2}^{(1)} &= 80+m_{W_2\to S_2,t}^{(0),\mathrm{sup}}=50. \end{aligned} \]

Here, \(-30\) is simply a teaching value showing that supply constraints lower a store’s internal state; it does not mean realized demand necessarily drops by 30 units.

Round 2: Spatial messages transmit the supply impact to the neighbouring store. The spatial function first computes the source store’s demand shortfall \(b_{S_j,t}-h_{S_j}^{(1)}\), and then uses \(p_{S_j,t}\) to gate transmission: when \(p=0\), the entire message is 0; when \(p=1\), this rule passes two-thirds of the shortfall to the neighbour. Because promoted store \(S_2\) has a shortfall of \(80-50=30\) units, it sends 20 to \(S_1\); unpromoted store \(S_1\) sends 0 to \(S_2\):

\[ \begin{aligned} m_{S_2\to S_1,t}^{(1),\mathrm{spa}} &= 1\times\frac{2}{3}(80-50)=20,\\ m_{S_1\to S_2,t}^{(1),\mathrm{spa}} &= 0\times\frac{2}{3}\max(0,100-100)=0,\\[4pt] h_{S_1}^{(2)} &= 100+m_{S_2\to S_1,t}^{(1),\mathrm{spa}}=120,\\ h_{S_2}^{(2)} &= 50+m_{S_1\to S_2,t}^{(1),\mathrm{spa}}=50. \end{aligned} \]

Finally, the prediction head reads \(h_{S_1}^{(2)}=120\) and \(h_{S_2}^{(2)}=50\) to output predicted values \(\hat y_{S_1,t+1}=113\) and \(\hat y_{S_2,t+1}=50\) (which sit close to the realized next-day labels \(115\) and \(52\) in the 2025-01-01 column of Section 2.3, leaving small prediction residuals of \(2\) units each). The spatial message of 20 is not the final forecast for \(S_1\); it first updates the internal representation of \(S_1\) from 100 to 120, which then alters the forecast through \(g_\omega\). This number represents an internal demand-spillover signal generated by the message function, not a guarantee that exactly 20 shoppers or 20 orders physically switched stores.

What Characteristics of a GNN Does This Example Illustrate? #

The point of this example is not merely that “the model sees more columns,” but how those columns participate in computation according to the graph topology. It highlights three core characteristics of a GNN:

• The graph determines where information can flow: Supply status at \(W_2\) first enters \(S_2\) along a supply edge and then influences \(S_1\) along a spatial edge. If either directed edge \(W_2\to S_2\) or \(S_2\to S_1\) were absent, this two-hop path would not exist. Here the graph acts as a computational constraint, not just a metadata label.

• Multi-round propagation composes multiple relationship types: Round 1 writes warehouse supply availability into the store representation, and Round 2 computes cross-store messages from store demand, promotion status, and that updated representation. As a result, the forecast for \(S_1\) reflects not only its own initial 100, but also a two-hop signal formed jointly by “upstream supply constraint + neighbour promotion.” Rather than simply placing two raw columns side by side, later rounds compute over the updated states from earlier rounds.

• Edges of the same relation type share computation rules: Both \(W_1\to S_1\) and \(W_2\to S_2\) share the supply message function, while \(S_1\to S_2\) and \(S_2\to S_1\) share the spatial message function. When the network scales to more warehouses and stores, the model reuses the rules learned for each relation type without needing a separate parameter set for every node.

These properties do not mean an ordinary neural network is mathematically incapable of using warehouse or neighbour information. Engineers can flatten neighbourhood fields into a wide table and feed them to a multilayer perceptron (MLP, a standard fully connected neural network), linear model, or tree model. What distinguishes a GNN is that it embeds “who can pass which type of information to whom” directly into the model architecture and reuses the same relational operations across nodes and edges. Whether that architecture beats a well-engineered traditional tabular model (often called a baseline model or reference model in machine learning—the foundational benchmark used to test whether a complex new model delivers genuine incremental value) still requires out-of-time validation (commonly called temporal backtesting in forecasting engineering—training strictly on historical periods and evaluating on future holdout periods without ever shuffling across time).

To highlight the two-hop path \(W_2\to S_2\to S_1\), this hand calculation processes only supply edges in Round 1 and only spatial edges in Round 2; a general heterogeneous GNN can generate both message types in every round and aggregate them separately. The diagram below places the node features, two propagation rounds, and prediction head along a single computational chain.

2.5 How Do GCN, GraphSAGE, and GAT Implement Different Neighbour Aggregations? #

The MPNN formulation in Section 2.4 specifies the general “generate message—aggregate messages—update representation” pipeline, without fixing which neighbours must be read in each round or how much weight each neighbour should receive. Graph Convolutional Networks (GCN), Graph Sample and Aggregate (GraphSAGE), and Graph Attention Networks (GAT) are three classic variants of this framework; they differ primarily in which neighbours they select and how they weight those neighbours during aggregation.

Before comparing aggregation rules, however, we must address node heterogeneity. Classical GCN, GraphSAGE, and GAT equations typically assume that all node representations entering a given layer already live in a shared hidden space and support compatible matrix operations. Our retail network does not satisfy that assumption out of the box: warehouse nodes carry only supply availability, whereas store nodes carry demand and promotion indicators, so their raw feature dimensions and semantics differ. Practical models therefore encode each node type separately and project them into a common hidden dimension first:

\[ \begin{aligned} \mathbf H_{\mathrm W}^{(0)}=f_{\mathrm W}(\mathbf X_t^{\mathrm W}), \qquad \mathbf H_{\mathrm S}^{(0)}=f_{\mathrm S}(\mathbf X_t^{\mathrm S}), \\[4pt] \mathbf H_{\mathrm W}^{(0)}\in \mathbb R^{|\mathcal V^{\mathrm W}|\times d_h}, \qquad \mathbf H_{\mathrm S}^{(0)}\in \mathbb R^{|\mathcal V^{\mathrm S}|\times d_h}. \end{aligned} \]

Here, \(f_{\mathrm W}\) and \(f_{\mathrm S}\) are type-specific input encoders for warehouses and stores, and \(d_h\) is the shared hidden dimension. Only after this projection can warehouse and store representations be combined inside subsequent message functions. Edge heterogeneity is a second issue: even when connected representations share the same hidden dimension, supply edges and spatial edges should not necessarily share the same message parameters. The classical formulas below illustrate core aggregation mechanisms; they do not automatically handle the two node types and two edge types in our retail graph on their own.

GCN (Graph Convolutional Network): Normalizing by graph topology and node degree. The classical GCN proposed by Kipf and Welling typically reads the full one-hop neighbourhood and mixes self and neighbour representations using a normalized adjacency matrix[2]. The formulation below assumes all rows already share the same hidden dimension and temporarily treats all connections as a single relation type:

\[ \mathbf H^{(k+1)}= \sigma\!\left( \widehat{\mathbf D}^{-1/2} \widehat{\mathbf A} \widehat{\mathbf D}^{-1/2} \mathbf H^{(k)} \mathbf W^{(k)} \right), \qquad \widehat{\mathbf A}=\mathbf A+\mathbf I. \]

Here, each row of \(\mathbf H^{(k)}\) is a store or warehouse representation at layer \(k\); \(\mathbf A\) records connections, \(\mathbf I\) adds self-loops, and \(\widehat{\mathbf D}\) normalizes by node degree (the number of edges connected to each node); \(\mathbf W^{(k)}\) and the nonlinear activation \(\sigma\) apply a learnable feature transformation. For two connected stores with scalar features 100 and 60, setting \(W^{(k)}=1\) and \(\sigma\) to the identity yields the simplified self-loop aggregation:

\[ \frac{1}{2}\times100+\frac{1}{2}\times60=80 \]

In a GCN, aggregation weights are determined primarily by graph connectivity and node degrees; during a forward pass at a given layer, the model does not automatically increase a neighbour’s structural weight just because that neighbour is running a promotion that day. Its simplicity makes it a natural graph baseline, but fixed structural normalization can be too rigid when neighbour influence shifts with business state.

GraphSAGE (Graph Sample and Aggregate): Sampling neighbours before applying a shared aggregator. Proposed by Hamilton et al., GraphSAGE randomly samples a fixed number of neighbours from each node’s neighbourhood and combines them using a mean, pooling (such as element-wise max), or other learnable aggregator[3]. The classical form below assumes the sampled neighbour representations are mutually comparable and share a single aggregator:

\[ \begin{aligned} \widetilde{\mathbf h}_{S_i}^{(k)} &= \operatorname{AGG}\!\left( \{\mathbf h_{S_j}^{(k)}:S_j\in\mathcal N(S_i)\} \right),\\ \mathbf h_{S_i}^{(k+1)} &= \sigma\!\left( \mathbf W^{(k)} [\mathbf h_{S_i}^{(k)}\,\|\,\widetilde{\mathbf h}_{S_i}^{(k)}] \right). \end{aligned} \]

Here, \(\mathcal N(S_i)\) denotes the neighbour subset sampled in the current round, and \(\|\) denotes vector concatenation (joining the target store’s own representation and its neighbourhood summary end to end). Sampling keeps per-node computation bounded on large graphs, while the shared aggregation and transformation functions can generate representations directly for newly opened stores that were never seen during training—a capability known in machine learning as inductive learning (learning generalizable aggregation rules that apply to unseen nodes, as opposed to transductive learning, which requires all test nodes to be present in the graph during training). This ability to handle unseen nodes still comes with conditions: a new store must have the required input features and a valid neighbourhood, and transfer can degrade if the new region’s distribution differs sharply from the training set. On a heterogeneous graph, engineers must also decide whether to sample separately by node and relation type or share a common sampling budget and aggregator.

GAT (Graph Attention Network): Learning neighbour weights from current features. Proposed by Veličković et al., GAT brings the attention mechanism into graph learning, dynamically computing attention scores across a target node’s neighbours and normalizing them (via Softmax) within that neighbourhood[4]. The classical form below assumes those neighbours already live in a comparable representation space and computes attention within a single relation type:

\[ \alpha_{S_iS_j}= \frac{\exp(e_{S_iS_j})} {\sum_{S_\ell\in\mathcal N(S_i)}\exp(e_{S_iS_\ell})} \]

Here, \(e_{S_iS_j}\) is a relevance score computed by the model from the current representations of target store \(S_i\) and neighbour \(S_j\). Suppose two neighbours of \(S_1\) receive scores \(e_{S_1S_2}=\ln 3\) and \(e_{S_1S_3}=0\); their attention weights are:

\[ \alpha_{S_1S_2}=\frac{3}{3+1}=0.75, \qquad \alpha_{S_1S_3}=\frac{1}{3+1}=0.25 \]

If the two neighbours pass scalar representations of 80 and 120, the aggregated result is \(0.75\times80+0.25\times120=90\). Unlike fixed degree normalization, GAT allows neighbour weights to vary dynamically with current features (such as whether a neighbour is on promotion today). Even so, attention weights only describe how the model allocates weight during that aggregation step; they do not prove that a specific neighbour causally drove the demand change. If both supply and spatial edges participate in attention on a heterogeneous graph, relation-specific projections or scoring functions are required so that scores across different relation types remain comparable.

Comparing the three side by side:

Model	Neighbours used per round	How the neighbourhood summary is formed	Primary practical consideration
GCN	Typically the full one-hop neighbourhood	Normalized by graph topology and node degree	Simple architecture; well suited as a graph baseline
GraphSAGE	Can sample a fixed number of neighbours	Shared mean, pooling, or learnable aggregator	Controls compute on large graphs and supports inductive generalization to unseen nodes
GAT	Typically the target node’s current neighbourhood	Dynamically learns attention weights from current representations	More flexible when neighbour influence varies with input state

All three follow the same “neighbours send messages—target node aggregates—node representation updates” pattern; they do not form a hierarchy from weak to strong performance. Choosing among them depends on graph scale, whether unseen nodes must be handled, whether neighbour effects vary with state, and out-of-time validation results. Note too that the classical equations above illustrate aggregation over a single relation type, and classical GCN is often applied to undirected or symmetrized graphs; for the heterogeneous directed retail graph in this article, an actual implementation must still preserve edge direction, distinguish supply from spatial edges, and assign relation-specific message or transformation functions as in Section 2.4.

2.6 Graph Neural Networks Are Not the Best Choice for Every Relational Forecasting Problem #

Understanding how a GNN exploits relationships does not mean every store network should be modelled with a GNN. An AI engineer’s goal is to deliver reliable forecasts sustainably at an acceptable cost, not to pick the most complex-looking architecture. For most retail demand forecasting projects, a sounder starting point is to translate business relationships into verifiable features and build a solid, cost-effective baseline demand forecasting model first.

For example, each “store–date” sample can be augmented with:

• sales at nearby stores over the past seven days;
• distance-weighted promotion intensity across neighbouring stores;
• inventory availability at the store’s upstream warehouse;
• delivery distance between the warehouse and the store.

These network features can enter a linear regression model or serve as exogenous regressors in Prophet; when interactions are distinctly nonlinear, Gradient Boosted Decision Trees (GBDT, most commonly implemented with XGBoost or LightGBM) are natural choices. Well-engineered traditional models often deliver an excellent cost–performance ratio: they train and infer faster, require less data, make error attribution easier to trace, and are simpler to explain to business stakeholders. A solid baseline demand forecasting model is not a perfunctory step before running GNN experiments; it is the yardstick that measures how much incremental value graph structure actually provides.

Model selection also cannot rely solely on error metrics from a single offline split. Production systems must account for robustness (a system’s ability to remain stable under noisy inputs or changing environments) and Machine Learning Operations (MLOps) return on investment (ROI): Does the model remain stable when edges drop out, warehouses or stores open and close, or relationship weights drift? Can the team monitor shifts in node features, edge features, and graph topology? When an anomalous forecast appears, can engineers tell whether the issue stems from raw data, graph construction, message passing, or the prediction head? If a GNN yields only marginal or unstable accuracy gains while significantly increasing data-pipeline, training, deployment, monitoring, and maintenance costs, its overall ROI may fall below that of a simpler model.

That does not negate the value of GNNs. When the forecasting target genuinely depends on multi-hop relationships, heterogeneous entities, and state-dependent neighbour interactions that are difficult to flatten stably into a fixed-width table, a GNN can directly reuse local computation rules on the graph and propagate information along supply and spatial edges. It can also share relational parameters across stores and continue handling variable neighbour counts as the network grows. These architectural assumptions baked directly into the network structure (known in machine learning as structural inductive biases—the built-in prior that connected entities interact through shared relational rules) are difficult for linear regression, Prophet, or gradient-boosted trees to express naturally through fixed tabular columns alone.

A pragmatic engineering sequence is therefore: build interpretable single-store and feature-engineered baselines first, evaluate the GNN’s incremental gain under strict out-of-time validation, and track latency, resource consumption, failure modes, and maintenance overhead alongside accuracy. Only when the gain from relational structure consistently outweighs the added complexity does a GNN merit production deployment. The next chapter places node states, relational propagation, and stochastic disturbances into a unified time-varying retail network framework to clarify where that advantage can arise.

3. From Stochastic Processes to Spatio-Temporal Graph Neural Networks: How Do States and Relationships Co-evolve? #

Chapter 1 introduced how a single entity’s state evolves over time, and Chapter 2 showed how information propagates across a graph within a single dated snapshot. Real retail demand spans both dimensions at once: predicting tomorrow requires knowing not only which warehouses and stores are connected today, but also how warehouse states, promotions, and store demand evolved from the past to the present. Spatio-temporal forecasting must therefore operate along two axes: reading consecutive snapshots across time, and exchanging information across the graph within each snapshot.

3.1 From Univariate Time Series to Dynamic Networks #

Consider a single store first. Its state on date \(t+1\) may depend jointly on its current state, known external inputs, and unexplained disturbances:

\[ x_{t+1}=f(x_t,u_t)+\epsilon_t \]

Here, \(x_t\) is the demand state or related features on date \(t\), \(u_t\) represents external inputs such as promotions, prices, and holidays, and \(\epsilon_t\) represents stochastic disturbances not explained by the model. This equation abstracts how the underlying business process evolves; it is not the training objective of a specific forecasting model.

When the forecasting scope expands to multiple stores and warehouses, each date is no longer represented by a single scalar \(x_t\), but by a graph snapshot as defined in Chapter 2:

\[ \mathcal Z_t=(G_t,\mathbf X_t,\mathbf R_t) \]

Here, \(G_t\) records the nodes and edges on date \(t\), \(\mathbf X_t\) stores the node features for that day, and \(\mathbf R_t\) stores the edge features for each relation type. If warehouse supply assignments and store locations remain unchanged over the short term, \(G_t\) can stay fixed while demand, promotions, supply availability, and travel times still cause \(\mathbf X_t\) and \(\mathbf R_t\) to vary from day to day.

What the forecasting model actually reads is a window of historical snapshots to output next-day demand across all stores:

\[ \widehat{\mathbf y}_{t+1} = F_\theta\!\left( \mathcal Z_{t-L+1},\ldots,\mathcal Z_t \right) \]

Here, \(L\) is the lookback window length, \(F_\theta\) is the prediction function learned from training data, and \(\widehat{\mathbf y}_{t+1}\) jointly contains the next-day forecasts for all stores. This formulation cleanly separates “how the system evolves” from “how the model uses history to predict,” without restricting \(F_\theta\) to any particular GNN or time-series architecture.

3.2 Why Can Graph Structure Help Time-Series Forecasting? #

Return to the four-node network from Chapter 2 (\(W_1\to S_1 \leftrightarrow S_2 \leftarrow W_2\)). Even when the graph connections stay fixed for days, the business signals carried along those edges can undergo a sign flip as upstream warehouse states shift:

• Date \(t\) (negative cannibalization): Upstream warehouse \(W_2\) supplies normally (\(a_{W_2,t}^{\mathrm W}=1\)) and neighbouring store \(S_2\) runs a promotion (\(p_{S_2,t}=1\)); well-stocked \(S_2\) draws price-sensitive shoppers away from \(S_1\), pushing \(S_1\)’s demand down.

• Date \(t+1\) (positive spillover): \(W_2\) suffers a stockout disruption (\(a_{W_2,t+1}^{\mathrm W}=0\)) while \(S_2\) continues its promotion (\(p_{S_2,t+1}=1\)); \(S_2\) quickly sells out, and disappointed shoppers flow along the very same spatial edge to well-stocked \(S_1\), pushing \(S_1\)’s demand up.

Notice the trap here: from the standpoint of focal store \(S_1\), observing the one-hop neighbour promotion flag \(p_{S_2}=1\) alone cannot tell the model whether the effect on \(S_1\) is negative or positive—its sign depends entirely on whether the two-hop upstream warehouse \(W_2\) is supplying normally. A single-store model looking only at \(S_1\)’s history is blind to both; even if an engineer manually appends one-hop neighbour promotions into a flat feature table, failing to condition on the two-hop upstream warehouse state causes positive and negative episodes to cancel out across training samples. A two-layer GNN, by contrast, follows the directed heterogeneous path \(W_2 \to S_2 \to S_1\), first combining “warehouse state \(\times\) store promotion” at \(S_2\) and then passing that supply-conditioned signal to \(S_1\).

Two indices must be kept distinct here. The subscripts \(t,t+1\) mark the progression of business calendar dates; the superscripts \(k,k+1\) mark another round of message passing within the same date snapshot. Running two rounds of graph propagation does not turn date \(t\) into date \(t+1\); it simply extends a store’s field of view on date \(t\) from “one-hop warehouses and neighbours” to “two-hop upstream warehouses of neighbours.” A spatio-temporal graph model must handle both dimensions before producing a next-day forecast.

3.3 Random Walks on Graphs, Feature Diffusion, and Spatio-Temporal Forecasting #

Chapter 1 used a row-stochastic matrix \(\mathbf P\) to describe how a single warehouse’s state transitions over time; the same matrix form can also characterize how features propagate across a store spatial graph within a single snapshot. Turning spatial proximity into a computable propagation operator involves three matrices:

• Spatial adjacency matrix \(\mathbf A\) (who is near whom, and how close): Non-negative entries \(A_{ij}\ge 0\) measure the raw spatial affinity between store \(i\) and store \(j\) (larger for closer stores, and \(0\) outside the local catchment area). Because \(A_{ij}\) is only a raw proximity score, each row need not sum to \(1\).
• Row-degree matrix \(\mathbf D=\operatorname{diag}(\mathbf A\mathbf 1)\) (each store’s total surrounding connectivity): The diagonal entry \(D_{ii}=\sum_j A_{ij}\) sums the raw affinities between store \(i\) and all of its neighbouring stores. A store inside a dense commercial district has a large \(D_{ii}\); an isolated store has a small one.
• Random-walk transition matrix \(\mathbf P=\mathbf D^{-1}\mathbf A\) (within-cluster allocation and aggregation weights): When every node has at least one outgoing edge, left-multiplying by \(\mathbf D^{-1}\) divides each row of \(\mathbf A\) by its row sum \(D_{ii}\), ensuring that every row sums strictly to \(1\):

\[ \mathbf P=\mathbf D^{-1}\mathbf A,\qquad P_{ij}=\frac{A_{ij}}{D_{ii}}. \]

We can extend the two-store example (\(S_1, S_2\)) from Chapter 2 by adding a third store \(S_3\) in the same district to see how these three matrices work together. Focusing on store \(S_1\), suppose it sits very close to \(S_2\) (raw affinity \(A_{12}=0.6\)) and slightly further from \(S_3\) (raw affinity \(A_{13}=0.2\)), with no self-loop (\(A_{11}=0\)). The first row of the raw adjacency matrix is \(\begin{bmatrix}0 & 0.6 & 0.2\end{bmatrix}\), and its row degree (total connectivity) is \(D_{11}=0.6+0.2=0.8\). Dividing by \(D_{11}\) yields the first row of the transition matrix \(\mathbf P\):

\[ P_{1\cdot} = \frac{1}{0.8}\begin{bmatrix}0 & 0.6 & 0.2\end{bmatrix} = \begin{bmatrix}0 & 0.75 & 0.25\end{bmatrix}. \]

In retail terms, this row normalization has a clear interpretation: among all neighbours around \(S_1\)—whether we are pulling neighbour promotion signals inward to \(S_1\) or estimating where unmet customer demand flows outward when \(S_1\) stocks out—adjacent store \(S_2\) carries 75% of the relational weight and slightly further store \(S_3\) carries 25%.

If each row of \(\mathbf H^{(0)}\) stores a store’s initial features, the simplest linear graph propagation repeatedly left-multiplies the feature matrix by \(\mathbf P\):

\[ \mathbf H^{(k+1)}=\mathbf P\mathbf H^{(k)},\qquad \mathbf H^{(K)}=\mathbf P^K\mathbf H^{(0)} \]

Under row normalization, \(P_{ij}\) is the one-step weight pointing from store \(i\) to neighbouring store \(j\); each matrix multiplication mixes neighbour features into the focal store according to those transition shares, \(\mathbf P^2\) composes two-step paths, and \(\mathbf P^K\) allows information to reach up to \(K\) hops away. In our three-store example above, if \(S_2\) and \(S_3\) have promotion depths of \(0.20\) and \(0.00\) on day \(t\), one step of propagation \(\mathbf P\mathbf H^{(0)}\) aggregates an effective neighbourhood promotion intensity of \(0.75\times 0.20 + 0.25\times 0.00 = 0.15\) for \(S_1\). This step count connects directly to GNN propagation rounds: one aggregation round reads one-hop neighbours, and multiple rounds progressively expand the receptive field (a term borrowed from computer vision denoting how many hops out in the network a focal store can “see” when making a prediction), while a GNN additionally applies learnable transformations, nonlinearities, residual connections, and relation-specific rules.

This formulation also explains why repeated propagation gradually weakens differences across nodes: in the simplified setting of a fixed linear propagation operator, as \(K\) grows, the matrix power \(\mathbf P^K\) can approach the stationary limit discussed in Section 1.2, making node features increasingly similar across the graph. That is why deep GNNs are prone to over-smoothing—excessive neighbourhood averaging washes out local store differences.

Putting the four-node toy network from Section 3.2 together with the matrix operators here clarifies the division of labour between the two edge types in our heterogeneous retail network: warehouse\(\to\)store supply edges \(\mathcal E^{\mathrm{sup}}\) transmit upstream inventory availability across layers (determining whether a store faces replenishment caps the next day), whereas store\(\leftrightarrow\)store spatial edges \(\mathcal E^{\mathrm{spa}}\) use the row-normalized matrix \(\mathbf P^{\mathrm{spa}}\) to govern customer allocation and feature aggregation among peer stores in the same layer. Note also the two distinct business meanings carried by the same row-stochastic matrix form: in the warehouse Markov chain \(\mathbf P^{\mathrm{wh}}\), indices represent operating states of a single warehouse and entries give transition probabilities among Normal, Congested, and Disrupted over time; in the store spatial graph \(\mathbf P^{\mathrm{spa}}\), indices represent distinct stores and \(P_{ij}\) serves as a within-snapshot catchment allocation weight. Chapter 4’s synthetic experiment reads the same \(\mathbf P^{\mathrm{spa}}\) in both directions: when aggregating promotional competition from neighbours toward focal store \(i\), we center on the focal store and read row entries \(P_{ij}\); when allocating unmet demand outward from stocked-out store \(j\) to store \(i\), we center on the source store and read \(P_{ji}\).

In continuous-time formulations, replacing discrete propagation steps \(K\) with continuous diffusion time \(\tau\) and a graph Laplacian matrix \(\mathbf L = \mathbf D - \mathbf A\) (an operator that measures how much each node’s value differs from its surrounding neighbours) yields the graph heat equation \(\frac{\mathrm d\mathbf H(\tau)}{\mathrm d\tau}=-\mathbf L\mathbf H(\tau)\) (analogous to heat diffusing naturally along pipes from warmer nodes to cooler neighbours). Classical Spatio-Temporal Graph Neural Network (ST-GNN) architectures build directly on this “graph diffusion + temporal evolution” paradigm: Li et al.’s DCRNN (Diffusion Convolutional Recurrent Neural Network) captures spatial dependencies via bidirectional diffusion on directed graphs paired with recurrent neural networks (RNNs) for temporal dynamics[5], while Yu et al.’s STGCN (Spatio-Temporal Graph Convolutional Network) interleaves graph and temporal convolutional blocks to model spatial and temporal dependencies jointly[6].

3.4 Why Is Transferring Models from Traffic Forecasting to Retail Non-Trivial? #

Road links in a traffic network have clear physical meaning, and vehicles physically travel along those roads. Edges in a retail network are closer to modelling hypotheses: geographic proximity does not guarantee cross-store shopping, sales correlation does not prove direct store-to-store influence, and warehouse supply links cannot share the same propagation rule as inter-store competition. Before transferring a traffic forecasting architecture to retail, a team must answer at least four questions:

• Are the edges credible? Do connections come from verifiable supply records and explicit distance rules, or merely from historical sales correlations?

• Does graph construction cause future data leakage? (In engineering, data leakage—or look-ahead bias—occurs when a model accidentally peeks at future information from the test period while training or predicting for today.) If edges are built from sales correlations computed across the entire dataset, has test-period information already leaked into the training graph?

• Does the graph change over time? When stores open or close, warehouse assignments shift, or delivery lead times change, do relationships learned during training still hold at inference time?

• Does the added complexity deliver incremental value? Can a conventional model augmented with neighbour sales, promotions, and warehouse states match the GNN’s accuracy?

Architectural diagrams alone cannot answer these questions. The next chapter constructs a fully controlled, reproducible synthetic retail network: we explicitly specify how warehouse states evolve and how demand propagates across relationships, then use an out-of-time test split to see where single-store models, feature-engineered models, and graph models diverge.

4. Synthetic Experiment: Does Graph Structure Provide Predictive Signal Beyond Traditional Reference Models? #

In our earlier article, Retail Demand Forecasting with Prophet: From Model Mechanics to Engineering Practice, we built a single-store synthetic dataset (2023-01-01 to 2025-12-31, 1,096 days in total, with the first 1,068 days for training and the final 28 days for testing). Under that single-store “local self-sufficiency” setting, Prophet achieved a Weighted Absolute Percentage Error (WAPE, total absolute error divided by total actual demand—lower is more accurate) of just 0.46% against the noise-free open-store true expected demand (denoted \(\mu_t^{\mathrm{open}}\)). Section 4.5 gives the full metric formulas used in this article.

We now expand that single-store benchmark into a retail network comprising 2 regional warehouses + 10 retail stores (across 3 commercial clusters) + asynchronous staggered promotion calendars. The data-generating process (DGP, the ground-truth mathematical rules we specify to simulate the dataset) combines store-opening and inventory hard constraints, warehouse Markov chains, and spatial-graph customer allocation mechanisms; forecasting models then attempt to recover these signals from observable features. The synthetic data generator, model comparisons, plotting scripts, and subsequent updates are available at GitHub: retail-gnn-demand-forecast.

A candid note on experimental design: a home-turf benchmark tailored to test GNNs. To put our cards on the table, this 2-warehouse, 10-store network is not a random sample from a live retail chain; it is a controlled “home-turf” experiment specifically designed to test and demonstrate GNN modelling capabilities—or, as a well-known Chinese idiom puts it, “wrapping a whole batch of dumplings (a synthetic dataset with cross-warehouse shocks and neighbour spillover) just for this dish of vinegar (evaluating GNN demand forecasting).” Because we deliberately baked warehouse disruptions, a 1-day logistics lead time, and two-hop customer spillover into the data-generating mechanism, the GNN will naturally come out ahead in the model comparison that follows. The real purpose of this controlled experiment is not to declare that “GNNs win across every retail scenario,” but to open the black box and show where single-store models and manual feature engineering hit their information ceiling when cross-entity propagation genuinely exists, and how graph structure recovers those missing signals. As for whether live retail systems satisfy these conditions—and when you should not deploy a GNN—Chapter 5 will examine those boundaries in detail.

4.1 Heterogeneous Retail Network Topology #

Following the notation in Sections 2.2 and 3.3, we define the retail network as a heterogeneous directed weighted graph \(G = (\mathcal V^{\mathrm W} \cup \mathcal V^{\mathrm S},\; \mathcal E^{\mathrm{sup}} \cup \mathcal E^{\mathrm{spa}})\), with warehouse nodes \(\mathcal V^{\mathrm W} = \{W_1, W_2\}\) and store nodes \(\mathcal V^{\mathrm S} = \{S_1, \dots, S_{10}\}\). Figure 1 illustrates the spatial layout of the network alongside its transition matrix:

As shown in Figure 1(a), the 10 stores are grouped into three commercial clusters. Most stores are supplied by their local regional warehouse, but inside the Eastern High-Density District (Cluster A) and the Central Corridor (Cluster B), we deliberately place two cross-warehouse bridge stores (\(S_4\) and \(S_6\)): geographically they sit right next to stores supplied by Eastern Hub \(W_1\), yet they receive their stock from Western Hub \(W_2\).

Why design the network this way? In real supermarket chains, convenience store networks, and dark-store (or front-warehouse) grids (small neighbourhood fulfilment centers closed to walk-in shoppers and dedicated exclusively to rapid online delivery), consumer shopping catchments and backend logistics service boundaries rarely line up neatly. Two stores across the street from each other may be assigned to different regional distribution centers (DCs, large central warehouses that pick and ship replenishment stock to stores across a region) because of municipal district lines, different acquisition or opening waves, trunk-route splitting, or capacity limits at a single DC; many retailers also deliberately assign “buffer stores” along warehouse boundary zones to a secondary DC to diversify stockout risk. This creates a common operational pattern: stores that compete for the exact same local foot traffic in the front end, but rely on different fulfilment hubs in the back end. Each store group’s coordinates, primary warehouse \(w(i)\), and spatial neighbourhood (\(d_{ij} \le 3.5\text{ km}\)) are summarized below:

Cluster and stores	Warehouse	Coordinates (km)	Spatial neighbourhood and topological role
A East High-Density: \(S_1, S_2, S_3\)	\(W_1\) (East Hub)	\((2.0, 8.0),\,(3.0, 8.4),\,(2.4, 7.0)\)	Pairwise adjacent (incl. \(S_4\)); core eastern stores constrained by \(W_1\) replenishment
A East High-Density: \(S_4\)	\(W_2\) (West Hub)	\((3.5, 7.3)\)	Adjacent to \(S_1, S_2, S_3\); eastern bridge store absorbing spillover when \(W_1\) stocks out
B Central Corridor: \(S_5, S_7\)	\(W_1\) (East Hub)	\((7.5, 5.0),\,(8.8, 4.6)\)	Flanking and adjacent to \(S_6\); constrained by \(W_1\) replenishment
B Central Corridor: \(S_6\)	\(W_2\) (West Hub)	\((8.3, 5.5)\)	Connects \(S_5, S_7\); central bridge store with two-way cannibalization and spillover
C West Center: \(S_8, S_9, S_{10}\)	\(W_2\) (West Hub)	\((14.0, 2.5),\,(15.0, 2.9),\,(14.4, 1.7)\)	Pairwise adjacent; western regional stores independently supplied by \(W_2\)

Echoing the division of labour between the two edge types in Section 3.3, the heterogeneous edges connecting these nodes map directly onto two physical chains in retail operations:

1. Warehouse–store supply edges \(\mathcal E^{\mathrm{sup}}\) (dashed lines in Figure 1a, representing the daily replenishment chain—the high-frequency retail cadence of ordering daily, picking overnight, and stocking shelves by morning): transmit the operating state of primary warehouse \(w(i)\) on day \(t-1\) to store \(i\), setting the store’s inventory cap \(V_{i,t}\) on day \(t\). How much stock store \(i\) can sell during day \(t\) depends on whether its upstream DC experienced congestion or stockout disruptions the night before (day \(t-1\)).

2. Store–store spatial edges \(\mathcal E^{\mathrm{spa}}\) (solid lines in Figure 1a, representing local catchment substitution): shoppers buying everyday groceries rarely travel across town when an item is on promotion or out of stock, but readily substitute within a local travel radius (set here to \(d_{ij} \le 3.5\text{ km}\), roughly a 10-minute drive or bike ride); within that radius, closer stores exert stronger promotional pull and stockout substitution. We therefore compute the raw spatial affinity matrix \(\mathbf A^{\mathrm{spa}}\) using a Gaussian distance-decay kernel (a smooth bell-shaped decay function that assigns values near \(1\) at very short distances and decays smoothly with distance, parameterized here by scale \(\sigma_d = 2.0\text{ km}\)) truncated at \(3.5\text{ km}\), and then—exactly as in Section 3.3—left-multiply by the inverse row-degree matrix \((\mathbf D^{\mathrm{spa}})^{-1}\) to turn raw proximity scores into the row-normalized spatial random-walk transition matrix \(\mathbf P^{\mathrm{spa}}\) whose rows sum to \(100\%\):

\[ \begin{aligned} A_{ij}^{\mathrm{spa}} &= \exp\!\left(-\frac{d_{ij}^2}{2\sigma_d^2}\right)\cdot \mathbb{I}(0 < d_{ij} \le 3.5\text{ km}), \\ \mathbf P^{\mathrm{spa}} &= (\mathbf D^{\mathrm{spa}})^{-1}\mathbf A^{\mathrm{spa}}. \end{aligned} \]

In the heatmap of Figure 1(b), \(\mathbf P^{\mathrm{spa}}\) exhibits a clean block-diagonal structure (all non-zero weights sit inside three square blocks along the main diagonal corresponding to the three districts, while off-diagonal entries across districts are \(0\)): store pairs more than \(3.5\text{ km}\) apart have no cross-district customer flow, whereas inside Cluster A (upper-left blue box) and Cluster B (middle orange box), stores supplied by Eastern Hub \(W_1\) and Western Hub \(W_2\) are tightly bound by high-weight spatial edges. Whenever \(W_1\) suffers a disruption, unmet demand at its affected stores naturally spills over along \(\mathbf P^{\mathrm{spa}}\) to bridge stores \(S_4\) and \(S_6\), which remain well stocked by \(W_2\).

4.2 Four-Layer Coupled Data-Generating Process (DGP) #

To stay strictly aligned with the single-store Prophet article, the simulation spans 2023-01-01 to 2025-12-31 (\(T=1,096\) days), using the first 1,068 days for training and the final 28 days (2025-12-04 to 2025-12-31) for testing. Daily observations are generated sequentially across four layers:

Layer 1: Discrete-time Markov chain for warehouse states (corresponding to Chapter 1)

Each warehouse \(w \in \{W_1, W_2\}\) has a daily operating state \(Z_{w,t}^{\mathrm{wh}} \in \{0\text{ (Normal)}, 1\text{ (Congested)}, 2\text{ (Disrupted)}\}\). The superscript \(\mathrm{wh}\) denotes warehouse state to avoid confusion with store nodes \(S_i\). It evolves according to the transition matrix:

\[ \mathbf P^{\mathrm{wh}} = \begin{bmatrix} 0.78 & 0.15 & 0.07 \\ 0.42 & 0.40 & 0.18 \\ 0.38 & 0.32 & 0.30 \end{bmatrix} \]

Because overnight store replenishment carries a 1-day logistics lead time, the effective inventory supply cap \(V_{i,t}\) at store \(i\) on day \(t\) is determined by its primary warehouse’s lagged state on day \(t-1\), \(Z_{w(i),t-1}^{\mathrm{wh}}\) (so when forecasting day-\(t\) demand, this state is an already observed historical input with zero future leakage):

\[ V_{i,t} = \begin{cases} 999\text{ (ample supply, uncensored)}, & \text{if } Z_{w(i),t-1}^{\mathrm{wh}} = 0 \\ 95\text{ (mild replenishment cap)}, & \text{if } Z_{w(i),t-1}^{\mathrm{wh}} = 1 \\ 50\text{ (severe disruption cap)}, & \text{if } Z_{w(i),t-1}^{\mathrm{wh}} = 2 \end{cases} \]

Layer 2: Asynchronous staggered promotions and single-store baseline demand (aligned with the Prophet article formula)

In real retail operations, stores frequently run staggered promotional schedules due to rotating flyer waves, store-format differences, or localized flash sales—creating day-to-day price differentials across neighbouring stores in the same catchment. We therefore assign each store an asynchronous promotion calendar with a phase offset \(\phi_i \in \{0, 1, \dots, 12\}\), combining regular promotions on a 14-day cycle (20% off, \(q_{i,t}=0.20\)) and deep flash promotions on a 21-day cycle (35% off, \(q_{i,t}=0.35\)). In the absence of cross-store interference, store \(i\)’s baseline self-demand \(\mu_{i,t}^{\mathrm{self}}\) follows the exact additive structure from the single-store article:

\[ \mu_{i,t}^{\mathrm{self}} = \beta_{0,i} + \beta_{1,i} \cdot t + \beta_{W,i} \cdot W_t + A_i \cdot \sin\left(\frac{2\pi(t-15)}{365.25}\right) + \beta_{H,i} \cdot H_t + \gamma_i \cdot q_{i,t} \]

The six terms on the right-hand side represent store \(i\)’s base daily foot traffic \(\beta_{0,i}\), long-term linear growth trend \(\beta_{1,i}t\), weekend shopping lift \(\beta_{W,i}W_t\), annual summer–winter seasonality \(A_i\sin(\cdot)\), public-holiday demand surge \(\beta_{H,i}H_t\), and direct sales lift from its own discount \(\gamma_i q_{i,t}\).

Layer 3: Sign-flipping bidirectional nonlinear coupling on the graph (corresponding to Sections 3.2 and 3.3)

Once each store’s baseline demand is set, stores within the same commercial cluster interact in both directions across the spatial transition matrix \(\mathbf P^{\mathrm{spa}}\) constructed in Section 4.1—notice how the matrix indices flip between \(P_{ij}^{\mathrm{spa}}\) and \(P_{ji}^{\mathrm{spa}}\) in the two equations below:

• One-hop neighbour promotion cannibalization \(C_{i,t}\) (reading row entries \(P_{ij}^{\mathrm{spa}}\), focal store’s shoppers drawn outward): When a neighbouring store \(j \in \mathcal{N}(i)\) offers a deeper discount than focal store \(i\) (\(q_{j,t} > q_{i,t}\)) and has ample shelf stock that day (\(V_{j,t} > 120\), i.e., its upstream warehouse is in state \(Z_{w(j),t-1}^{\mathrm{wh}} = 0\)), price-comparing shoppers who would normally visit focal store \(i\) are drawn away to neighbour \(j\) according to the outgoing weight \(P_{ij}^{\mathrm{spa}}\):

\[ C_{i,t} = 125.0 \sum_{j \in \mathcal{N}(i)} P_{ij}^{\mathrm{spa}} \cdot \max(0, q_{j,t} - q_{i,t}) \cdot \mathbb{I}(V_{j,t} > 120) \]

• Two-hop warehouse state transmission and neighbour stockout spillover \(L_{i,t}\) (reading entries \(P_{ji}^{\mathrm{spa}}\), stocked-out neighbour’s shoppers flowing inward): Conversely, when two-hop upstream warehouse congestion or disruption caps neighbour \(j\)’s inventory at \(V_{j,t}\) so that its net demand exceeds available stock, neighbour \(j\) sells out and leaves \(U_{j,t}\) units of unmet demand. Assuming \(15\%\) of those disappointed shoppers abandon the purchase while the remaining \(\eta = 85\%\) look for an alternative store in the same catchment, that unmet demand flows outward from stocked-out store \(j\) into focal store \(i\) according to share \(P_{ji}^{\mathrm{spa}}\) (using \(L\) for spillover inflow into store \(i\) to avoid confusion with store IDs \(S_i\)):

\[ \begin{aligned} U_{j,t} &= \max\!\big(0,\; \mu_{j,t}^{\mathrm{self}} - C_{j,t} - V_{j,t}\big), \\ L_{i,t} &= \eta \sum_{j \in \mathcal{N}(i)} P_{ji}^{\mathrm{spa}} \cdot U_{j,t},\qquad \eta=0.85. \end{aligned} \]

Layer 4: Sales censoring from store-opening and inventory observation rules

Combining these three layers yields the open-store conditional mean \(\mu_{i,t}^{\mathrm{open}}\) and noisy latent true demand \(D_{i,t}^{\mathrm{sim}}\) (where \(\epsilon_{i,t} \sim \mathcal{N}(0, 4.5^2)\)). Finally, subject to the store-opening guard variable \(O_{i,t} \in \{0, 1\}\) (all stores close on Good Friday and Christmas Day, setting \(O_{i,t}=0\)) and the inventory cap \(V_{i,t}\), the transaction system records observed sales \(y_{i,t}^{\mathrm{sales}}\) and stockout indicators \(\operatorname{stockout}_{i,t}\):

\[ \begin{aligned} \mu_{i,t}^{\mathrm{open}} &= \max\!\big(10,\; \mu_{i,t}^{\mathrm{self}} - C_{i,t} + L_{i,t}\big), \\ D_{i,t}^{\mathrm{sim}} &= \max\!\big(0,\; \mu_{i,t}^{\mathrm{open}} + \epsilon_{i,t}\big), \\ y_{i,t}^{\mathrm{sales}} &= O_{i,t} \cdot \operatorname{round}\!\big(\min(D_{i,t}^{\mathrm{sim}}, V_{i,t})\big), \\ \operatorname{stockout}_{i,t} &= \mathbb{I}\!\big(D_{i,t}^{\mathrm{sim}} > V_{i,t} \land O_{i,t} = 1\big). \end{aligned} \]

4.3 What Does Store-Level Sales Data Look Like? Network-Wide Trends and Single-Store Decomposition #

With the data-generating mechanism defined, let us inspect the resulting sales series across stores.

Unlike the isolated, smooth single-store series in Retail Demand Forecasting with Prophet: From Model Mechanics to Engineering Practice, daily store sales inside a multi-warehouse, multi-store network reflect both regular local cycles and structural shocks from the surrounding network.

Figure 2 plots daily observed sales and latent demand across all 10 stores during the second half of 2025 (July to December, 184 days):

Three patterns stand out in Figure 2:

1. Tiered store scale and asynchronous promotion rhythms: Average daily sales across the 10 stores range between 85 and 103 units. The yellow-shaded 14-day and 21-day promotion windows are staggered across stores, preventing the entire network from running identical promotions on the same day.

2. Widespread stockout censoring (red triangles): Whenever an upstream warehouse enters a congested state (cap of 95 units) or disrupted state (cap of 50 units), stores supplied by that warehouse see their observed sales \(y_{i,t}^{\mathrm{sales}}\) clipped at a horizontal ceiling—even during peak promotions—opening a wide gap between the grey dashed line (latent demand \(D_{i,t}^{\mathrm{sim}}\)) and realized sales.

3. Distinct behaviour at cross-warehouse bridge stores (S4 and S6, orange lines): Stores S4 and S6 sit next to stores supplied by \(W_1\), but receive their own stock from Western Hub \(W_2\). When \(W_1\) suffers a disruption, S4 and S6 are not bound by \(W_1\)’s inventory cap and absorb unmet demand from their neighbours in the simulation, creating sharp demand spikes.

To unpack how individual store series are formed, Figure 3 compares the demand components of core eastern store S1 and central bridge store S6:

Figure 3 separates three distinct effects acting on a single store:

• Neighbour promotion cannibalization dips (blue/yellow shading): When the focal store is not on promotion while an amply stocked neighbour runs a deep discount, the focal store’s sales fall below its own baseline demand \(\mu^{\mathrm{self}}\).

• Cross-store stockout spillover peaks (green shading): When a promoted neighbour hits an upstream warehouse stockout, unmet demand flows along spatial edges to the well-stocked focal store, pushing its open-store conditional mean \(\mu^{\mathrm{open}}\) well above its own baseline demand.

• Inventory-cap censoring loss (pink shading): When the focal store itself faces an upstream warehouse disruption, latent demand cannot fully convert into realized sales, clipping observed outcomes at the inventory cap.

A single store’s sales series therefore mixes its own calendar cycles, troughs caused by neighbour promotions, and spikes caused by neighbour stockouts. Looking only at the store’s own history makes it hard to disentangle these sources—which is why the next section compares neighbourhood features and graph models.

4.4 Leakage-Safe Model Comparison: Single-Store, Neighbourhood Features, and a Two-Layer Analytical Graph Model #

Following the modelling rules in Sections 2.6 and 3.4, every evaluated model is restricted during training and inference to observed sales \(y_{i,t}^{\mathrm{sales}}\), stockout indicators \(\operatorname{stockout}_{i,t}\), scheduled promotions \(q_{i,t}\), and one-day-lagged warehouse states \(Z_{w,t-1}^{\mathrm{wh}}\), with zero access to unobserved true latent demand \(D_{i,t}^{\mathrm{sim}}\).

Why use closed-form NumPy Ordinary Least Squares (OLS—computing the exact error-minimizing parameter vector in a single linear-algebra step without iterative gradient descent) to compare the four model families here before introducing deep-learning toolkits like PyG and LibCity in Section 4.6? In machine learning benchmarks, when an end-to-end trained deep GNN beats a traditional statistical baseline, a natural question arises: “Did error drop because the graph structure provided genuine information, or because neural-network hyperparameter tuning, nonlinear layers, or lucky random seeds gave it an edge?” To run a clean information-boundary ablation study (a controlled experiment that holds the solver fixed and adds one layer of relational information at a time to isolate each layer’s true contribution), we deliberately strip away gradient-descent variance and hyperparameter tuning so that all four models use the exact same closed-form least-squares solver—changing only the information horizon visible to the model (from single-store history \(\to\) one-hop flat neighbourhood features \(\to\) two-hop heterogeneous graph propagation).

Specifically, Models 1A, 1B, and 2 fit parameters independently for each store:

\[ \hat y_{i,t}^{(m)}=\mathbf{x}_{i,t}^{(m)\mathsf T}\hat{\boldsymbol\beta}_{i}^{(m)},\qquad \hat{\boldsymbol\beta}_{i}^{(m)}=\arg\min_{\boldsymbol\beta} \sum_{t\in\mathcal T_i^{(m)}} \left(y_{i,t}^{\mathrm{sales}}-\mathbf{x}_{i,t}^{(m)\mathsf T}\boldsymbol\beta\right)^2, \]

where \(m\in\{1A,1B,2\}\) and \(\mathcal T_i^{(m)}\) denotes the training dates available to model \(m\) for store \(i\). Model 3, mirroring the GNN parameter sharing in Section 2.3, pools uncensored training observations across all stores to estimate a shared relational coefficient vector:

\[ \hat y_{i,t}^{(3)}=\mathbf{x}_{i,t}^{(3)\mathsf T}\hat{\boldsymbol\beta}^{(3)},\qquad \hat{\boldsymbol\beta}^{(3)}=\arg\min_{\boldsymbol\beta} \sum_{(i,t)\in\mathcal T^{(3)}} \left(y_{i,t}^{\mathrm{sales}}-\mathbf{x}_{i,t}^{(3)\mathsf T}\boldsymbol\beta\right)^2. \]

Expanding the information horizon step by step across the four models creates a four-rung ladder for testing the value of graph structure:

1. Model 1A: Naive Single-Store Analytical Baseline (Naive Local OLS, 0-hop + censored data): Uses only store \(i\)’s own calendar and promotion features without filtering out historical stockout days, fitting directly on observed sales \(y_{i,t}^{\mathrm{sales}}\) across all open store-days.

2. Model 1B: Uncensored Single-Store Analytical Baseline (Uncensored Local OLS, 0-hop + de-censored): Uses the exact same input features as 1A, but filters training data using observed \(\operatorname{stockout}_{i,t} = 0\) to fit only on uncensored open days—representing a clean single-store baseline that makes the best possible use of the store’s own history.

3. Model 2: First-Order Neighbourhood Feature Engineering (1-Hop Flat Features, corresponding to Section 2.6): Extends Model 1B by flattening one-hop neighbourhood signals into a wide table—concatenating spatially weighted neighbour promotion intensity \(\sum_j P_{ij}^{\mathrm{spa}} q_{j,t}\), lagged neighbour sales \(\sum_j P_{ij}^{\mathrm{spa}} y_{j,t-1}^{\mathrm{sales}}\), and the focal store’s own lagged warehouse state \(Z_{w(i),t-1}^{\mathrm{wh}}\). Crucially, however, it cannot see two-hop “neighbour upstream warehouse states” \(Z_{w(j),t-1}^{\mathrm{wh}}\) and does not condition neighbour promotions on neighbour supply availability.

4. Model 3: Two-Layer Analytical Graph Model (corresponding to the two-hop path in Sections 2.4 and 3.2): Faithfully reproduces the two-hop “neighbour warehouse \(\to\) neighbour store \(\to\) focal store” message-passing path described in Section 3.2 (using strictly observable features and never touching \(D_{i,t}^{\mathrm{sim}}\) or \(\mu_{i,t}^{\mathrm{open}}\)):

   • Hop 1 (Warehouse \(\to\) Store, along supply edges \(\mathcal E^{\mathrm{sup}}\)): At each neighbouring store \(j\), combines its lagged upstream warehouse state \(Z_{w(j),t-1}^{\mathrm{wh}}\) with its same-day promotion \(q_{j,t}\) and calendar features to form a supply-conditioned first-layer representation \(\mathbf{h}_{j,t}^{(1)}\) that distinguishes “well-stocked promotion” from “stockout-capped promotion.”

   • Hop 2 (Store \(\to\) Store, along spatial edges \(\mathcal E^{\mathrm{spa}}\)): Aggregates those first-layer neighbour representations \(\mathbf{h}_{j,t}^{(1)}\) into focal store \(i\) via the spatial transition matrix \(\mathbf P^{\mathrm{spa}}\), allowing the focal store to separate one-hop cannibalization from two-hop cross-warehouse spillover.

Figure 4 illustrates the dynamic coupling between two adjacent stores in the Central Corridor across the 28-day test window:

• In Figure 4(a), store \(S_5\) (supplied by Eastern Hub \(W_1\)) experiences an upstream disruption from 12-08 to 12-12; although its latent demand reaches 170–180 units, observed sales \(y_{5,t}^{\mathrm{sales}}\) (grey dotted line) are clipped at \(V_{5,t}=50\) units, leaving the pink-shaded unmet demand \(U_{5,t}\).

• In Figure 4(b), cross-warehouse bridge store \(S_6\)—located right next to \(S_5\) and \(S_7\) but supplied by healthy warehouse \(W_2\)—absorbs this spillover traffic, pushing its latent demand to 201 units on 12-08. The single-store model (orange dashed line) and first-order neighbourhood feature model (purple dash-dotted line) underpredict this peak by more than 70 units, whereas the analytical graph model (green solid line) recovers most of the spillover peak and subsequent trough. This gap reflects the two-hop cross-warehouse transmission mechanism explicitly built into the synthetic generator.

4.5 Quantitative Evaluation and Ablation Across Three Business Regimes #

Across the held-out 28-day test set (270 open store-days), we evaluate predictions \(\hat y_n\) against targets \(y_n\) using three metrics:

\[ \operatorname{WAPE}=\frac{\sum_n|\hat y_n-y_n|}{\sum_n|y_n|},\qquad \operatorname{MAE}=\frac{1}{N}\sum_n|\hat y_n-y_n|,\qquad \operatorname{Bias}=\frac{\sum_n(\hat y_n-y_n)}{\sum_n y_n}. \]

Here, \(N\) is the number of observations in the evaluation subset; because all demand and sales targets in this study are non-negative with positive subset totals, the denominators of WAPE and Bias are well defined. WAPE (Weighted Absolute Percentage Error) measures total absolute error relative to total target volume, MAE (Mean Absolute Error) retains the intuitive unit of “average error in units per store-day,” and signed Bias (relative systematic bias) indicates overall overprediction (positive) or underprediction (negative). Comparing the same predictions against the noise-free conditional mean, noisy latent demand, and observed sales yields the table below (each cell reports WAPE / MAE / Bias):

Evaluated model	Noise-free mean \(\mu_{i,t}^{\mathrm{open}}\)	Latent demand \(D_{i,t}^{\mathrm{sim}}\)	Observed sales \(y_{i,t}^{\mathrm{sales}}\)
1A. Naive Local OLS	17.62% / 19.97 / -13.97%	18.14% / 20.54 / -13.87%	21.14% / 20.45 / +0.78%
1B. Uncensored Local OLS	12.99% / 14.72 / -5.74%	13.78% / 15.59 / -5.63%	20.02% / 19.37 / +10.43%
2. 1-Hop Flat Features	11.79% / 13.36 / -5.17%	12.75% / 14.43 / -5.06%	14.40% / 13.93 / +0.17%
3. Two-Layer Analytical Graph	1.92% / 2.17 / -1.01%	3.79% / 4.29 / -0.90%	3.18% / 3.08 / -0.21%

The table highlights the mechanism behind each model family: Model 1A suffers from historical stockout censoring and systematically underpredicts latent demand (Bias -13.87%); Model 1B reduces own-stockout censoring bias by filtering out stockout days, yet remains blind to neighbour cannibalization and cross-store spillover; Model 2 adds first-order linear neighbourhood features but lowers latent-demand WAPE by only about 1 percentage point over 1B; and Model 3, which matches the generator’s two-hop interaction structure, reduces latent-demand WAPE to 3.79%.

To directly answer the question raised in Section 2.6—“When is a single-store model already sufficient, and when is a graph-structured model worth evaluating?”—we slice the test set into three business regimes (because regime conditions can overlap, row sample sizes do not sum to the full test set; metrics are evaluated against latent true demand \(D_{i,t}^{\mathrm{sim}}\)):

Business regime (Test Regimes)	Sample size	1A. Naive Local WAPE (Bias)	1B. Uncensored Local WAPE (Bias)	2. 1-Hop Flat WAPE (Bias)	3. Analytical Graph WAPE (Bias)
Normal Regime	57 store-days	13.50% (-13.50%)	3.73% (+0.05%)	4.06% (+1.72%)	2.74% (+0.11%)
Cannibalization Regime	142 store-days	10.13% (+0.35%)	11.40% (+5.72%)	10.15% (+4.79%)	4.24% (+0.62%)
Spillover Regime	88 store-days	29.23% (-29.23%)	22.35% (-22.24%)	21.11% (-20.65%)	3.95% (-3.17%)

Combining the regime table with the systematic bias in Figure 5(b) yields a clear engineering takeaway: in the Normal Regime without cross-store shocks, the uncensored single-store model (1B) already reaches 3.73% WAPE, leaving limited incremental room for a graph model. During shock regimes, however, the effect of a neighbour promotion \(q_{j,t}\) flips sign depending on upstream supply—when the neighbour’s warehouse is normal, it causes negative cannibalization and produces a +5.72% overprediction bias if unmodelled; when the neighbour’s warehouse is disrupted, it causes positive spillover and produces a -22.24% underprediction bias if unmodelled. Including neighbour promotion only as a first-order linear feature (Model 2) causes these opposite-sign effects to partially cancel across training samples. Our two-layer construction first forms the “warehouse state \(\times\) neighbour promotion” interaction and then aggregates across the spatial graph, allowing it to separate two mechanisms that operate under different conditions with opposite signs. This advantage follows from the prescribed data-generating process and does not imply that only graph models can represent such interactions; in real-world projects, gradient-boosted decision trees (such as LightGBM and XGBoost) equipped with explicit interaction features should also be included as strong baselines. To keep our experimental narrative focused on the progression from time-series models to graph propagation, we omit separate tree-model benchmarks here—readers interested in tree-based tabular forecasting and single-store Prophet modelling can refer to our earlier articles, From a Single Tree to a Forest: Mathematical Principles, Algorithmic Evolution, and Engineering Practices of Tree Models and Retail Demand Forecasting with Prophet: From Model Mechanics to Engineering Practice.

Up to this point, however, all four comparisons have relied on hand-written NumPy closed-form Ordinary Least Squares (OLS). In day-to-day engineering, practitioners rarely write custom closed-form solvers; they reach for mature open-source or production forecasting libraries. That brings us to the second stage of our controlled experiment: if we hold the three de-censored information horizons fixed (Model 1B single-store, Model 2 one-hop neighbourhood, and Model 3 two-hop spatio-temporal graph) and replace the underlying solver with mainstream time-series and graph-learning toolkits, can single-store or one-hop models close the accuracy gap through library-level modelling power? And what accuracy and runtime cost do genuine graph neural network frameworks deliver?

4.6 Panorama Benchmark Across Mainstream Time-Series and Graph-Learning Toolkits #

Whereas Sections 4.4–4.5 held the solver fixed (closed-form OLS) and varied the information horizon, this section holds each information horizon fixed and swaps the underlying software toolkit and solver mechanism. Using the exact same leakage-safe data split in the companion repository, we feed the three information tiers corresponding to Model 1B, Model 2, and Model 3 into four mainstream open-source toolkits alongside the closed-form analytical baselines.

Before examining the scorecard, let us briefly introduce the four representative toolkits evaluated here and their core mechanisms:

• Meta Prophet: Meta’s open-source decomposable statistical time-series library, which models a series as the sum of trend, multi-period Fourier seasonality (superposing sine and cosine waves to fit weekly and annual cycles), holiday effects, and extra regressors (add_regressor), using the probabilistic programming language Stan’s quasi-Newton optimizer (L-BFGS) for fast Maximum A Posteriori (MAP) parameter estimation—well suited to business series with strong calendar structure and interpretability requirements.
• Darts: An open-source Python library by Unit8 for unified time-series forecasting, wrapping classical statistical models, linear regression, and deep forecasting architectures behind a standardized TimeSeries abstraction for multi-series training and backtesting.
• PyTorch Geometric (PyG) / PyG Temporal: A leading PyTorch-based graph neural network library with native support for heterogeneous graphs (HeteroData) and message-passing operators (MessagePassing), making it straightforward to implement custom two-hop spatio-temporal message passing along warehouse\(\to\)store supply edges and store\(\leftrightarrow\)store spatial edges.
• LibCity: An open-source unified benchmarking library for spatio-temporal data mining (widely used in traffic flow, urban grid demand, and multi-store retail forecasting), packaging spatio-temporal graph dataset slicing, standardized ST-GNN training pipelines, and evaluation metrics into an end-to-end framework.

Panorama Toolkit Benchmark Scorecard #

The test set again spans 2025-12-04 to 2025-12-31 across 270 open store-days. With each toolkit’s solver mechanism introduced above, the table below groups the results by the three de-censored information tiers from Section 4.4 (1B Single-store, 2 1-Hop, and 3 2-Hop Graph), reporting single-run wall-clock runtime and regime-specific accuracy on identical hardware (evaluated against true latent demand \(D_{i,t}^{\mathrm{sim}}\)):

Tier and toolkit	Time	Full / Normal WAPE	Cannib. WAPE (Bias)	Spillover WAPE (Bias)
1B Single · Formula OLS	0.005s	13.78% / 3.73%	11.40% (+5.7%)	22.35% (-22.2%)
1B Single · Meta Prophet	25.4s	13.98% / 3.84%	11.73% (+5.6%)	22.43% (-22.4%)
1B Single · Darts	12.1s	13.78% / 3.73%	11.40% (+5.7%)	22.35% (-22.2%)
2 1-Hop · 1-Hop Flat OLS	0.010s	12.75% / 4.06%	10.15% (+4.8%)	21.11% (-20.7%)
2 1-Hop · Prophet add_regressor	38.6s	12.89% / 4.09%	10.48% (+4.9%)	21.08% (-20.7%)
3 2-Hop Graph · Analytical ST-GNN	0.015s	3.79% / 2.74%	4.24% (+0.6%)	3.95% (-3.2%)
3 2-Hop Graph · PyG Temporal	4.8s	3.82% / 2.74%	4.24% (+0.6%)	3.98% (-3.2%)
3 2-Hop Graph · LibCity	18.5s	3.80% / 2.74%	4.24% (+0.6%)	3.96% (-3.2%)

Accuracy Tiers and Runtime Overhead: Why Swapping Libraries Cannot Fix an Information Blind Spot #

Reading the scorecard alongside Figure 6 highlights two engineering takeaways—one vertical (accuracy across information tiers) and one horizontal (compute cost across toolkits):

1. Across accuracy tiers: the information boundary sets the ceiling, and swapping libraries cannot fix an information blind spot. Within the same information tier, whether we run a 0.005s NumPy closed-form OLS solver, a 25.4s Meta Prophet model, or a 12.1s Darts pipeline, single-store full-test WAPE stays locked between 13.78% and 13.98%—and between 22.35% and 22.43% in the stockout spillover regime. Likewise, in the one-hop neighbourhood tier, 1-Hop Flat OLS and Prophet add_regressor produce virtually identical errors (roughly 12.8% full WAPE and 21.1% spillover WAPE). Only when the model is given access to the two-hop “warehouse\(\to\)neighbour store\(\to\)focal store” interaction structure do Analytical ST-GNN, PyG message passing, and LibCity all drop together to 3.79%–3.82% full WAPE (and roughly 4.0% spillover WAPE). This confirms our central thesis: if a model never observes the two-hop interaction between upstream warehouse states and neighbour promotions, no outer software wrapper can conjure that missing signal out of thin air.

2. Across runtime and ecosystem layers: compute overhead must match the required level of engineering abstraction. Figure 6(a) shows that single-run runtime spans nearly four orders of magnitude—from 0.005 seconds for Formula OLS to 38.6 seconds for Prophet with exogenous regressors—while Figure 6(b) organizes the tools into four abstraction layers: NumPy / SciPy analytical solutions provide a zero-overhead reference for verifying mechanisms and metrics; Prophet / Darts offer turnkey time-series engineering interfaces; PyG / DGL provide customizable heterogeneous graph message-passing operators; and LibCity offers standardized spatio-temporal benchmarking pipelines. Because each layer serves a different engineering purpose, they cannot be ranked by a single runtime figure or WAPE number alone. In production systems, adopting a graph framework entails not only seconds or minutes of extra training time, but also tighter requirements on feature latency, fault isolation, and MLOps maintenance—precisely the engineering boundaries we examine in Chapter 5.

5. The Gap Between Theory and Practice: Why a Graph Algorithm That Wins on Paper May Be Unexecutable in Production #

Seeing the experimental results in Chapter 4—where the two-hop heterogeneous graph model slashes spillover-regime forecast error from 22.35% down to 3.95%—the most natural impulse is to greenlight a production GNN project immediately. In industrial machine learning, however, few illusions are more dangerous than equating a paper victory in a controlled benchmark with an executable solution in production.

Chapter 4 was able to showcase the power of graph structure so cleanly because we operated from a “God’s-eye view,” taking for granted four premises that are routine in a laboratory but rarely hold in live retail operations: perfectly punctual data, noise-free relational topology, isolated node failures, and frictionless execution of forecasts. Once those four laboratory guardrails are removed, a graph neural network that looks unimpeachable in mathematical derivations and offline backtests may prove unsafe—or outright impossible—to execute in a live replenishment pipeline. Bridging the gap between theory and practice requires overcoming three real-world hurdles in sequence.

5.1 Hurdle 1: When Theoretical “Perfect States” Collide with Real-World Data Distortion and Latency Inversion #

In the equations of stochastic processes and graph neural networks, observed historical sales \(y_{i,t}^{\mathrm{sales}}\), lagged warehouse states \(Z_{w,t-1}^{\mathrm{wh}}\), and the spatial transition matrix \(\mathbf P^{\mathrm{spa}}\) are treated as crisp, deterministic inputs that arrive right on schedule. On the retail floor, however, engineers must first contend with a messy, distorted data foundation:

1. The Illusion of Aggregate Metrics and the Censoring Feedback Loop In theoretical modelling, we routinely rank models by their average MAE or WAPE across the full test set; in operations, an aggregate average is often a smokescreen that hides severe local failures. Recall Model 1A in Section 4.5 (the single-store baseline fitted directly on all observed sales): its overall sales Bias across the test set is just +0.78%, which looks almost perfectly unbiased on a dashboard; yet during the critical cross-store stockout spillover regime, its true Bias against latent demand plunges to -29.23%. Worse still, real retail systems rarely have shelf sensors that record the exact minute an item sells out. If a model trains directly on stockout-censored sales and feeds forecasts that underpredict demand by nearly 30% into an automated replenishment engine, it triggers a self-reinforcing deadlock: stockout \(\to\) censored sales \(\to\) depressed forecast \(\to\) smaller replenishment order \(\to\) even faster stockout. The model’s average error on executive reports may look reassuringly low while physical store shelves grow emptier by the week.

2. Phantom Inventory and Latency Inversion (Phantom Inventory & Data Latency) In the synthetic experiment of Chapter 4, the upstream warehouse state \(Z_{w,t-1}^{\mathrm{wh}}\) is observed without measurement error and is guaranteed to land in the feature store before the nightly forecasting job starts. In real retail operations, book inventory persistently drifts away from shelf-available stock: fresh spoilage, shrinkage, shoplifting, missed checkout scans, items sitting in shoppers’ carts, and misplaced backroom pallets all create phantom (or ghost) inventory—stock that exists in the database while the physical shelf is bare. At the same time, upstream sorting congestion and trunk-line transit delays are often logged only after store staff finish unloading and scanning pallets the next morning, or during end-of-month reconciliation. When the cross-node signals a graph model relies on arrive after the inference cutoff (latency inversion), even the most elegant two-hop graph operator is left multiplying missing or stale features at prediction time.

3. Dilution of Spatial Edges by Unobserved Confounders Our synthetic generator assumes that customer switching between stores follows a smooth Gaussian decay with geographic distance. Real foot-traffic shifts are simultaneously buffeted by dozens of confounders that never appear in enterprise tables: sudden roadworks blocking a store entrance, localized afternoon thunderstorms, a competitor convenience store across the street running a grand-opening giveaway, nearby office lease vacancies, or school holidays. These unobserved confounders severely dilute the signal-to-noise ratio of cross-store edges: a trained GNN may learn genuine, repeatable substitution patterns, or it may overfit a two-week coincidence caused by street construction into a spurious edge weight. Once the roadworks end and external conditions shift, that learned “pseudo-edge” will keep injecting systematic errors into production forecasts.

5.2 Hurdle 2: When Graph “Message Passing” Turns into a Production “Fault Amplifier” #

In textbooks and research papers, message passing is portrayed as an elegant cooperative mechanism: each node draws complementary signals from its neighbours to enrich its representation. From the perspective of site reliability and production engineering, however, every cross-node aggregation step intentionally breaches the system’s fault-isolation boundary.

1. Cascading Amplification from Single-Node Glitch to Subgraph Contamination In a conventional single-store architecture (such as per-store Prophet or LightGBM pipelines), each store’s forecast job acts as a natural firewall: even if store \(S_5\) reports wildly corrupted data because of a Point-of-Sale (POS, store checkout register) network outage, a faulty handheld barcode scanner during stocktaking, or a stalled Extract-Transform-Load (ETL, data ingestion and cleansing) pipeline (for instance, inventory suddenly flipping negative or sales spiking a hundredfold), the blast radius (a reliability-engineering term for the maximum scope of damage a single failure can cause) stays strictly confined to store \(S_5\), leaving replenishment across the other nine stores untouched. Switching to a two-hop graph model changes the failure mode fundamentally: a neighbour’s corrupted input is no longer a private local issue, but propagates within milliseconds through matrix multiplications with \(\mathbf P^{\mathrm{spa}}\) and \(\mathbf P^{\mathrm{sup}}\) to contaminate every connected store in the subgraph. If the node reporting corrupted data happens to be a cross-warehouse bridge node like \(S_4\) or \(S_6\), a localized sensor glitch can be amplified by the graph propagation operator into a district-wide forecasting avalanche. Without node-level anomaly circuit breakers (which automatically sever outgoing message propagation from a node whose input features fail sanity checks, much like an electrical fuse) and local fallback policies (automatically reverting a store to its standalone single-store forecast when neighbourhood signals are unreliable), a production GNN acts as a high-risk fault amplifier.

2. Store Heterogeneity, Combinatorial SKU Explosion, and Over-Smoothing Our synthetic benchmark contains just 10 homogeneous stores and 1 abstract product; a large retail chain manages thousands of stores—ranging from 50-square-meter corner convenience shops and 500-square-meter neighbourhood grocers to 10,000-square-meter warehouse clubs—alongside tens of thousands of active Stock Keeping Units (SKUs, individual product items distinguished to the finest level of package size, flavour, and variant) with vastly different customer demographics, shelf facings (the number of identical product packages displayed side by side on the front edge of a shelf facing the aisle—more facings mean higher visual exposure), and price elasticities. If a team attempts to pack store proximity, warehouse supply, and cross-SKU product substitution into a single massive heterogeneous graph:

   • On the infrastructure side, edge counts explode combinatorially, causing graph sampling, GPU memory footprints, and relational feature-alignment pipelines to blow past compute budgets;
   • On the modelling side, stacking deeper graph convolution layers easily triggers over-smoothing, averaging distinct neighbourhood shops and flagship hypermarkets into nearly indistinguishable node embeddings and erasing the local store characteristics that matter most.

3. The Operational Quagmire of Dynamic Topology and Cold Starts Retail networks never stand still: new store openings, remodel closures, pop-up holiday dark stores, and quarterly trunk-route reassignments continually rewrite the node set \(\mathcal V_t\) and edge set \(\mathcal E_t\). When a single-store model encounters a new store, it simply falls back to a standard single-store cold-start template (making initial forecasts for a newly opened store or newly launched item before sufficient historical data accumulates). When a graph model encounters a topology change, engineers must not only recalibrate neighbourhood normalization weights (otherwise a change in neighbour degree immediately shifts message scales), but also verify that aggregation weights learned on the old topology remain stable on the rewired graph.

5.3 Hurdle 3: Why “High Forecast Accuracy” Still Does Not Make Algorithmic Commands Directly Executable #

Suppose an engineering team overcomes all of the hurdles above—implementing stockout uncensoring, real-time feed validation, and node-level circuit breakers so that the graph model consistently delivers lower forecast error than baseline models in production. Can the model’s forecasts or transshipment recommendations now be wired directly into operational systems for automatic execution?

Still no—because probabilistic forecasting and deterministic business execution inhabit two fundamentally different mathematical worlds.

Statistical and machine learning models optimize for soft constraints and expected value: even when a graph model achieves a 3.8% WAPE, it is still answering the probabilistic question, “Under the historical distribution, what is the expected demand at store \(S_6\) tomorrow?”—and it is allowed to be wrong on individual instances. Core enterprise transaction and fulfilment systems—including Enterprise Resource Planning (ERP, managing procurement and financial ledgers), Warehouse Management Systems (WMS, managing bin locations, picking, and dispatch), and Order Management Systems (OMS, managing order routing and fulfilment status)—are governed by hard constraints and state-machine invariants (non-negotiable business rules that must hold before and after every state transition, such as physical sellable inventory never dropping below zero): even if an algorithm is right 99.9% of the time, a single command that violates a physical or financial invariant can trigger a severe operational incident.

Consider a concrete scenario: our graph forecasting module detects congestion at Eastern Hub \(W_1\) alongside tomorrow’s promotion at \(S_5\), and accurately predicts a +45% spillover demand surge at cross-warehouse neighbour \(S_6\). A downstream decision agent immediately issues what looks like a textbook command: “Tonight, transfer 80 cases from store \(S_4\) (also supplied by Western Hub \(W_2\)) to bolster store \(S_6\).” If the system lacks strict execution boundaries and lets the algorithm call inventory APIs directly, tomorrow morning will run headlong into a wall of operational constraints:

• Physical and fulfilment hard constraints: Of the 100 cases shown on store \(S_4\)’s book ledger, 35 cases are already locked and reserved by pending Online-to-Offline (O2O) orders (orders placed on mobile delivery apps to be picked directly from the local store and delivered within the hour), and another 10 cases sit shrink-wrapped at the back of a walk-in cooler awaiting breakdown—leaving fewer than 80 cases of genuinely transferable stock. Forcing the deduction directly causes an inventory breach (overselling / stockout default) on \(S_4\)’s committed online orders, triggering order cancellations and customer penalties.
• Logistics and scheduling boundaries: Nightly inter-store cold-chain trucks may require a minimum dispatch load of 100 cases, or store \(S_6\)’s night receiving dock may already be fully booked, leaving no crew to unload an unscheduled transfer.
• System safety and idempotency (a core distributed-systems property guaranteeing that if the exact same command is sent multiple times due to network retries, the system executes the stock deduction and truck dispatch once and only once): If the Warehouse Management System (WMS) fails to acknowledge the request within 3 seconds due to a network blip, will an automatic retry from the AI agent be recognized as the same transaction, or will it create a duplicate transfer order and put 160 cases on the road?

As large language models (LLMs) and AI agents rapidly enter enterprise supply-chain workflows, this distinction—forecasting is forecasting, execution is execution—is more critical than ever. A highly accurate demand forecast is never a license to bypass inventory reservation locks and oversell stock; a plausible chain of agent reasoning is never a license to bypass approval guards or idempotency checks and mutate financial ledgers directly. A production-grade intelligent replenishment architecture must strictly separate probabilistic algorithms from deterministic business rules across five distinct layers:

System layer	Core responsibility	Governing mathematical / engineering principle	Representative techniques and safeguards
1. Observation & Cleansing	Identify store closures, stockout censoring, and corrupted feeds	Restore true observation boundaries; block dirty data from entering the graph	Censoring flags, cycle-count calibration, node-level anomaly circuit breakers
2. Forecasting (Focus of this post)	Estimate future states or uncensored latent demand distributions	Probabilistic inference and conditional expectation (bounded error permitted)	Prophet, tree models, heterogeneous ST-GNN, local per-node fallbacks
3. Decision & Planning	Solve for replenishment / transfer plans under cost and capacity limits	Operations research and constrained optimization (balancing risk vs. return)	Mixed-integer programming (MIP), safety-stock policies, transfer heuristics
4. Business Execution	Execute cross-service state transitions reliably; enforce hard red lines	Deterministic finite-state machines (FSM) and hard invariants	On-hand reservation checks, authorization guards, ACID transactions, idempotency keys
5. Monitoring & Auditing	Track live regime-specific errors and system health	Closed-loop feedback and distribution-drift detection	Regime WAPE / Bias tracking, topology-drift alerts, immutable audit logs

Only when the Layer 2 graph model rests securely atop Layer 1’s data circuit breakers—and operates strictly inside Layer 4’s deterministic state-machine guardrails—can the offline accuracy gains demonstrated in Chapter 4 be converted safely into real operational value.

Conclusion: From Seeing Relationships to Respecting Boundaries #

Stepping back across the entire journey—from single-store time series through stochastic processes, graph propagation operators, and heterogeneous spatio-temporal graph networks, all the way to live retail engineering—every chapter has really been wrestling with the same foundational systems question: To face an uncertain tomorrow, what information about the past and the surrounding environment must a system retain, and in what form?

• When a store evolves self-sufficiently from its own history, a single-store time series is the most economical and resilient state record;
• When upstream fulfilment regimes switch randomly over time, Markov chains provide the probabilistic language to characterize state transitions;
• When a single node’s own history can no longer determine its future conditional distribution—because neighbour promotions and cross-warehouse disruptions are reshaping local demand across the network—a graph neural network uses message passing to encode the states of connected warehouses and neighbouring stores into the focal node’s hidden representation \(\mathbf h_{i,t}^{(k)}\);
• And when an algorithm finally reaches out to issue commands in the physical world, deterministic state machines and hard business invariants draw the non-negotiable safety boundaries around probabilistic predictions.

Looking ahead at the future of graph learning in industrial and retail supply-chain systems, the era of chasing leaderboard gains simply by stacking deeper graph convolution layers on fixed static graphs is fading. Three more durable directions are taking shape:

1. From Correlation Graphs to Causal and Dynamic Sparse Graphs (Causal & Dynamic Graph Learning) Next-generation industrial graph models are moving beyond fitting surface correlations on static distance matrices, increasingly integrating causal inference and sparse structure learning. In real environments full of unobserved confounders—weather shocks, roadworks, competitor promos—causal graph learning aims to prune away transient, spurious “pseudo-edges” and retain only cross-entity links backed by genuine causal propagation mechanisms, adapting dynamically as stores open, close, and rewire their supply lines.

2. From Monolithic End-to-End Graphs to Time-Series Foundation Models + Lightweight Graph Adapters As time-series foundation models and strong tabular tree pipelines mature, production forecasting architectures will rarely need to train massive end-to-end GNNs from scratch across every SKU and store. A far more pragmatic paradigm is emerging: use pre-trained time-series foundation models or gradient-boosted trees to anchor local store trends and seasonality, then attach a lightweight, low-hop, interpretable Graph Adapter on top dedicated specifically to correcting structural residuals caused by promotion cannibalization, stockout spillover, and cross-warehouse fulfilment shocks.

3. From Isolated Accuracy Metrics to End-to-End “Predict-Then-Optimize” and Guardrailed Agentic Workflows The value of a graph model will increasingly be judged not by an isolated offline WAPE number, but by its role as the relational sensing layer of a supply-chain digital twin (a virtual computational replica of the physical warehouse–store network that evolves in sync with real-world operations). Its graph-structured uncertainty estimates will feed directly into downstream operations-research solvers (End-to-End Predict-Then-Optimize) and, under the strict supervision of deterministic state machines and operational guardrails, empower AI agents to close the loop safely and idempotently from anomaly detection to collaborative cross-store fulfilment.

Seeing the true relationships among entities tells us when to introduce graph structure to break through the single-store information ceiling; respecting real-world data noise and execution boundaries tells us when to stop piling on model complexity and let disciplined systems engineering guard the bottom line. Holding fast to both yardsticks is what allows an engineer to neither be boxed in by isolated single-series thinking nor swept away by the hype of ever-larger models.

Reproduction note: The companion GitHub repository contains the Chapter 4 synthetic data generator, multi-toolkit benchmarks, runtime logs, and plotting scripts, where subsequent experiments and verification updates will continue to be maintained.

References #

[1] Gilmer, J., et al. Neural Message Passing for Quantum Chemistry. Proceedings of the 34th International Conference on Machine Learning, 2017. https://proceedings.mlr.press/v70/gilmer17a.html

[2] Kipf, T. N., & Welling, M. Semi-Supervised Classification with Graph Convolutional Networks. International Conference on Learning Representations, 2017. https://arxiv.org/abs/1609.02907

[3] Hamilton, W. L., Ying, R., & Leskovec, J. Inductive Representation Learning on Large Graphs. Advances in Neural Information Processing Systems, 2017. https://proceedings.neurips.cc/paper/2017/hash/5dd9db5e033da9c6fb5ba83c7a7ebea9-Abstract.html

[4] Veličković, P., et al. Graph Attention Networks. International Conference on Learning Representations, 2018. https://arxiv.org/abs/1710.10903

[5] Li, Y., Yu, R., Shahabi, C., & Liu, Y. Diffusion Convolutional Recurrent Neural Network: Data-Driven Traffic Forecasting. International Conference on Learning Representations, 2018. https://arxiv.org/abs/1707.01926

[6] Yu, B., Yin, H., & Zhu, Z. Spatio-Temporal Graph Convolutional Networks: A Deep Learning Framework for Traffic Forecasting. Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, 2018. https://doi.org/10.24963/ijcai.2018/505