A Critical Audit of Spatiotemporal Forecasting Benchmark Datasets and Baselines

# A Critical Audit of Spatiotemporal Forecasting Benchmark Datasets and Baselines
Source: [https://arxiv.org/html/2608.20980](https://arxiv.org/html/2608.20980)
Kenneth MartinThanks:Equal contribution\.Affiliation:Imperial College LondonAffiliation:London, United KingdomEmail:[k\.martin24@imperial\.ac\.uk](mailto:)Simon Heilig11footnotemark:1Asja FischerAffiliation:Ruhr University BochumAffiliation:Bochum, GermanyMichel F\. C\. HaddadAffiliation:Queen Mary University of LondonAffiliation:London, United KingdomAdam M\. SykulskiAffiliation:Imperial College LondonAffiliation:London, United KingdomMoshe EliasofAffiliation:Ben\-Gurion University of the NegevAffiliation:Beersheba, Israel

###### Abstract

Graph neural networks \(GNNs\) are routinely employed for short\-range forecasting on multivariate time series with a spatial graph structure\. Despite the availability of many alternative datasets, method innovations within this domain are predominantly assessed against a rather limited set of benchmark datasets, most notably Chickenpox, PedalMe, WikiMaths, METR\-LA, and PEMS\-BAY\. The evaluation protocols contain baselines spanning from historical averages to classical machine learning approaches\. These baselines often show competitive performance compared to GNNs\. In the present work, we take a step back and analyse the benchmark datasets via classical time series methods to uncover why spatially\-unaware linear models pose a stronger competitor than previously reported, casting further doubt on the discriminative reliability of the aforementioned widely adopted datasets\. Our statistical analysis provides a toolset for identifying significant spatial and temporal correlations, while revealing a structural bias introduced by first\-order differenced datasets\. We therefore recommend reducing the over\-reliance on such datasets for method comparison, and instead advocate for more rigorous statistical evaluation\. By applying the results of our analysis to a simple hybrid model, we show how our methodology can lead to novel ways of developing GNN models\.

## 1Introduction

Graph neural networks \(GNNs\) have seen great promise in applications requiring spatial information processing\[[37](https://arxiv.org/html/2608.20980#bib.bib35),[6](https://arxiv.org/html/2608.20980#bib.bib36),[9](https://arxiv.org/html/2608.20980#bib.bib37)\]due to their inductive bias stemming from the message passing mechanism, aggregating information over irregular spatial measurements\. With the advancement of filters for signals defined on a graph\[[7](https://arxiv.org/html/2608.20980#bib.bib34),[19](https://arxiv.org/html/2608.20980#bib.bib33),[26](https://arxiv.org/html/2608.20980#bib.bib32)\], generalizing the convolution neural network \(CNN\) beyond grid\-structured data, and more recently the development of neural ODE\-inspired GNNs\[[29](https://arxiv.org/html/2608.20980#bib.bib31),[4](https://arxiv.org/html/2608.20980#bib.bib30)\], a natural link to processing data spread across both space and time has emerged\[[21](https://arxiv.org/html/2608.20980#bib.bib9),[43](https://arxiv.org/html/2608.20980#bib.bib11),[20](https://arxiv.org/html/2608.20980#bib.bib14),[11](https://arxiv.org/html/2608.20980#bib.bib27),[10](https://arxiv.org/html/2608.20980#bib.bib41),[3](https://arxiv.org/html/2608.20980#bib.bib38)\]\. The common taxonomy differentiates between static graphs, i\.e\., fixed spatial structure, and dynamic graphs, i\.e\., the spatial interactions are allowed to vary in time\[[15](https://arxiv.org/html/2608.20980#bib.bib1)\]\. In the present work, we focus on static graphs equipped with time series features, representing a multivariate forecasting problem based on sparse spatial interactions given by a graph defined in the data scenario, e\.g\., traffic congestion or power\-grid load forecasting\[[5](https://arxiv.org/html/2608.20980#bib.bib13)\]\.

At the core of the evaluation protocol, multiple works investigated simple baselines given by historical averages, linear regression or kernel machine learning\[[21](https://arxiv.org/html/2608.20980#bib.bib9),[25](https://arxiv.org/html/2608.20980#bib.bib28),[44](https://arxiv.org/html/2608.20980#bib.bib29)\]\. Expected to be outperformed by non\-linear models, DLinear\[[44](https://arxiv.org/html/2608.20980#bib.bib29)\]—a trend and seasonality decomposition\-based approach—presented a baseline which was able to compete at a performance level of transformers\. In the same fashion,[13](https://arxiv.org/html/2608.20980#bib.bib40)pointed at a strong baseline for continuous time data using a state\-space motivated approach building on classical moving average and persistent forecasting\. Apart from baselines, standardized protocols play a crucial role in the reliable advancement of methods demonstrating progress on benchmark datasets\. Facilitated by popular libraries111Re\-released and available in[PyTorch Geometric Temporal](https://pytorch-geometric-temporal.readthedocs.io/)\., datasets like Chickenpox, PedalMe, WikiMaths, METR\-LA, and PEMS\-BAY\[[31](https://arxiv.org/html/2608.20980#bib.bib25),[21](https://arxiv.org/html/2608.20980#bib.bib9)\], sourced from a range of real\-world forecasting problems including internet and road traffic and public health, are frequently employed as the benchmark datasets\[[46](https://arxiv.org/html/2608.20980#bib.bib10),[12](https://arxiv.org/html/2608.20980#bib.bib26),[22](https://arxiv.org/html/2608.20980#bib.bib15),[16](https://arxiv.org/html/2608.20980#bib.bib39),[3](https://arxiv.org/html/2608.20980#bib.bib38),[5](https://arxiv.org/html/2608.20980#bib.bib13)\]\.

In particular, many studies have directly compared spatiotemporal GNNs to baseline models using these benchmark datasets\. For example,[25](https://arxiv.org/html/2608.20980#bib.bib28)evaluated a constant mean predictor and node\-wise ridge regression, which demonstrated early evidence that nonlinear models were underperforming on Chickenpox, PedalMe, and WikiMaths\. The first simple yet well\-performing approaches on METR\-LA and PEMS\-BAY were based on time and space identifiers\[[34](https://arxiv.org/html/2608.20980#bib.bib22)\]\. While[21](https://arxiv.org/html/2608.20980#bib.bib9)included ARIMA \(auto\-regressive integrated moving average\) alongside seasonal historical averages, linear support vector regression, and LSTM \(long short\-term memory\) as baseline methods, they come short in hyperparameter tuning ablations and have fallen behind nonlinear GNNs ever since\[[5](https://arxiv.org/html/2608.20980#bib.bib13)\]\. Considering the great interest in these datasets, we are taking a step back to more deeply scrutinize the data sources, their evaluation protocols, and the employed baselines\.

Main Contributions\.

- •A statistical analysis of lagged temporal and spatial correlations on all five widely adopted benchmark datasets: Chickenpox, PedalMe, WikiMaths, METR\-LA, and PEMS\-BAY\.
- •An analysis of existing evaluation protocols, revealing important shortcomings such as the running of models and comparisons on differenced rather than raw time series, which yields inferior predictive performance across the board\.
- •A re\-evaluation of baselines, motivating the use of spatially uninformed models such as SARIMA \(seasonal ARIMA\) on Chickenpox and WikiMaths\.
- •Motivated by these re\-evaluations, we propose training of GNNs on SARIMA residuals as a training target that can raise predictive performance to the level of state\-of\-the\-art methods on METR\-LA and PEMS\-BAY\.
- •A synthetic experiment showing how performance increases from linear models come partially from better adaptiveness to heterogeneity in node\-wise time series properties\.

## 2Benchmark Datasets and Evaluation Protocols

The data summaries for the benchmark datasets— Chickenpox, WikiMaths, PedalME, METR\-LA, PEMS\-BAY—are shown in Table[1](https://arxiv.org/html/2608.20980#S2.T1)\. We shall first perform a correlation analysis of these datasets in[Section2\.1](https://arxiv.org/html/2608.20980#S2.SS1)to understand their key properties, before then going deeper and looking at evaluation protocols for the various datasets in[Section2\.2](https://arxiv.org/html/2608.20980#S2.SS2), where in particular the Chickenpox and PedalME datasets are differenced by default in the PyTorch Geometric Temporal library \(showing weekly differences rather than absolute numbers\), and we study the impact of this on forecast method comparison\. First, however, we establish some notation for the paper\.

#### Notation\.

Let𝒢=\(𝒱,ℰ,𝐀\)\\mathcal\{G\}=\(\\mathcal\{V\},\\mathcal\{E\},\\mathbf\{A\}\)be a weighted undirected graph withN=\|𝒱\|N=\|\\mathcal\{V\}\|nodes and\|ℰ\|\|\\mathcal\{E\}\|edges, where𝐀∈ℝN×N\\mathbf\{A\}\\in\\mathbb\{R\}^\{N\\times N\}is the adjacency matrix withAi​j\>0A\_\{ij\}\>0if\(i,j\)∈ℰ\(i,j\)\\in\\mathcal\{E\}and zero otherwise\. Each nodei∈𝒱i\\in\\mathcal\{V\}is associated with a scalar observationyt\(i\)∈ℝy\_\{t\}^\{\(i\)\}\\in\\mathbb\{R\}at discrete timet∈\{1,…,T\}t\\in\\\{1,\\ldots,T\\\}\. Stacking over nodes, we write the system state at timettas𝐲t∈ℝN\\mathbf\{y\}\_\{t\}\\in\\mathbb\{R\}^\{N\}, and the full dataset as𝐘∈ℝT×N\\mathbf\{Y\}\\in\\mathbb\{R\}^\{T\\times N\}\.

Table 1:Dataset statistics for spatiotemporal forecasting benchmarks\.A forecasting task is defined as follows\. Given a history window of lengthHH, i\.e\.𝐘t−H\+1:t∈ℝH×N\\mathbf\{Y\}\_\{t\-H\+1:t\}\\in\\mathbb\{R\}^\{H\\times N\}, and the graph𝒢\\mathcal\{G\}, the goal is to predict the nextkksteps𝐘t\+1:t\+k∈ℝk×N\\mathbf\{Y\}\_\{t\+1:t\+k\}\\in\\mathbb\{R\}^\{k\\times N\}\. We denote the neighborhood of nodeiias𝒩⁡\(i\)=\{j:\(i,j\)∈ℰ\}\\mathcal\{N\}\(i\)=\\\{j:\(i,j\)\\in\\mathcal\{E\}\\\}, and write the univariate time series for nodeiias𝐲\(i\)=\(y1\(i\),…,yT\(i\)\)⊤∈ℝT\\mathbf\{y\}^\{\(i\)\}=\(y^\{\(i\)\}\_\{1\},\\ldots,y^\{\(i\)\}\_\{T\}\)^\{\\top\}\\in\\mathbb\{R\}^\{T\}\.

### 2\.1Correlation Analysis

#### Setup\.

Based on established statistical tools\[[38](https://arxiv.org/html/2608.20980#bib.bib4)\], we analyzespatialandtemporalcorrelation probed at a number of temporal lags equal to3×3\\timesthe input window in each dataset, examining whether default window sizes cover properties such as seasonality\. By spatial correlation, we refer to the linear relationship between nodeiiat timettand its neighborsj∈𝒩⁡\(i\)∖\{i\}j\\in\\mathcal\{N\}\(i\)\\setminus\\\{i\\\}at timeτ≤t\\tau\\leq t, weighted by the \(normalized\) adjacency matrixAi​jA\_\{ij\}\. To compute this we first define the weighted neighborhood signal

y¯τ\(i\)=∑j∈𝒩⁡\(i\)∖\{i\}Ai​j​yτ\(j\),\\bar\{y\}^\{\(i\)\}\_\{\\tau\}\\;=\\;\\sum\_\{j\\in\\mathcal\{N\}\(i\)\\setminus\\\{i\\\}\}A\_\{ij\}\\,y^\{\(j\)\}\_\{\\tau\},\(1\)and then compute the \(Pearson\) spatial correlation as

ρisp\(τ\)=Corr\(yt\(i\),y¯t−τ\(i\)\),τ=0,1,2,…\.\\rho^\{\\mathrm\{sp\}\}\_\{i\}\(\\tau\)\\;=\\;\\mathrm\{Corr\}\\\!\\left\(y^\{\(i\)\}\_\{t\},\\,\\bar\{y\}^\{\(i\)\}\_\{t\-\\tau\}\\right\),\\quad\\tau=0,1,2,\\ldots\.\(2\)By temporal correlation we refer to the autocorrelation of each node’s time series, i\.e\.

ρitemp\(τ\)=Corr\(yt\(i\),yt−τ\(i\)\),τ=0,1,2,…\.\\rho^\{\\mathrm\{temp\}\}\_\{i\}\(\\tau\)\\;=\\;\\mathrm\{Corr\}\\\!\\left\(y^\{\(i\)\}\_\{t\},\\,y^\{\(i\)\}\_\{t\-\\tau\}\\right\),\\quad\\tau=0,1,2,\\ldots\.\(3\)
Furthermore, we generally expect that in a spatiotemporal dataset, the state of nodeiiand its neighbours at timettshare mutual information\. To disentangle this contemporaneous spatial relationship, and inspired by the standard usage of the partial autocorrelation function\[[38](https://arxiv.org/html/2608.20980#bib.bib4)\], we perform a procedure which we callpartial spatialcorrelation, where we estimate the correlation remaining between the state of neighbouring nodes at timet−τt\-\\tauand the present state of nodeiiat timett,conditionalon the state of nodeiiat timet−τt\-\\tau\. This allows for an estimate of how predictive \(in a linear sense\) is the lagged spatial information versus a time\-only model\. Formally, lety^t\(i\)=α​yt−τ\(i\)\+β\\hat\{y\}^\{\(i\)\}\_\{t\}=\\alpha\\,y^\{\(i\)\}\_\{t\-\\tau\}\+\\betaandy¯^t−τ\(i\)=γ​yt−τ\(i\)\+δ\\hat\{\\bar\{y\}\}^\{\(i\)\}\_\{t\-\\tau\}=\\gamma\\,y^\{\(i\)\}\_\{t\-\\tau\}\+\\deltabe the least\-squares projections ofyt\(i\)y^\{\(i\)\}\_\{t\}andy¯t−τ\(i\)\\bar\{y\}^\{\(i\)\}\_\{t\-\\tau\}ontoyt−τ\(i\)y^\{\(i\)\}\_\{t\-\\tau\}, respectively, wherey¯t−τ\(i\)\\bar\{y\}^\{\(i\)\}\_\{t\-\\tau\}is the weighted neighbourhood signal at lagτ\\tauas defined in \([1](https://arxiv.org/html/2608.20980#S2.E1)\)\. Thepartial spatialcorrelation is then

ρipsp\(τ\)=Corr\(yt\(i\)−y^t\(i\),y¯t−τ\(i\)−y¯^t−τ\(i\)\),τ=1,2,…\.\\rho^\{\\mathrm\{psp\}\}\_\{i\}\(\\tau\)\\;=\\;\\mathrm\{Corr\}\\\!\\left\(y^\{\(i\)\}\_\{t\}\-\\hat\{y\}^\{\(i\)\}\_\{t\},\\;\\bar\{y\}^\{\(i\)\}\_\{t\-\\tau\}\-\\hat\{\\bar\{y\}\}^\{\(i\)\}\_\{t\-\\tau\}\\right\),\\quad\\tau=1,2,\\ldots\.\(4\)
The observed temporal, spatial, and partial spatial correlations for each of the five benchmark datasets of interest \(computed over the full time series\) are displayed in Figure[1](https://arxiv.org/html/2608.20980#S2.F1)\(Chickenpox, PedalMe, WikiMaths\) and Figure[2](https://arxiv.org/html/2608.20980#S2.F2)\(METR\-LA, PEMS\-BAY\) respectively\. For each lagτ\\tauconsidered, we report the median and 5th/95th percentiles ofρitemp​\(τ\)\\rho^\{\\mathrm\{temp\}\}\_\{i\}\(\\tau\),ρisp​\(τ\)\\rho^\{\\mathrm\{sp\}\}\_\{i\}\(\\tau\), andρipsp​\(τ\)\\rho^\{\\mathrm\{psp\}\}\_\{i\}\(\\tau\)across all nodesi∈𝒱i\\in\\mathcal\{V\}\.

#### Limitations\.

The caveat to the above metrics is that they are limited to capturing linear dependency structures\. Nonetheless, since our focus lies in re\-evaluating the power of linear baselines for the benchmark datasets in Section[3](https://arxiv.org/html/2608.20980#S3), these statistics serve as a sufficient tool to inform us about appropriate baseline models and evaluation protocols\.

![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/chickenpox_acf.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/pedalme_acf.png)\(b\)
![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/wikimaths_acf.png)\(c\)

Figure 1:Lagged temporal \(blue\), spatial \(orange\), and partial spatial \(green\) correlation of the Chickenpox \(a\), PedalMe \(b\), and WikiMaths \(c\) datasets\. The line indicates the median correlation over nodes, while the shaded region indicates the 5th and 95th percentiles\.
#### Findings \(Chickenpox\)\.

Examining Figure[1\(a\)](https://arxiv.org/html/2608.20980#S2.F1.sf1), the Chickenpox dataset shows very weak correlation across all lags considered, except for the lag 1 temporal correlation, which is strongly negative\. In time series analysis, this pattern can often indicateoverdifferencing— differencing a time series more than is necessary to achieve stationarity\. Indeed, this dataset has been released in the PyTorch Geometric Temporal library by default as a differenced time series in the library \(giving the difference in weekly reported cases rather than the absolute number\), the impact of which we will investigate further in Section[2\.2](https://arxiv.org/html/2608.20980#S2.SS2)\. The median partial spatial correlation is very close to zero at all lags, indicating almost no spatial information\.

#### Findings \(PedalMe\)\.

As expected from its small number of timesteps, the correlation estimates are very volatile, as presented in Figure[1\(b\)](https://arxiv.org/html/2608.20980#S2.F1.sf2)\. It is therefore difficult to assess with certainty any consistent properties\. As with Chickenpox, the data does appear to have been overdifferenced, considering the observed negative temporal correlation at lag 1, and indeed this dataset has also been released in the PyTorch Geometric Temporal library to show weekly differences in demand rather than absolute numbers\. Unlike with Chickenpox, however, it appears to potentially have more spatial information remaining after controlling for time, though it is difficult to make such a claim with certainty\. In particular, when applying the Ljung\-Box test\[[24](https://arxiv.org/html/2608.20980#bib.bib42)\]on the node\-wise spatial and temporal correlation functions, the test failed to distinguish them from white noise\.

#### Findings \(WikiMaths\)\.

In Figure[1\(c\)](https://arxiv.org/html/2608.20980#S2.F1.sf3)observe a clear weekly \(lag 7\) cycle\. Considering the dataset’s loader comes with a size of 8 time steps, this means the default settings of the loader can effectively cover the seasonality that is present\. Even after controlling for time, the partial spatial correlation function still shows some amount of spatial information also on a weekly cycle\.

#### Findings \(METR\-LA, PEMS\-BAY\)\.

Examining the traffic datasests in Figure[2](https://arxiv.org/html/2608.20980#S2.F2), we see a very different picture\. Due to the high sampling rate of the data \(5 minutes\), we can see a strong, slowly decaying process with a long memory in both the temporal and spatial correlations\. However, when accounting for temporal correlation, the partial spatial correlation becomes much smaller\. Yet it also remains very stable over its median, decaying slowly in both datasets\.

Overall, the traffic datasets show the strongest evidence of residual spatial information \(in a linear sense\)\. Despite that, the benchmarks are, in general, characterized by much weaker spatial than temporal correlation\. Since a node’s own temporal correlation is a much more compact and stable target, this might suggest that over these benchmarks, models with higher complexity and more weight on the time axis could be the most successful\.

![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/metrla_acf.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/pemsbay_acf.png)\(b\)

Figure 2:Lagged temporal \(blue\), spatial \(orange\), and partial spatial \(green\) correlation of the METR\-LA \(a\), and PEMSBAY \(b\) datasets\. The line indicates the median correlation over nodes, while the shaded region indicates the 5th and 95th percentiles\.

### 2\.2Should models be run on differenced data?

Differencing is often performed on time series to make the series more stationary and potentially easier to model\[[17](https://arxiv.org/html/2608.20980#bib.bib2)\]\. First order differencing, i\.e\. the operationZt=∇Yt=Yt−Yt−1Z\_\{t\}=\\nabla Y\_\{t\}=Y\_\{t\}\-Y\_\{t\-1\}, removes a linear trend\. However, differencing more than is necessary \(overdifferencing\) can instead make it much harder to fit models for forecasting\. This is due to the fact that differencing acts as a high\-pass filter, causing meaningful low frequency dynamics to dampen, and high frequency noise to strengthen; effectively removing information from the data\. To see why, we present an adapted proof from[38](https://arxiv.org/html/2608.20980#bib.bib4)which considers how the spectrum of a differenced time series compares to the original time series:

###### First order differencing is a high pass filter\.

Let the time series of interest beYtY\_\{t\}, and apply the first order differencing operatorZt=∇YtZ\_\{t\}=\\nabla Y\_\{t\}\. We will assume thatYtY\_\{t\}is weakly stationary and admits a spectrumfy​\(ω\)f\_\{y\}\(\\omega\)\. If we define the spectrum ofZtZ\_\{t\}asfz​\(ω\)f\_\{z\}\(\\omega\), andA⁡\(ω\)A\(\\omega\)as the frequency response of a linear filter, it is a well known result that the spectrum of the output of a linear filter isfz​\(ω\)=\|A⁡\(ω\)\|2f\_\{z\}\(\\omega\)=\|A\(\\omega\)\|^\{2\}fy​\(ω\)f\_\{y\}\(\\omega\)\. The frequency response of the operationZt=∇YtZ\_\{t\}=\\nabla Y\_\{t\}isA⁡\(ω\)=1−exp⁡\(−i​ω\)A\(\\omega\)=1\-\\exp\(\-i\\omega\)\. Therefore\|1−exp⁡\(−2​π​i​ω\)\|2=2−2​cos⁡\(ω\)=4​sin2⁡\(ω/2\)\|1\-\\exp\(\-2\\pi i\\omega\)\|^\{2\}=2\-2\\cos\(\\omega\)=4\\sin^\{2\}\(\\omega/2\)\. Hence for smallω\\omega,fz​\(ω\)≈ω2​fy​\(ω\)f\_\{z\}\(\\omega\)\\approx\\omega^\{2\}f\_\{y\}\(\\omega\), leading to greater dampening asω\\omegaapproaches00\. For largeω\\omegaapproaching the Nyquist frequency,π\\pi, thenfz​\(ω\)≈4​fy​\(ω\)f\_\{z\}\(\\omega\)\\approx 4f\_\{y\}\(\\omega\), and the filter will relatively amplify high\-frequency components\. ∎

It follows then, if measurement error and noise is high\-frequency concentrated, and the true predictive signal of interest is low\-frequency concentrated, that models trained on overdifferenced datasets will in general tend to overfit or misfit the data, concentrating on the noise instead of the signal\.

As discussed already, the Chickenpox dataset released from the PyTorch Geometric Temporal library is differenced by default\. If we instead examine the correlations of the undifferenced dataset222available at[https://archive\.ics\.uci\.edu/dataset/580/hungarian\+chickenpox\+cases](https://archive.ics.uci.edu/dataset/580/hungarian+chickenpox+cases)in Figure[3](https://arxiv.org/html/2608.20980#S2.F3), we can see that differencing has substantially changed the correlation properties of the Chickenpox dataset\. By undoing the differencing, the temporal correlation function has changed from being non\-existent \(as in Figure[1](https://arxiv.org/html/2608.20980#S2.F1)\), to being a persistent function with correlations lasting until at least lag 10\. Furthermore, its spatial signal also gained much more linear informativeness, remaining much larger even after computing partial spatial correlations\. The PedalMe dataset remains noisy after undifferencing due to its low sample size, though its lag 1 and lag 2 temporal correlation function’s median appears much more stable, albeit with large intervals around the median\.

Given the above rationale, we hypothesize that training on the undifferenced chickenpox dataset should lead to models that overfit less, since it will have much more interpretable signal around its low frequency components, compared to the differenced dataset, where noise is amplified\. To test this hypothesis, we ran a suite of one\-step ahead forecasting models on the undifferenced \(raw\) and differenced \(loader\) datasets, the results of which are displayed in Appendix Table[5](https://arxiv.org/html/2608.20980#A2.T5)\. We can indeed see that the ratio of mean squared errorstestMSE\\text\{test\}\_\{\\text\{MSE\}\}totrainMSE\{\\text\{train\}\_\{\\text\{MSE\}\}\}becomes smaller on the undifferenced \(raw\) chickenpox dataset, indicating better generalization of models\. Furthermore, all models trained on the undifferenced \(raw\) dataset show improved performance on the test set even after their predictions are differenced \(denotedtestMSEDS\{\\text\{test\}^\{\\text\{DS\}\}\_\{\\text\{MSE\}\}\}\) when compared to the models that were trained on differenced \(loader\) datasets\. This identifies a significant shortcoming of many existing evaluation protocols in the literature which have trained and tested models on differenced data\.

![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/raw_chickenpox_acf.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/raw_pedalme_acf.png)\(b\)

Figure 3:Lagged temporal \(blue\), spatial \(orange\), and partial spatial \(green\) correlation of the undifferenced Chickenpox \(a\), and PedalMe \(b\) datasets\. The line indicates the median correlation over nodes, while the shaded region indicates the 5th and 95th percentiles\.

## 3Spatially Uninformed Baselines

Motivated by our analyses in Section[2\.1](https://arxiv.org/html/2608.20980#S2.SS1)revealing that the analyzed spatio\-temporal graph benchmark datasets are strongly driven by the temporal signals, we now investigate simple temporal models as spatially unaware baselines\. As will be demonstrated, even simple linear models have very strong performance, in some cases almost reaching state of the art\. Although some simple linear benchmarks, such as historical average and linear regression for Chickenpox, PedalMe, and WikiMaths\[[25](https://arxiv.org/html/2608.20980#bib.bib28)\]or seasonal averages and ARIMA \(auto\-regressive integrated moving average\) for METR\-LA and PEMS\-BAY\[[21](https://arxiv.org/html/2608.20980#bib.bib9)\], have previously been used in the literature, we now investigate the potential for improved baseline formulation\. Specifically, we noticed the following shortcomings: \(i\) inconsistent context windows across baselines; \(ii\) fixed hyperparameters for ARIMA\(pp,dd,qq\) for all nodes; and \(iii\) no GPU\-optimized backends supporting node\-wise fitting such as AutoARIMA\[[18](https://arxiv.org/html/2608.20980#bib.bib7)\]or Cuml\[[30](https://arxiv.org/html/2608.20980#bib.bib8)\]\.

In Section[3\.1](https://arxiv.org/html/2608.20980#S3.SS1)we shall revisit some of the previously used linear baselines in the literature on three smaller benchmarks\. We suggest SARIMA as an improved linear benchmark in Section[3\.2](https://arxiv.org/html/2608.20980#S3.SS2), showing it tends to outperform many larger models on smaller datasets\. By combining SARIMA with a spatial GNN model on the larger traffic datasets METR\-LA and PEMS\-BAY, we reach the level of state\-of\-the\-art in Section[3\.3](https://arxiv.org/html/2608.20980#S3.SS3)\.

### 3\.1Naive Baselines

Motivated by prior works pointing to simple baselines such as persistent forecasting\[[13](https://arxiv.org/html/2608.20980#bib.bib40)\]or DLinear\[[44](https://arxiv.org/html/2608.20980#bib.bib29)\], we compile the following list of baselines on Chickenpox, PedalMe and WikiMaths:

- •Persistence\.𝐲^t\+1=𝐲t\\hat\{\\mathbf\{y\}\}\_\{t\+1\}=\\mathbf\{y\}\_\{t\}\. Reuses the most recent observation per node; the simplest reference whenever the signal is persistent\.
- •Historical average\.Per\-node mean of the target conditioned on a season bin: \(snapshot indexmodP\\bmod\\ P\) for the social/web datasets \(P=7P=7for WikiMaths,P=52P=52for Chickenpox/PedalMe\)\. Captures pure periodic structure with no temporal regression\.
- •AR\(H\)\(H\)per node\.A separate ridge\-regularized autoregressor of orderHHper nodeii, fit in closed form:y^\(i\)t\+1=𝐰i⊤𝐲\(i\)t−H\+1:t\+bi\\hat\{y\}^\{\(i\)\}\_\{t\+1\}=\\mathbf\{w\}\_\{i\}^\{\\top\}\\mathbf\{y\}^\{\(i\)\}\_\{t\-H\+1:t\}\+b\_\{i\}\. This is a spatially\-unaware temporal model\.
- •RidgeVAR\.A single ridge regression that maps the flattened lag matrix𝐘t−H\+1:t∈ℝH×N\\mathbf\{Y\}\_\{t\-H\+1:t\}\\in\\mathbb\{R\}^\{H\\times N\}of all nodes to all next\-step targets𝐲t\+1∈ℝN\\mathbf\{y\}\_\{t\+1\}\\in\\mathbb\{R\}^\{N\}\. Captures arbitrary linear cross\-node dynamics; effectively a fully connected linear graph filter\.
- •DLinear\.A trend/seasonal decomposition baseline\[[44](https://arxiv.org/html/2608.20980#bib.bib29)\]: the input series is split into a moving\-average trend and the residual seasonal component, each mapped to the next step by a shared linear layer\. Trained with Adam \(MSE loss\)\.

#### Results\.

As shown in Table[2](https://arxiv.org/html/2608.20980#S3.T2), among the best baselines for Chickenpox and PedalMe are AR\(HH\) per node and DLinear, despite them being spatially unaware\. With WikiMaths, RidgeVAR is coming first and even placing above the state\-of\-the\-art TDE\-GNN\[[11](https://arxiv.org/html/2608.20980#bib.bib27)\]\.

Table 2:Linear baseline test MSE on Chickenpox, WikiMaths, and PedalMe compared to state\-of\-the\-art TDE\-GNN\[[11](https://arxiv.org/html/2608.20980#bib.bib27)\]\. We evaluate one\-step\-ahead prediction on the standard 90/10 train/test split\. Best per dataset in bold\.

### 3\.2SARIMA

SARIMA is a well known class of linear time series models that can predict future time stepsyt\+hy\_\{t\+h\}as a seasonally\-informed linear combination of past observationsyτy\_\{\\tau\},τ≤t\\tau\\leq t, and past noise innovationsϵτ\\epsilon\_\{\\tau\}\. After naive linear forecasts, they form the next step up in complexity\. The order of a SARIMA decides how many lags and coefficients it needs, acting as hyperparameters, after which the coefficients can be fitted using maximum likelihood optimization\. Further background details on SARIMA, including our implementation procedures on the benchmark datasets, can be found in Appendix[A](https://arxiv.org/html/2608.20980#A1)\.

#### Results\.

In Table[3](https://arxiv.org/html/2608.20980#S3.T3)we include SARIMA as a baseline in comparison with various state\-of\-the\-art models\. On WikiMaths then SARIMA outperforms existing state\-of\-the\-art models\. With Chickenpox, SARIMA is only outperformed by TDE\-GNN, perhaps because of the much weaker temporal signal on the differenced Chickenpox dataset\. However, ARIMA with temporal encodings is able to outperform TDE\-GNN\. Only on PedalMe do ARIMA\-based methods do poorly\. When we diagnosed the fitted model orders, we found that due to the small number of samples, ARIMAs were not able to fit much more than a single MA term\. On traffic data \(shown in Table[4](https://arxiv.org/html/2608.20980#A1.T4)\), the fixed order SARIMA does not perform that well in comparison to many GNNs\. In part this can be explained by the large size of the time series, the imposed uniform order on the SARIMA, and the stable spatial signal present in the time series\. However, it still is able to beat the previously best linear benchmark, DLinear\.

### 3\.3Predicting SARIMA residuals is a meaningful task

Since many of the benchmarks have strong temporal correlation, and we saw how effective SARIMAs are at capturing that temporal correlation, we hypothesize that residual errors from a SARIMA might create a better target for the training of GNNs, as, the lessened temporal information remaining in the data could potentially allow the models to focus more on the remnant graph parts of the data\.

#### Setup\.

For each benchmark, SARIMA training is done with the aforementioned appropriate protocol\. For each time steptt, we add an additional feature channel that includes past observed SARIMA errorsy^t−τ=yt−τ,pred−yt−τ,true\\hat\{y\}\_\{t\-\\tau\}=y\_\{t\-\\tau,\\text\{pred\}\}\-y\_\{t\-\\tau,\\text\{true\}\}forτ∈\[1,H\]\\tau\\in\[1,H\], set the SARIMA erroryt−τy\_\{t\-\\tau\}as the training target, and addyt\+1,predy\_\{t\+1,\\text\{pred\}\}as an additional covariate\.

This augmented data is then stacked onto the original dataset, and fed into a follow up model that can better use spatial information to correct SARIMA errors\. In our case, we mostly used GCRN\-GRU\. For comparison with a spatially unaware model, we used LSTMs\.

#### Results\.

On Chickenpox \(Table[3](https://arxiv.org/html/2608.20980#S3.T3)\), this residual learning model was able to marginally improve upon the SARIMA with time encoded Fourier features\. Despite our hypothesis that it would be spatial features driving this progress, on Chickenpox the pure time series SARIMA \+ LSTM still has the best performance\. On PedalMe we can see that the models do not improve that much upon previous results, with SARIMA GCRN\-GRU actually decreasing performance\. On WikiMaths these models were able to improve significantly upon vanilla SARIMA models, further pushing the state of the art on the benchmark\.

The best improvements, however, came from the residual learning model on METR\-LA and PEMSBAY \(Table[4](https://arxiv.org/html/2608.20980#A1.T4)\)\. Perhaps because the signal is spatially richer, the residual forecaster is able to exceed previously reported metrics to become state of the art on both traffic datasets\.

Table 3:Average test MSE and standard deviation \(↓\\downarrow\) of 10 experimental repetitions on Chickenpox, PedalMe, and WikiMaths\. Data split \(9:1\)\. Baseline results are reported from[31](https://arxiv.org/html/2608.20980#bib.bib25),[12](https://arxiv.org/html/2608.20980#bib.bib26),[11](https://arxiv.org/html/2608.20980#bib.bib27)\.Bold= best;underline= second best;italic= third best per dataset\.

## 4Outlook

### 4\.1Suggestions for benchmarks and baselines

In our analysis in Section[2\.1](https://arxiv.org/html/2608.20980#S2.SS1), we demonstrated a repeatable procedure to analyze the respective spatial and temporal properties of benchmarks\. Based on that analysis, we recommended that the PedalMe and Chickenpox datasets should remain undifferenced when used for benchmarking purposes\.

Additionally, in Section[3\.2](https://arxiv.org/html/2608.20980#S3.SS2)we showed how SARIMA models can form effective baseline models even in the spatiotemporal graph regression setting\. Furthermore, finding stronger baseline models can inspire new methods, as we have done in Section[3\.3](https://arxiv.org/html/2608.20980#S3.SS3)by integrating SARIMA with GNNs\.

### 4\.2What makes linear models so effective on space\-time graphs?

A partial answer can be found in how well linear models and GNNs adapt to heterogeneity in nodal time series properties\. It has been previously pointed out in the node\-classification literature that GNNs might not do so well with heterophily\[[47](https://arxiv.org/html/2608.20980#bib.bib6)\]\. Given this, we hypothesize that part of what makes ARIMAs and other linear models so effective at a cheap price is that they do not implicitly assume that the data they are modelling have a homogenous spatiotemporal response to the past signal, unlike GNNs, which might perhaps do, even if outside of node classification tasks\.

To test this hypothesis, we set up a synthetic experiment\. Consider the following discrete random graph process:

Xi,t=ai​Xi,t−1\+bi​λ​∑j≠iwi,j​Xj,t−1\+ϵi,t−1\.X\_\{i,t\}=a\_\{i\}X\_\{i,t\-1\}\+b\_\{i\}\\lambda\\sum\_\{j\\neq i\}w\_\{i,j\}X\_\{j,t\-1\}\+\\epsilon\_\{i,t\-1\}\.\(5\)In other words, an AR\(1\) process with a graph aggregating operation\. The coefficientaia\_\{i\}controls the temporal memory of the process,bib\_\{i\}the node\-specific spatial dependence,λ\\lambdais the global graph spatial weighting\.

After simulating500500time steps with1010nodes, we separately fit a GRU\-GCN\[[32](https://arxiv.org/html/2608.20980#bib.bib24)\]and ARIMA on the space time graphXXwith a 90/10 split\. We repeat this procedure, while varying the variance ofaa,bb, andww, and changing the size ofλ\\lambda\. For full details, see Appendix[C](https://arxiv.org/html/2608.20980#A3)\. We perform a regression with test MSE gap between GRUGCN and ARIMA as the regression target\. The dependent variables are: variance node temporal sensitivity ,Var⁡\(\{ai\}i=1N\)\\operatorname\{Var\}\\left\(\\\{a\_\{i\}\\\}\_\{i=1\}^\{N\}\\right\); variance of node spatial sensitivity,Var⁡\(\{bi\}i=1N\)\\operatorname\{Var\}\\left\(\\\{b\_\{i\}\\\}\_\{i=1\}^\{N\}\\right\); the average variance of node neighbour weights,1N​∑iVar⁡\(\{wi,j\}j=1N\)\\frac\{1\}\{N\}\\sum\_\{i\}\\operatorname\{Var\}\\left\(\\\{w\_\{i,j\}\\\}\_\{j=1\}^\{N\}\\right\); and the global spatial sensitivityλ\\lambda\.

As we would expect, the regression results reveal that the variance of the temporal coefficientsVar⁡\(\{ai\}i=1N\)\\operatorname\{Var\}\\left\(\\\{a\_\{i\}\\\}\_\{i=1\}^\{N\}\\right\)has a positive influence on the quantityGRUGCNMSE−ARIMAMSE\\text\{GRUGCN\}\_\{\\text\{MSE\}\}\-\\text\{ARIMA\}\_\{\\text\{MSE\}\}\(see Table[6](https://arxiv.org/html/2608.20980#A3.T6)\)\. On the other hand, the variance of the spatial coefficients was not found to be significant in the regression\. Most of the spatial variation seems to be contained within changes inλ\\lambdaand changes in the mean nodal variance\. The global spatial coefficient also seems to be the only parameter that has a positive effect on GNN performance, perhaps because it is global, and hence homogeneous\. Despite being spatial, variance in the spatial neighbor weights seems to hinder the GNN more than help it\.

## 5Conclusion

In this work, we analyze datasets commonly used for benchmarking spatiotemporal regression—Chickenpox, PedalMe, WikiMaths, METR\-LA, and PEMS\-BAY— by means of a statistical correlation analysis, revealing most of the datasets are dominated by node\-wise temporal information\. We find that the official data\-loader from the PyTorch Geometric Temporal library releases first\-order differenced data for Chickenpox and PedalMe, which obscures the low frequency signal and weakens benchmark quality\. By comparing to the performance of spatially uninformed linear models, notably SARIMA models, we uncover a significant performance gap of state\-of\-the\-art GNNs, and show how targetting the residuals of linear models based on SARIMA can improve performance of GNN models\.

## Acknowledgments

The work of KM and AMS was supported by the Engineering and Physical Sciences Research Council \(EPSRC\) \[grant number EP/Y03533X/1\]\. SH and AF are supported by DFG project TRR 391 Spatio\-temporal Statistics for the Transition of Energy and Transport \(project 520388526\)\. ME acknowledges support from the Israeli Ministry of Innovation, Science & Technology\.

## References

- Baiet al\.\(2021\)J\. Bai, J\. Zhu, Y\. Song, L\. Zhao, Z\. Hou, R\. Du, and H\. LiA3t\-gcn: attention temporal graph convolutional network for traffic forecasting\.ISPRS International Journal of Geo\-Information10\(7\),pp\. 485\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1)\.
- Baiet al\.\(2020\)L\. Bai, L\. Yao, C\. Li, X\. Wang, and C\. WangAdaptive graph convolutional recurrent network for traffic forecasting\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.33,pp\. 17804–17815\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.7.1),[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1)\.
- Ceniet al\.\(2025\)A\. Ceni, A\. Gravina, C\. Gallicchio, D\. Bacciu, C\. Schonlieb, and M\. EliasofMessage\-passing state\-space models: improving graph learning with modern sequence modeling\.arXiv preprint arXiv:2505\.18728\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1),[§1](https://arxiv.org/html/2608.20980#S1.p2.1)\.
- Choiet al\.\(2023\)J\. Choi, S\. Hong, N\. Park, and S\. ChoGREAD: graph neural reaction\-diffusion networks\.InProceedings of the 40th International Conference on Machine Learning,A\. Krause, E\. Brunskill, K\. Cho, B\. Engelhardt, S\. Sabato, and J\. Scarlett \(Eds\.\),Proceedings of Machine Learning Research, Vol\.202,pp\. 5722–5747\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Corradiniet al\.\(2025\)F\. Corradini, F\. Gerosa, M\. Gori, C\. Lucheroni, M\. Piangerelli, and M\. ZannottiA systematic literature review of spatio\-temporal graph neural network models for time series forecasting and classification\.Neural Networks,pp\. 108269\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1),[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[§1](https://arxiv.org/html/2608.20980#S1.p3.1)\.
- Corsoet al\.\(2024\)G\. Corso, H\. Stark, S\. Jegelka, T\. Jaakkola, and R\. BarzilayGraph neural networks\.Nature Reviews Methods Primers4\(1\),pp\. 17\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Defferrardet al\.\(2016\)M\. Defferrard, X\. Bresson, and P\. VandergheynstConvolutional neural networks on graphs with fast localized spectral filtering\.Advances in Neural Information Processing Systems29\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Denget al\.\(2023\)J\. Deng, J\. Deng, D\. Yin, R\. Jiang, and X\. SongTts\-norm: forecasting tensor time series via multi\-way normalization\.ACM Transactions on Knowledge Discovery from Data18\(1\),pp\. 1–25\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.10.1)\.
- Duet al\.\(2024\)H\. Du, J\. Wang, J\. Hui, L\. Zhang, and H\. WangDenseGNN: universal and scalable deeper graph neural networks for high\-performance property prediction in crystals and molecules\.npj Computational Materials10\(1\),pp\. 292\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Eliasofet al\.\(2025\)M\. Eliasof, A\. Gravina, A\. Ceni, C\. Gallicchio, D\. Bacciu, and C\. SchönliebGraph adaptive autoregressive moving average models\.InForty\-second International Conference on Machine Learning,Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Eliasofet al\.\(2024\)M\. Eliasof, E\. Haber, E\. Treister, and C\. SchönliebOn the temporal domain of differential equation inspired graph neural networks\.InProceedings of The 27th International Conference on Artificial Intelligence and Statistics \(AISTATS\),Vol\.238,pp\. 1792–1800\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1),[§1](https://arxiv.org/html/2608.20980#S1.p1.1),[§3\.1](https://arxiv.org/html/2608.20980#S3.SS1.SSS0.Px1.p1.1),[Table 2](https://arxiv.org/html/2608.20980#S3.T2),[Table 2](https://arxiv.org/html/2608.20980#S3.T2.4),[Table 3](https://arxiv.org/html/2608.20980#S3.T3),[Table 3](https://arxiv.org/html/2608.20980#S3.T3.7)\.
- Erricaet al\.\(2023\)F\. Errica, A\. Gravina, D\. Bacciu, and A\. MicheliHidden Markov models for temporal graph representation learning\.InProceedings of the 31st European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning \(ESANN\),Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[Table 3](https://arxiv.org/html/2608.20980#S3.T3),[Table 3](https://arxiv.org/html/2608.20980#S3.T3.7)\.
- Feldmanet al\.\(2026\)O\. Feldman, K\. S\. I\. Mantri, M\. Eliasof, and C\. BaskinRevisting node affinity prediction in temporal graphs\.InThe Fourteenth International Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[§3\.1](https://arxiv.org/html/2608.20980#S3.SS1.p1.1)\.
- Garzaet al\.\(2022\)A\. Garza, M\. Mergenthaler Canseco, C\. Challú, and K\. G\. OlivaresStatsForecast: lightning fast forecasting with statistical and econometric models\.Note:PyCon Salt Lake City, Utah, US 2022External Links:[Link](https://github.com/Nixtla/statsforecast)Cited by:[Appendix A](https://arxiv.org/html/2608.20980#A1.SS0.SSS0.Px1.p1.1)\.
- Gravina and Bacciu \(2024\)A\. Gravina and D\. BacciuDeep learning for dynamic graphs: models and benchmarks\.IEEE Transactions on Neural Networks and Learning Systems35\(9\),pp\. 11788–11801\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Gravinaet al\.\(2024\)A\. Gravina, D\. Zambon, D\. Bacciu, and C\. AlippiTemporal graph odes for irregularly\-sampled time series\.InProceedings of the Thirty\-Third International Joint Conference on Artificial Intelligence,IJCAI ’24\.External Links:ISBN 978\-1\-956792\-04\-1,[Link](https://doi.org/10.24963/ijcai.2024/445)Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p2.1)\.
- Hyndman and Athanasopoulos \(2021\)R\. J\. Hyndman and G\. AthanasopoulosForecasting: principles and practice\.Third edition edition,OTexts,Melbourne, Australia\(eng\)\.Note:Preface – 1 Getting started – 2 Time series graphics – 3 Time series decomposition – 4 Time series features – 5 The forecaster’s toolbox – 6 Judgmental forecasts – 7 Time series regression models – 8 Exponential smoothing – 9 ARIMA models – 10 Dynamic regression models – 11 Forecasting hierarchical and grouped time series – 12 Advanced forecasting methods – 13 Some practical forecasting issues – Appendix: Using R – Appendix: For instructors – Appendix: ReviewsThis open access resource is available through OTexts Includes bibliographical referencesExternal Links:ISBN 978\-0\-9875071\-3\-6Cited by:[Appendix A](https://arxiv.org/html/2608.20980#A1.SS0.SSS0.Px1.p5.1),[Definition A\.1](https://arxiv.org/html/2608.20980#A1.Thmtheorem1.p1.1),[§2\.2](https://arxiv.org/html/2608.20980#S2.SS2.p1.1)\.
- Hyndman and Khandakar \(2008\)R\. J\. Hyndman and Y\. KhandakarAutomatic time series forecasting: the forecast package for r\.Journal of Statistical Software27\(3\),pp\. 1–22\.External Links:[Link](https://www.jstatsoft.org/index.php/jss/article/view/v027i03),[Document](https://dx.doi.org/10.18637/jss.v027.i03)Cited by:[Appendix A](https://arxiv.org/html/2608.20980#A1.SS0.SSS0.Px1.p1.1),[§3](https://arxiv.org/html/2608.20980#S3.p1.1)\.
- Kipf and Welling \(2017\)T\. N\. Kipf and M\. WellingSemi\-supervised classification with graph convolutional networks\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=SJU4ayYgl)Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Liet al\.\(2024\)C\. Li, Y\. Zhao, and Z\. ZhangAMGCN: adaptive multigraph convolutional networks for traffic speed forecasting: c\. li et al\.\.Applied Intelligence54\(3\),pp\. 2594–2613\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Liet al\.\(2018\)Y\. Li, R\. Yu, C\. Shahabi, and Y\. LiuDiffusion convolutional recurrent neural network: data\-driven traffic forecasting\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=SJiHXGWAZ)Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.4.1),[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1),[§1](https://arxiv.org/html/2608.20980#S1.p1.1),[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[§1](https://arxiv.org/html/2608.20980#S1.p3.1),[§3](https://arxiv.org/html/2608.20980#S3.p1.1)\.
- Liuet al\.\(2023a\)H\. Liu, Z\. Dong, R\. Jiang, J\. Deng, J\. Deng, Q\. Chen, and X\. SongSpatio\-temporal adaptive embedding makes vanilla transformer sota for traffic forecasting\.InProceedings of the 32nd ACM international conference on information and knowledge management,pp\. 4125–4129\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p2.1)\.
- Liuet al\.\(2023b\)H\. Liu, Z\. Dong, R\. Jiang, J\. Deng, J\. Deng, Q\. Chen, and X\. SongSpatio\-temporal adaptive embedding makes vanilla transformer SOTA for traffic forecasting\.InProceedings of the 32nd ACM International Conference on Information and Knowledge Management \(CIKM\),pp\. 4125–4129\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.14.1)\.
- Ljung and Box \(1978\)G\. M\. Ljung and G\. E\. BoxOn a measure of lack of fit in time series models\.Biometrika65\(2\),pp\. 297–303\.Cited by:[§2\.1](https://arxiv.org/html/2608.20980#S2.SS1.SSS0.Px4.p1.1)\.
- Micheli and Tortorella \(2022\)A\. Micheli and D\. TortorellaDiscrete\-time dynamic graph echo state networks\.Neurocomputing496,pp\. 85–95\.External Links:ISSN 0925\-2312,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.neucom.2022.05.001),[Link](https://www.sciencedirect.com/science/article/pii/S092523122200532X)Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[§1](https://arxiv.org/html/2608.20980#S1.p3.1),[§3](https://arxiv.org/html/2608.20980#S3.p1.1)\.
- Ortegaet al\.\(2018\)A\. Ortega, P\. Frossard, J\. Kovačević, J\. M\. F\. Moura, and P\. VandergheynstGraph signal processing: overview, challenges, and applications\.Proceedings of the IEEE106\(5\),pp\. 808–828\.External Links:[Document](https://dx.doi.org/10.1109/JPROC.2018.2820126)Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Panagopouloset al\.\(2021\)G\. Panagopoulos, G\. Nikolentzos, and M\. VazirgiannisTransfer graph neural networks for pandemic forecasting\.InProceedings of the AAAI conference on artificial intelligence,Vol\.35,pp\. 4838–4845\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1)\.
- Parejaet al\.\(2020\)A\. Pareja, G\. Domeniconi, J\. Chen, T\. Ma, T\. Suzumura, H\. Kanezashi, T\. Kaler, T\. Schardl, and C\. LeisersonEvolvegcn: evolving graph convolutional networks for dynamic graphs\.InProceedings of the AAAI conference on artificial intelligence,Vol\.34,pp\. 5363–5370\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1)\.
- Poliet al\.\(2019\)M\. Poli, S\. Massaroli, J\. Park, A\. Yamashita, H\. Asama, and J\. ParkGraph neural ordinary differential equations\.arXiv preprint arXiv:1911\.07532\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Raschkaet al\.\(2020\)S\. Raschka, J\. Patterson, and C\. NoletMachine learning in python: main developments and technology trends in data science, machine learning, and artificial intelligence\.arXiv preprint arXiv:2002\.04803\.Cited by:[Appendix A](https://arxiv.org/html/2608.20980#A1.SS0.SSS0.Px1.p1.1),[§3](https://arxiv.org/html/2608.20980#S3.p1.1)\.
- Rozemberczkiet al\.\(2021\)B\. Rozemberczki, P\. Scherer, Y\. He, G\. Panagopoulos, A\. Riedel, M\. Astefanoaei, O\. Kiss, F\. Beres, G\. Lopez, N\. Collignon, and R\. SarkarPyTorch Geometric Temporal: spatiotemporal signal processing with neural machine learning models\.InProceedings of the 30th ACM International Conference on Information and Knowledge Management \(CIKM\),pp\. 4564–4573\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1),[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[Table 3](https://arxiv.org/html/2608.20980#S3.T3),[Table 3](https://arxiv.org/html/2608.20980#S3.T3.7)\.
- Seoet al\.\(2018\)Y\. Seo, M\. Defferrard, P\. Vandergheynst, and X\. BressonStructured sequence modeling with graph convolutional recurrent networks\.InInternational Conference on Neural Information Processing,pp\. 362–373\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1),[§4\.2](https://arxiv.org/html/2608.20980#S4.SS2.p3.1)\.
- Shanget al\.\(2021\)C\. Shang, J\. Chen, and J\. BiDiscrete graph structure learning for forecasting multiple time series\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.9.1)\.
- Shaoet al\.\(2022a\)Z\. Shao, Z\. Zhang, F\. Wang, W\. Wei, and Y\. XuSpatial\-temporal identity: a simple yet effective baseline for multivariate time series forecasting\.InProceedings of the 31st ACM International Conference on Information and Knowledge Management \(CIKM\),pp\. 4454–4458\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.12.1),[§1](https://arxiv.org/html/2608.20980#S1.p3.1)\.
- Shaoet al\.\(2022b\)Z\. Shao, Z\. Zhang, F\. Wang, and Y\. XuPre\-training enhanced spatial\-temporal graph neural network for multivariate time series forecasting\.InProceedings of the 28th ACM SIGKDD conference on knowledge discovery and data mining,pp\. 1567–1577\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.13.1)\.
- Shaoet al\.\(2022c\)Z\. Shao, Z\. Zhang, W\. Wei, F\. Wang, Y\. Xu, X\. Cao, and C\. S\. JensenDecoupled dynamic spatial\-temporal graph neural network for traffic forecasting\.InProceedings of the VLDB Endowment,Vol\.15,pp\. 2733–2746\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.11.1)\.
- Shenet al\.\(2022\)Y\. Shen, J\. Zhang, S\. Song, and K\. B\. LetaiefGraph neural networks for wireless communications: from theory to practice\.IEEE Transactions on Wireless Communications22\(5\),pp\. 3554–3569\.Cited by:[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Shumway and Stoffer \(2025\)R\. H\. Shumway and D\. S\. StofferTime Series Analysis and Its Applications: With R Examples\.Springer Texts in Statistics,Springer Nature Switzerland,Cham\(en\)\.External Links:ISBN 978\-3\-031\-70583\-0 978\-3\-031\-70584\-7,[Link](https://link.springer.com/10.1007/978-3-031-70584-7),[Document](https://dx.doi.org/10.1007/978-3-031-70584-7)Cited by:[§2\.1](https://arxiv.org/html/2608.20980#S2.SS1.SSS0.Px1.p1.1),[§2\.1](https://arxiv.org/html/2608.20980#S2.SS1.SSS0.Px1.p2.1),[§2\.2](https://arxiv.org/html/2608.20980#S2.SS2.p1.1)\.
- Taheri and Berger\-Wolf \(2019\)A\. Taheri and T\. Berger\-WolfPredictive temporal embedding of dynamic graphs\.InProceedings of the 2019 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining,pp\. 57–64\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1)\.
- Tsay \(2010\)R\. S\. TsayAnalysis of financial time series\.3rd ed edition,Wiley series in probability and statistics,Wiley,Cambridge, Mass\.Note:Financial time series and their characteristics – Linear time series – Volatility models – Nonlinear models and their applicationsExternal Links:ISBN 978\-0\-470\-41435\-4Cited by:[Appendix C](https://arxiv.org/html/2608.20980#A3.p2.2)\.
- Wuet al\.\(2020\)Z\. Wu, S\. Pan, G\. Long, J\. Jiang, X\. Chang, and C\. ZhangConnecting the dots: multivariate time series forecasting with graph neural networks\.InProceedings of the 26th ACM SIGKDD Conference on Knowledge Discovery and Data Mining \(KDD\),pp\. 753–763\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.8.1)\.
- Wuet al\.\(2019\)Z\. Wu, S\. Pan, G\. Long, J\. Jiang, and C\. ZhangGraph wavenet for deep spatial\-temporal graph modeling\.InProceedings of the 28th International Joint Conference on Artificial Intelligence \(IJCAI\),pp\. 1907–1913\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.5.1)\.
- Yuet al\.\(2018\)B\. Yu, H\. Yin, and Z\. ZhuSpatio\-temporal graph convolutional networks: a deep learning framework for traffic forecasting\.InProceedings of the 27th International Joint Conference on Artificial Intelligence,IJCAI’18,pp\. 3634–3640\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.3.1),[§1](https://arxiv.org/html/2608.20980#S1.p1.1)\.
- Zenget al\.\(2023\)A\. Zeng, M\. Chen, L\. Zhang, and Q\. XuAre transformers effective for time series forecasting?\.InProceedings of the AAAI conference on artificial intelligence,Vol\.37,pp\. 11121–11128\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.15.1),[§1](https://arxiv.org/html/2608.20980#S1.p2.1),[5th item](https://arxiv.org/html/2608.20980#S3.I1.i5.p1.1),[§3\.1](https://arxiv.org/html/2608.20980#S3.SS1.p1.1)\.
- Zhaoet al\.\(2019\)L\. Zhao, Y\. Song, C\. Zhang, Y\. Liu, P\. Wang, T\. Lin, M\. Deng, and H\. LiT\-gcn: a temporal graph convolutional network for traffic prediction\.IEEE transactions on intelligent transportation systems21\(9\),pp\. 3848–3858\.Cited by:[Appendix B](https://arxiv.org/html/2608.20980#A2.SS0.SSS0.Px1.p1.1)\.
- Zhenget al\.\(2020\)C\. Zheng, X\. Fan, C\. Wang, and J\. QiGman: a graph multi\-attention network for traffic prediction\.InProceedings of the AAAI conference on artificial intelligence,Vol\.34,pp\. 1234–1241\.Cited by:[Table 4](https://arxiv.org/html/2608.20980#A1.T4.8.6.1),[§1](https://arxiv.org/html/2608.20980#S1.p2.1)\.
- Zhuet al\.\(2020\)J\. Zhu, Y\. Yan, L\. Zhao, M\. Heimann, L\. Akoglu, and D\. KoutraBeyond Homophily in Graph Neural Networks: Current Limitations and Effective Designs\.arXiv\.Note:arXiv:2006\.11468 \[cs\]Comment: Accepted to NeurIPS 2020; version with full appendixExternal Links:[Link](http://arxiv.org/abs/2006.11468),[Document](https://dx.doi.org/10.48550/arXiv.2006.11468)Cited by:[§4\.2](https://arxiv.org/html/2608.20980#S4.SS2.p1.1)\.

## Appendix ABackground details on SARIMA

Consider a time series of interest, observed up to a fixed timett:\{yτ\}τ=1t\\\{y\_\{\\tau\}\\\}\_\{\\tau=1\}^\{t\},τ∈ℤ\\tau\\in\\mathbb\{Z\}\. Additionally let\{ϵτ\}τ=1t\\\{\\epsilon\_\{\\tau\}\\\}\_\{\\tau=1\}^\{t\}denote the series of noise innovations, which is, as a modelling assumption, defined to be a zero\-mean white noise process\. For our purposes, we are interested in predicting the next valueyt\+1y\_\{t\+1\}as a linear combination of previous values ofyty\_\{t\}and previous noise innovationsϵt\\epsilon\_\{t\}\.

###### Definition A\.1\(Backshift and differencing operators\)\.

Following notation from[17](https://arxiv.org/html/2608.20980#bib.bib2), we define the backshift operatorBnB^\{n\}asBn​yt=yt−nB^\{n\}y\_\{t\}=y\_\{t\-n\}, and thedd\-th order differencing operator∇d\\nabla^\{d\}as\(1−B\)d\\left\(1\-B\\right\)^\{d\}\.

The order of an ARIMA is defined byppas the autoregressive \(AR\) order,ddas the differencing order, andqqas the moving average \(MA\) order\. If we defineϕj\\phi\_\{j\}as thejj\-th autoregressive coefficient, andθj\\theta\_\{j\}as thejj\-th moving average coefficient, then we can define the AR and MA polynomials as

ϕ⁡\(B\)=1−∑j=1pϕj​Bj,θ⁡\(B\)=1\+∑j=1qθj​Bj\.\\displaystyle\\phi\(B\)=1\-\\sum\_\{j=1\}^\{p\}\\phi\_\{j\}B^\{j\},\\quad\\theta\(B\)=1\+\\sum\_\{j=1\}^\{q\}\\theta\_\{j\}B^\{j\}\.
This allows us to write out the full definition of an ARIMA\(p,d,q\)\(p,d,q\)model as:

ϕ⁡\(B\)​∇dyt=θ⁡\(B\)​ϵt\.\\phi\(B\)\\nabla^\{d\}y\_\{t\}=\\theta\(B\)\\epsilon\_\{t\}\.\(6\)
To model seasonality, we assume the time series has a certain seasonal dependence at a fixed periodss\. For example, we might say that traffic data has a daily seasonality ofs=24s=24hours\. In those cases a SARIMA model can be used, which extends the ARIMA model with additional seasonal autoregressive and moving average terms\.

###### Definition A\.2\(Seasonal backshift and differencing operators\)\.

In the case of seasonality, the seasonal backshift operatorBsnB\_\{s\}^\{n\}is defined asBsn​yt=yt−n​sB\_\{s\}^\{n\}y\_\{t\}=y\_\{t\-ns\}, and the seasonal differencing operator∇sD\\nabla\_\{s\}^\{D\}is defined as\(1−Bs\)D\\left\(1\-B\_\{s\}\\right\)^\{D\}\.

If we letPPbe the seasonal autoregressive order,DDbe the seasonal differencing order, andQQbe the seasonal moving average order, and define their coefficients as before, we can define the seasonal AR and MA polynomials as

Φ⁡\(Bs\)=1−∑j=1PΦj​Bsj,Θ⁡\(Bs\)=1\+∑j=1QΘj​Bsj\.\\displaystyle\\Phi\(B\_\{s\}\)=1\-\\sum\_\{j=1\}^\{P\}\\Phi\_\{j\}B^\{j\}\_\{s\},\\quad\\Theta\(B\_\{s\}\)=1\+\\sum\_\{j=1\}^\{Q\}\\Theta\_\{j\}B^\{j\}\_\{s\}\.
###### Definition A\.3\.

The SARIMA\(p,d,q\)​\(P,D,Q\)s\(p,d,q\)\(P,D,Q\)\_\{s\}model is defined as the model that satisfies the following equation:

Φ⁡\(Bs\)​ϕ​\(B\)​∇d​∇sDyt=Θ⁡\(Bs\)​θ​\(B\)​ϵt\.\\Phi\(B\_\{s\}\)\\phi\(B\)\\nabla^\{d\}\\nabla\_\{s\}^\{D\}y\_\{t\}=\\Theta\(B\_\{s\}\)\\theta\(B\)\\epsilon\_\{t\}\.\(7\)

#### Implementation Details\.

By using an efficient implementation of AutoARIMA\[[18](https://arxiv.org/html/2608.20980#bib.bib7)\]in Python packages StatsForecast\[[14](https://arxiv.org/html/2608.20980#bib.bib5)\]and CuML\[[30](https://arxiv.org/html/2608.20980#bib.bib8)\], which selects the optimal order of the SARIMA based on a criterion that is minimized, we were able to scale the fitting of these even to large datasets of thousands of data points and hundreds of nodes\.

To evaluate SARIMA fairly on benchmarks which have a set default input window size, we restrict the maximum order size accordingly, such that for a benchmark with input window of sizeHH, we at most used theHH\-lagged observationyt−Hy\_\{t\-H\}\.

The maximal lag of the model can be found by collecting the maximal order of backshift operators in the differencing operators∇\\nablaand the AR and MA polynomials, which follows from the additive property of the backshift operator:Φ⁡\(Bs\)\\Phi\(B^\{s\}\)contributes a backshift operator of at most orders​PsP,ϕ⁡\(B\)\\phi\(B\)of at most orderpp,∇d\\nabla^\{d\}of at most orderdd, and∇sD\\nabla\_\{s\}^\{D\}of at most orders​DsD\. Hence, the largest lag in the model iss⁡\(P\+D\)\+p\+ds\(P\+D\)\+p\+d\.

Therefore by restricting the search space of SARIMA orders such thats⁡\(P\+D\)\+p\+d≤Hs\(P\+D\)\+p\+d\\leq H, there is a fair equivalence to the windows provided by the loaders\. We also train ARIMAs in the unrestricted regime, where we put no restrictions on their maximum size\.

For the small datasets—Chickenpox, PedalMe, and WikiMaths—model selection is done by the StatsForecast implementation of AutoARIMA\. In the restricted regime, a grid search is performed over permissible orders, with the optimization criterion being the Akaike Information Criterion \(AIC\)\. In the unrestricted regime, AutoARIMA was freely able to choose any reasonable order\. This made a practical difference only with WikiMaths, where this allowed for the fitting of a SARIMA \(seasonal ARIMA\) with period 7\. For Chickenpox, in the unrestricted regime we also encoded long period yearly seasonal features using Fourier Features withK=6K=6terms\. This is a standard trick in time series forecasting \(see[17](https://arxiv.org/html/2608.20980#bib.bib2)\) to encode seasonalities that would be too large for a SARIMA to deal with\.

Over the larger traffic datasets, for the sake of efficiency, we fit SARIMAs with a fixed order of\(4,1,4\)×\(4,1,4,12\)\(4,1,4\)\\times\(4,1,4,12\), which as the optimal fixed order we found over a grid search with limited scope\. We only evaluated on an "unfair" protocol with no restrction on ARIMA orders, to match the setup of DLinear\. The seasonality was used to allow the ARIMA to capture the slower decay of the traffic processes, while avoiding having to fit a parameter for every lag\.

Table 4:Traffic forecasting performance at the60\-minute horizon\(horizon 12\) on METR\-LA \(207 sensors\) and PEMS\-BAY \(325 sensors\)\. Metrics: MAE / RMSE / MAPE \(%\) of test split \(7:1:2\)\.Bold= best per dataset;underline= second best within the unified\-pipeline block\.

## Appendix BRe\-Evaluating Chickenpox

#### Setup\.

Based on the evaluation of[31](https://arxiv.org/html/2608.20980#bib.bib25), we train each model A3T\-GCN\[[1](https://arxiv.org/html/2608.20980#bib.bib45)\], AGCRN\[[2](https://arxiv.org/html/2608.20980#bib.bib17)\], DCRNN\[[21](https://arxiv.org/html/2608.20980#bib.bib9)\], DyGrAE\[[39](https://arxiv.org/html/2608.20980#bib.bib44)\], EvolveGCN\-H\[[28](https://arxiv.org/html/2608.20980#bib.bib43)\], EvolveGCN\-O\[[28](https://arxiv.org/html/2608.20980#bib.bib43)\], GC\-LSTM, GConvGRU\[[32](https://arxiv.org/html/2608.20980#bib.bib24)\], GConvLSTM\[[32](https://arxiv.org/html/2608.20980#bib.bib24)\], MPNN\-LSTM\[[27](https://arxiv.org/html/2608.20980#bib.bib46)\], TDE\-GNN\[[11](https://arxiv.org/html/2608.20980#bib.bib27)\], T\-GCN\[[45](https://arxiv.org/html/2608.20980#bib.bib47)\]on the officially available data in Pytorch Geometric Temporalbeing first\-order differenced, denoted asloader, and on the initially released data, denoted asraw\. The models are initialized over ten random seeds and the data split follows the standard temporal 90/10 train/test split\. We used Python 3\.12, Nvidia\-A40 46GB, torch\-geometric\-temporal 0\.56\.2\. Results are reported in Table[5](https://arxiv.org/html/2608.20980#A2.T5)\.

Table 5:Chickenpox Train/test MSE±\\,\\pm\\,std overn=10n\{=\}10seeds and generalization gap/tr = \(test−\-train\) / train MSE in both native and differenced space units\.*Native units*differ across modes \(raw = z\-scored level; loader = z\-scored first difference\): MSE is comparablewithina mode but notacrossmodes\.*Diff\-space*\(DS\) columns project raw\-mode predictions\.

## Appendix CSynthetic Graph Experiment

The node\-level equation for our synthetic process is:

Xi,t=ai​Xi,t−1\+bi​λ​∑j≠iwi,j​Xj,t−1\+ϵi,t−1\.X\_\{i,t\}=a\_\{i\}X\_\{i,t\-1\}\+b\_\{i\}\\lambda\\sum\_\{j\\neq i\}w\_\{i,j\}X\_\{j,t\-1\}\+\\epsilon\_\{i,t\-1\}\.\(8\)Let

Xt\\displaystyle X\_\{t\}=\(X1,tX2,tXN,t\),\\displaystyle=\\begin\{pmatrix\}X\_\{1,t\}\\\\ X\_\{2,t\}\\\\ \\vdots\\\\ X\_\{N,t\}\\end\{pmatrix\},ϵt−1\\displaystyle\\epsilon\_\{t\-1\}=\(ϵ1,t−1ϵ2,t−1ϵN,t−1\)\.\\displaystyle=\\begin\{pmatrix\}\\epsilon\_\{1,t\-1\}\\\\ \\epsilon\_\{2,t\-1\}\\\\ \\vdots\\\\ \\epsilon\_\{N,t\-1\}\\end\{pmatrix\}\.\(9\)Define

TA=diag⁡\(a1,…,aN\),TB=diag⁡\(b1,…,bN\),T\_\{A\}=\\operatorname\{diag\}\(a\_\{1\},\\ldots,a\_\{N\}\),\\qquad T\_\{B\}=\\operatorname\{diag\}\(b\_\{1\},\\ldots,b\_\{N\}\),and𝐀\\mathbf\{A\}the adjacency matrix, whereAi,i=0​∀iA\_\{i,i\}=0\\;\\;\\forall i, and∑jAi,j=1\\sum\_\{j\}A\_\{i,j\}=1\.

Then we can define the transition matrixT=TA\+λ​TB​𝐀T=T\_\{A\}\+\\lambda T\_\{B\}\\mathbf\{A\}, and write our synthetic process matrix form as:

Xt=T​Xt−1\+ϵt−1X\_\{t\}=TX\_\{t\-1\}\+\\epsilon\_\{t\-1\}\(10\)We assume thatϵt−1∼𝒩⁡\(0,σ2​IN\)\\epsilon\_\{t\-1\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}I\_\{N\}\)\. For weak stationarity of discrete graph process, the coefficients need to satisfy a spectral stability condition\. By noticing that the above form matches that of a VAR\(1\) process, and provided that the covariance matrix ofϵt−1\\epsilon\_\{t\-1\}is positive definite, a necessary and sufficient condition is for the spectral radius ofTTto satisfy\|ρ⁡\(T\)\|<1\|\\rho\(T\)\|<1\[[40](https://arxiv.org/html/2608.20980#bib.bib3)\]\. For our synthetic experiment, we letTA​∼iid​U​\(0\.5−δa,0\.5\+δa\)T\_\{A\}\\overset\{\\mathrm\{iid\}\}\{\\sim\}U\(0\.5\-\\delta\_\{a\},0\.5\+\\delta\_\{a\}\),TB​∼iid​U​\(0\.5−δb,0\.5\+δb\)T\_\{B\}\\overset\{\\mathrm\{iid\}\}\{\\sim\}U\(0\.5\-\\delta\_\{b\},0\.5\+\\delta\_\{b\}\)\. For each nodeii, we sample the off\-diagonal entries of𝐀\\mathbf\{A\}as𝐚i,−i​∼iid​Dirichlet​\(α​𝟏N−1\)\\mathbf\{a\}\_\{i,\-i\}\\overset\{\\mathrm\{iid\}\}\{\\sim\}\\text\{Dirichlet\}\(\\alpha\\mathbf\{1\}\_\{N\-1\}\), withAi,i=0A\_\{i,i\}=0\. The vector𝐚i,−j\\mathbf\{a\}\_\{i,\-j\}represents theii\-th row, with thejj\-th column excluded\. Theα\\alphaparameter describes the concentration of draws from the Dirichlet distribution\.

In our experimental set\-up, we simulate 500 time steps with 10 nodes\. We vary the Dirichlet concentration parameterα\\alphain the set\{0\.3,1\.0,2\.0,5\.0\}\\\{0\.3,1\.0,2\.0,5\.0\\\}and the spatial global multiplierλ∈\{0\.025,0\.05,0\.10,0\.15,0\.2\}\\lambda\\in\\\{0\.025,0\.05,0\.10,0\.15,0\.2\\\}\. For each\(α,λ\)\(\\alpha,\\lambda\)combination, we use combinations ofδa\\delta\_\{a\}andδb\\delta\_\{b\}chosen from 10 equally spaced values in the interval\[0\.015,0\.15\]\[\\sqrt\{0\.015\},\\sqrt\{0\.15\}\]\. This procedure grants a variety of node\-level temporal and spatial heterogeneity, providing a total of 2000 samples\.

Table 6:Linear regression of the GRUGCN–ARIMA test MSE gap with respect to measures of heterogeneity\.![Refer to caption](https://arxiv.org/html/2608.20980v1/figures/heterophily_heatmap.png)Figure 4:2D heatmap of the response of the OLS predicted test gapARIMAMSE−GNNMSE\\text\{ARIMA\}\_\{\\text\{MSE\}\}\-\\text\{GNN\}\_\{\\text\{MSE\}\}\. The y\-axis is the global spatial sensitivityλ\\lambda, the x\-axis is the variance of temporal coefficients across all nodes in the graph\. These were the 2 most significant coefficients from the linear regression, which had anR2R^\{2\}of 0\.591\.
## Appendix DExperimental settings

The large scale experiments performed on traffic datasets METRLA and PEMSBAY were ran with a 25 CPU cluster to distribute ARIMA grid searches\. These took under 30 minutes\. The GPU used was a Nvidia\-A40 46GB\. The small scale experiments were ran on a home laptop, with a GTX 1060 GPU\. In addition, the primary version numbers required were Python 3\.12 and torch\-geometric\-temporal 0\.56\.2\.