arXiv CSopen access

GraphToolbox: A Configurable Python Framework for Graph Neural Network Forecasting

1

GraphToolbox: A Configurable Python Framework for Graph Neural Network Forecasting

arXiv:2609.24609v1 [cs.LG] 21 Sep 2026

Eloi Campagne1,2 , Yvenn Amara-Ouali2,3 , Yannig Goude2,3 , Argyris Kalogeratos1 1 Centre Borelli, ENS Paris-Saclay, Gif-sur-Yvette, France 2 EDF Lab, Palaiseau, France 3 Laboratoire de Mathématiques d’Orsay, Université Paris-Saclay, Orsay, France

Abstract—Electricity forecasting often involves spatially related signals observed over regions, substations, and feeders, and Graph Neural Networks (GNNs) provide a natural way to represent these relations. Building a complete GNN forecasting experiment is nonetheless laborious, because graph construction, model selection, training, aggregation, and interpretation sit in incompatible tools. We present GraphToolbox, an open-source Python framework that unifies these stages in one configurationdriven pipeline built on PyTorch Geometric. It offers data-driven graph construction, an adapter that instantiates and trains 51 of the 65 PyTorch Geometric convolutions together with the recurrent cells of PyTorch Geometric Temporal, online expert aggregation, forecasting interpretability, and significance testing on cached forecasts. We evaluate the pipeline in two case studies. On French regional load, the 48 convolutions included in the complete forecasting sweep fall in a band from 1.14% to 1.60% error, online aggregation lowers this to 0.98%, and the graph models improve on classical additive and boosting baselines. On net-load, direct graph models are less accurate than a classical additive model, while forecasting each physical component separately improves them without closing that gap. Both comparisons use the same experimental interface, illustrating the role of GraphToolbox in systematic architectural evaluation. Index Terms—graph neural networks, electricity load forecasting, open-source software, spatio-temporal modeling, expert aggregation

I. I NTRODUCTION The operation of modern power systems rests on accurate short-term forecasts of electricity demand, which inform market decisions, unit commitment, and grid balancing. Historically, this task was addressed at the level of an aggregated national signal, where generalized additive models and, more recently, boosting and deep learning methods have proven remarkably effective [1]. Aggregation describes the measurement scale and does not remove the underlying spatial organization of the system. Electricity demand is distributed across interconnected regions, substations, and feeders whose components co-vary in ways that a single aggregated series cannot express. The decentralization of production, the integration of intermittent renewables, and the deployment of smart metering make this structure increasingly observable [2]. Graph Neural Networks offer a natural inductive bias for this setting. By restricting information exchange to graph neighborhoods, they encode the assumption that spatial or statistical proximity governs load co-variation, an assumption Corresponding author: [email protected].

that has driven their success in traffic forecasting and, increasingly, in power systems [3], [4]. A rigorous empirical comparison nevertheless requires several methodological and software choices. A practitioner must first decide how to turn a collection of time series into a graph, then choose among a rapidly growing family of convolution operators, integrate those operators into a training loop with appropriate metrics, combine several models to gain robustness, and finally interpret the resulting predictions. Each of these steps is supported by a different library, or by no library at all, and the integration code is often specific to a single project. GraphToolbox provides a common interface for these stages. The framework is an open-source Python package that organizes the entire forecasting workflow behind a small number of configuration dictionaries. Changing the graph, the convolution, or the aggregation strategy then requires modifying a configuration field while keeping the rest of the pipeline fixed. Its contributions are the following. First, it provides a menu of data-driven graph-construction methods that derive an adjacency structure from spatial coordinates or from the load signals themselves, together with an empirical selection procedure. Second, it introduces a convolution adapter that lifts 51 of the 65 PyTorch Geometric operators, and through a parallel adapter the recurrent cells of PyTorch Geometric Temporal, into a common forecasting model, which allows broad architectural comparisons within the same implementation. Third, it integrates online expert aggregation, which combines the forecasts of several models with weights that adapt to past performance. Fourth, it supplies interpretability tools tailored to forecasting, namely accumulated local effects for feature attribution and explanation graphs that display the spatial connections associated with a prediction. Fifth, it provides significance testing on cached forecasts, from pairwise Diebold–Mariano tests to the Model Confidence Set, which associates architectural comparisons with measures of statistical uncertainty. We describe the design of the framework, report the coverage of its convolution benchmark, and demonstrate its use on regional load and net-load forecasting tasks drawn from European electricity data. II. R ELATED TOOLS PyTorch Geometric supplies convolution operators and the sparse primitives of message passing [5], and PyTorch Geometric Temporal adds recurrent spatio-temporal cells and

2

TABLE I D ESIGN EMPHASIS OF G RAPH T OOLBOX RELATIVE TO P Y T ORCH G EOMETRIC T EMPORAL (P Y G-T) AND T ORCH S PATIOTEMPORAL ( TSL ), SCORED AGAINST THE RELEASED MODULES OF EACH LIBRARY AS OF S EPTEMBER 2026. T HE COMPARISON REFLECTS PRIMARY DESIGN FOCUS RATHER THAN AN EXHAUSTIVE AUDIT.

Capability

GraphToolbox

PyG-T

tsl

Shared foundation Built on PyTorch Geometric

✓

✓

✓

GraphToolbox emphasis Adapter over arbitrary convs Data-driven graph construction Empirical graph selection Online expert aggregation Forecasting interpretability Forecast significance testing Rolling-origin evaluation Integrated hyperparameter search

✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓

∼ – – – – – ∼ –

∼ ∼ – – – – ✓ ∼

Complementary strengths Prebuilt benchmark datasets Curated ST model zoo Missing-data imputation

∼† ∼† –†

✓ ✓ –

✓ ✓ ✓

✓: built-in capability, ∼: partial capability, –: no relevant capability, †: capability planned for the next GraphToolbox release.

models built from them [6]. GraphToolbox sits one level above and reimplements neither: an operator or recurrent cell is named in a configuration dictionary, lifted into a single metamodel, and carried from graph construction to prediction. The lower-level libraries leave open how these blocks assemble into a forecasting experiment, and that assembly is what GraphToolbox fixes. Closest in spirit is Torch Spatiotemporal, which provides data processing and model prototyping for neural forecasting over sensor networks [7]. It and PyTorch Geometric Temporal remain general, supplying layers and data structures for arbitrary spatio-temporal signals, whereas GraphToolbox makes graph construction, architectural benchmarking, aggregation and interpretability explicit stages of one workflow. Toolkits built around global or foundation models forecast collections of series without a relational structure [8]. Table I compares these stages. We score a capability as built-in, partial, or absent by whether a released module provides it directly. The adapter row requires instantiating and training an arbitrary PyTorch Geometric operator without operator-specific code, which the compared libraries offer only for a curated subset of layers. Data-driven graph construction requires deriving the adjacency from the series, so a library that accepts only a supplied or coordinate graph scores absent. Rolling-origin evaluation requires an expanding-window trainer exposed as a module. Forecasting interpretability is scored on the accumulated-localeffects module, the intrinsically additive graph model and the rendering of edge weights on the map, the fit of a post-hoc explainer being left to the user. The final block is scored the other way, collecting what PyG-T and tsl expose as firstclass modules while GraphToolbox only wraps or omits them, namely dataset zoos, curated spatio-temporal models and imputation, planned for the next release rather than claimed here.

III. F RAMEWORK DESIGN GraphToolbox is organized into modules that mirror the stages of a forecasting experiment (Figure 1). The graph, the architecture and the training regime are configuration dictionaries consumed by fixed execution code, so one factor varies at a time without bespoke scripts; the whole load study of Section V is driven by the four dictionaries of Figure 1. A. Data handling and graph construction The DataClass reads training and test tables, applies the temporal split, and manages calendar encodings and lagged features. The GraphDataset turns the tables into graph snapshots and fits feature and target scalers on the training window only; its field out_channels sets the forecasting horizon. The GraphBuilder implements the adjacency constructions proposed for load forecasting [9], each returning positive edge weights that grow with the strength of the relation between two nodes. A spatial graph thresholds a Gaussian kernel of geodesic distance. Data-driven graphs use dynamic time warping [10], correlation and precision matrices, a spline distance between temperature-to-load responses, or GL3SR, which learns a Laplacian from graph signals under a smoothness criterion [11]. Since no construction dominates across datasets, each candidate graph can be used to train the same model, and the graph minimizing validation error is retained. B. Models and the convolution adapter Every static architecture is expressed through one metamodel, myGNN: a linear encoder lifts node features to a fixed latent width, residual message-passing blocks in the pre-activation DeepGCN arrangement [12] refine it, and a linear read-out returns the multi-step forecast at every node in one pass. Each block couples a convolution with layer normalization and a rectified linear activation, and multi-head convolutions split the latent width across heads so that the head count inflates neither capacity nor computation. The convolution itself is supplied by name. myGNN receives a convolution class together with its keyword arguments and instantiates it through a ConvAdapter, the component that makes a catalog-wide sweep possible. The operators of PyTorch Geometric obey no common calling convention, some expecting scalar edge weights, others vector edge attributes, a relation type, or a positional input, and many expecting none of these. The adapter reads the constructor signature of the operator it is given, forwards only the arguments that operator accepts, and fills in the defaults it requires, whether a Chebyshev order, a relation count for the relational convolutions, or the small perceptrons that the isomorphism and edge convolutions expect. It then reads the forward signature, routes to the operator only the graph tensors it consumes, and supplies a neutral placeholder for any input the operator demands yet the forecasting graph does not carry. When an operator returns a width different from the latent one, as the multi-hop and signed convolutions do, a linear projection restores it, and the few operators with CPU-only kernels are pinned there with their tensors ferried across the device boundary transparently.

3

aggregation DataClass

GraphBuilder

GraphDataset

CSV / splits

graph construction

scaling, lags

myGNN ConvAdapter ConvAdapterTemporal

Trainer

metrics

RollingTrainer

significance data_kwargs

dataset_kwargs

conv_class, conv_kwargs

model_kwargs,

Optuna spaces

interpretability

Fig. 1. Overview of the GraphToolbox pipeline. Solid arrows carry data from raw series to forecasts and diagnostics, while dashed boxes represent the configuration dictionaries that parameterize each stage. Swapping the graph, the convolution operator, or the training regime amounts to editing a configuration field. from torch_geometric.nn.conv import GATConv from graphtoolbox.models import myGNN from graphtoolbox.training import Trainer model = myGNN( in_channels=ds_train.num_node_features, num_layers=2, hidden_channels=256, out_channels=48, # forecast horizon conv_class=GATConv, # <- swap operator here conv_kwargs={"heads": 2}) trainer = Trainer(model, ds_train, ds_val, ds_test, batch_size=16, model_kwargs={"lr": 1e-3}) pred, target, edge_index, attn = trainer.train()

Listing 1. Swapping a convolution is a one-field change. The same Trainer and GraphDataset serve any operator exposed through ConvAdapter.

The same path can expose the coefficients of the attention operators, which the interpretability module later renders as edge weights. Swapping a graph convolutional network for an attention network or a Chebyshev filter is then the one-field change of Listing 1. Operator-specific cases are implemented once and maintained centrally. Argument routing occupies about 240 lines of adapter code, and the significance module a further 340. The supported families include degree-normalized spectral convolutions [13], inductive neighborhood aggregators [14], personalized propagation through APPNP [15], attention operators from GATConv to GATv2Conv and TransformerConv [16]–[18], and polynomial spectral filters [19]. A parallel adapter, ConvAdapterTemporal, carries the same idea to the recurrent graph cells of PyTorch Geometric Temporal [6]. It reads the constructor signature to distinguish a steprecurrent cell, which consumes one time step at a time while carrying a graph-aware hidden state, from a windowed cell, which ingests the whole window and must be told the number of periods in advance. Its companion model, TemporalGNN, plays the part of myGNN for these cells. It separates the ordered lag sequence of the target from the static side features, unrolls that sequence through the cell into a hidden state, concatenates the static features, and projects the result to the forecast horizon. A numerically verified fast unroll of the Chebyshev gated cell folds the input-side gate transforms over the whole sequence and computes the graph normalization once, reverting to the reference step loop whenever an equivalence check fails. Since TemporalGNN exposes the attributes the Trainer reads, a recurrent cell trains, reconciles, and caches exactly as a static convolution does, which is what allows the case study of Section V to place the two families under one budget.

C. Training, optimization, and aggregation Training is handled by the Trainer, which manages the optimization loop, early stopping, and the evaluation of the standard forecasting metrics, namely mean absolute error, its normalized variant, mean absolute percentage error, root mean squared error, and bias. Optimization relies on Adam [20]. A RollingTrainer extends this loop to a rolling-origin protocol, in which the training window expands to absorb each past test window before predicting the next, which mirrors the way a forecasting system is retrained in operation. Hyperparameter search is delegated to an Optimizer built on Optuna, so that the number of layers, the hidden width, the learning rate, the batch size, and even the graph become tunable through a declarative search space [21]. Because different architectures carry different inductive biases, combining them is often more robust than trusting any single one. The framework provides a native Aggregation module that performs sequential expert aggregation either online or in batch, following the aggregation rules of the opera package [22], whose R and Python1 reference implementations informed our design. Its default rule is the parameterfree MLpol algorithm, which maintains a weight per expert and updates it according to past forecasting performance, so that the aggregated prediction is a time-varying convex combination of the experts, with exponentially weighted and Bernstein online variants also available. A uniform average is provided as a strong and parameter-free baseline, and an ensemble of independently initialized models supplies a first estimate of predictive spread. Aggregation can be applied either across architectures or across the initializations of a single architecture, and it can be arranged top-down at the aggregate level or bottom-up at the node level. D. Interpretability and visualization The interpretability module provides feature- and edge-level diagnostics for individual forecasts. Accumulated local effects estimate the marginal influence of a feature while accounting for its correlation with other inputs, which makes them more reliable than partial dependence in the highly collinear setting of weather and calendar drivers [23]. For the spatial structure, edge attributions can be rendered on the graph, whether they come from an explainer of PyTorch Geometric, which the user runs, or from attention weights when the architecture provides them [24]. These diagnostics are correlational, and we make 1 https://github.com/Dralliag/opera-python

4

no claim that they establish a causal link between a connection and a prediction. They are also no more stable than the attribution they are given, and Section V shows that stability is not to be assumed. Beyond these post-hoc attributions, the package also includes an intrinsically interpretable model, an additive graph network that encodes each feature group with its own subnetwork and sums their contributions, so that a prediction decomposes into per-group terms, each of which can be read on its own. A visualization module places node errors, learned graphs, and edge weights on maps of the corresponding territory. E. Significance testing An accuracy ranking without a measure of its reliability leaves open whether a gap reflects a genuine difference or the noise of a finite test window. The evaluation module operates entirely on cached forecasts, allowing the statistical analysis to run on stored predictions without retraining the models. Pairwise predictive accuracy is assessed with the Diebold– Mariano test [25] under a squared, absolute, or percentage loss, with a Newey-West long-run variance [26] and the smallsample correction of Harvey, Leybourne, and Newbold [27]. A Holm-Bonferroni adjustment controls the family-wise error rate across the pairwise comparisons. Sampling uncertainty on a metric is quantified by a moving-block bootstrap [28], whose blocks of consecutive steps preserve the autocorrelation of forecast errors that an ordinary resampling would destroy, with a default block of one day and 2000 resamples. To summarize an entire sweep, the module reports the Model Confidence Set [29], the subset of architectures that contains the best one with a prescribed probability, which gives a precise meaning to a group of operators the data cannot separate. These instruments address the sampling uncertainty of a metric on a finite test window, and they leave a second source untouched, the stochasticity of training itself. Because the sweeps of Section V report one run per architecture, we retrained the strongest operators over 5 seeds at their table configurations. Table II sets the resulting means against the entries of the main tables, which sit below them by 58 to 147 MW, so those entries are best read as one draw rather than as an expectation. Two conclusions remain stable across seeds. At the seed mean, the load model remains 148 MW below the recurrent baseline. For net-load, the paired improvement from direct to decomposed forecasting averages 200 MW with a standard deviation of 58 MW. The fine ordering within each architectural band is not stable, which is why we claim no unique best convolution. F. Scale On the 12-node French graph, 9 of the 10 representative convolutions train in 130 to 300 ms per epoch, so a full convergence run takes about a minute; only the relationalattention RGATConv is an order of magnitude slower. These timings were measured on the CPU of a MacBook Pro with an Apple M4 Pro (12 cores), which outperforms the MPS backend at this graph size.

TABLE II T RAINING - SEED VARIABILITY OF THE LEADING MODEL OF EACH STUDY. NATIONAL RMSE (MW), MEAN ± STANDARD DEVIATION OVER 5 SEEDS AT THE CONFIGURATION OF THE MAIN TABLES , AGAINST THE SINGLE RUN THOSE TABLES REPORT. Task and model

5 seeds

table entry

Load, LEConv Net-load direct, GCN2Conv Net-load decomposed, GatedGraphConv

930 ± 43 2185 ± 65 1985 ± 43

834 2038 1927

TABLE III C OST OF THE SAME CODE PATH AT TWO GRAPH SCALES , ON THE 12 F RENCH REGIONS AND ON 500 WEAVE-UK FEEDERS , MEASURED ON THE A PPLE M4 P RO CPU AT HIDDEN 256 AND 2 LAYERS . PARAMETER COUNTS DIFFER BETWEEN THE TWO COLUMNS BECAUSE THE INPUT DIMENSION DOES (187 FEATURES AGAINST 50). WALL - CLOCK TRAINING TIME IS OMITTED , THE F RENCH MODELS BEING READ FROM A CHECKPOINT WHILE THE WEAVE ONES ARE FITTED HERE , WHICH MAKES THE TWO INCOMPARABLE .

12 nodes

500 nodes

Operator

params (k)

infer (ms)

params (k)

infer (ms)

GraphSAGE GAT APPNP

32.0 24.4 62.0

0.57 3.81 3.04

23.3 15.6 26.9

6.87 31.4 360.2

Table III carries the more informative comparison. It repeats the measurement on a graph of 500 low-voltage feeders, the WEAVE-UK collection, through the same code path and on the same machine. Inference per forecast window rises from well under a millisecond to a few milliseconds for the neighborhood aggregator and from 4 to 31 ms for the attention operator, so a 40-fold increase in node count corresponds to roughly one order of magnitude higher latency for these two operators. APPNP behaves differently because it propagates over the dense correlation graph used for WEAVE in the absence of coordinates. Its latency reaches 360 ms, indicating the importance of edge density in addition to node count. These measurements show that sweeps on graphs of this size remain feasible on a workstation. They do not characterize computational cost on graphs with thousands of nodes. IV. C ONVOLUTION COVERAGE In practice, the breadth of a forecasting framework is measured by the fraction of available operators it can actually run. We instantiated every convolution in torch_geometric.nn.conv inside myGNN and attempted an end-to-end training pass on a homogeneous graph with standard node features. Table IV reports the outcome by family. 51 of the 65 tested operators, that is 78.5%, run without modification. A single operator, the fused attention kernel, requires an optional dependency that is unavailable on some platforms. The remaining 13 presuppose heterogeneous graphs, pointcloud inputs, or device-specific libraries, and therefore fall outside the homogeneous node-level regime that electricity forecasting occupies. This benchmark measures whether an operator instantiates and completes a training pass, and does not by itself certify forecasting accuracy. Of the 51 operators supported by the adapter, 48 enter the load-forecasting sweep

5

TABLE IV C OVERAGE OF P Y T ORCH G EOMETRIC CONVOLUTION OPERATORS WHEN RUN END - TO - END INSIDE myGNN . S KIPPED OPERATORS REQUIRE HETEROGENEOUS GRAPHS , POINT- CLOUD INPUTS , OR DEVICE - SPECIFIC LIBRARIES . Operator family

Working

Examples

GCN / spectral Attention-based MPNN / aggregation GIN-style (MLP) Edge-conditioned Recurrent / gated Residual / deep Spectral / poly Dynamic aggregators Relational Graph-level

8 6 8 2 6 3 5 5 3 3 2

GCN, Cheb, SG, GCN2, FA GAT, GATv2, Transformer SAGE, GEN, Graph, LE GIN, GINE NN, CG, GMM, General GatedGraph, ARMA, TAG FiLM, ResGated, PDN MixHop, FeaSt, PAN PNA, Edge, DynamicEdge RGCN, RGAT, FastRGCN WL, Signed

Working total Optional dependency Skipped (out of scope) Tested total

51 1 13 65

(78.5%) FusedGAT hetero, point-cloud, CUDA

and its significance analysis. Out of the supported operators, 4 are excluded because their required inputs do not match this forecasting setting, while the personalized-propagation layer is added. The results report a representative selection, including the 10 static convolutions of Table V and 4 temporal cells. V. C ASE STUDIES A. Regional load forecasting The primary case study forecasts French electricity load at the granularity of the 12 administrative regions supplied by the transmission operator, using half-hourly RTE open data [30] enriched with calendar features and Météo-France SYNOP weather observations [31]. The nodes of the graph are the regions, each positioned at the spatial coordinates of a major city of that region. A model reads its features at midnight and returns the 48 half-hours of the day that follows, so every input it consumes, the 48 half-hourly load lags included, predates that origin. The sweep holds the graph at the dynamic-time-warping construction and the architecture at 2 layers of width 64, which makes the comparison in Table V a direct product of the configuration mechanism rather than of separate implementations. Three operators in the full sweep, LEConv, GATConv, and APPNP, use configurations returned by a model-and-graph hyperparameter search; Table V includes the first two, whose entries carry that advantage. Table V collects the benchmark on the 2019 test set. The 48 convolutions in the sweep range from 1.14% to 1.60% MAPE, with a median of 1.31%; LEConv is best at 1.14% and 834 MW. The nonlinear temperature features built during preprocessing account for the margin over the numbers of a prior study on the same data [9]. Every operator of Table V improves on every reference: the national and regional additive models reach 1200 and 1248 MW, gradient boosting 1416 MW, a covariateconditioned recurrent network 1078 MW, and Chronos-2 1647 MW. Chronos-2 and Chronos-Bolt belong to the Chronos family of time-series foundation models [8]. The recurrent

network is the strongest reference; its gap to the graph models ranges from 143 MW at the weakest listed operator to 244 MW at the best. That gap is not by itself evidence for the spatial inductive bias, since a graph model also differs from a recurrent baseline in its per-node parameterization and its multi-step read-out. Isolating message passing requires holding everything else fixed. We therefore retrain the same 10 operators in dedicated paired runs at a strictly common configuration on the DTW and identity graphs; these runs are not directly comparable to the entries of Table V. Averaging their individual RMSEs gives 951 MW on the DTW graph and 984 MW on the identity graph. The median paired reduction is 23 MW, and 8 of the 10 operators improve. The effect is architecture-dependent, with the identity graph improving 2 operators. Message passing therefore contributes to, but does not account for, the overall gain. Online expert aggregation over the whole sweep tightens the result, from 1.11% for a uniform average to 0.98% for the MLpol mixture, at a national root mean squared error of 750 MW, 84 MW below the best single convolution. The block-bootstrap intervals of Table V do not resolve the ordering inside the sweep, and neither does a statistical test: a Model Confidence Set at the 90% level, computed on the same cached forecasts, retains 36 of the 48 operators, the 10 among those in the table, and a Holm-corrected Diebold– Mariano test separates none of those 10 from the best. The mixture is nevertheless significantly more accurate than the best single convolution under a paired Diebold–Mariano test (p = 0.001), despite the overlap of their marginal bootstrap intervals. Aggregation therefore removes the need to identify the best operator in advance, although it still requires running the whole sweep. We also evaluate MinT reconciliation and exclude it from the reported forecasts, because it degrades them. MinT projects the 12 regional forecasts onto a national forecast, here produced by the gradient-boosting model the framework fits on the aggregate. That model reads the load of one day and one week before the forecast origin, two exponentially weighted averages frozen at that origin, and calendar encodings of the half-hour, the weekday and the month, so its information set respects the same day-ahead constraint as the graph models. That national forecast reaches an error of 2385 MW over the test year, well above the regional models it would correct, and the projection accordingly degrades them: LEConv moves from 834 to 943 MW and the MLpol mixture from 750 to 913, and the mixture no longer separates from the best single operator (p = 0.10). Table V therefore reports the unreconciled forecasts, and the reconciliation is worth its cost only where the top level is the better predictor. The harness also runs T-GCN, A3T-GCN, DCRNN, and GConvGRU [32]–[35], either one-shot, updating once on the stacked lag channels as a static convolution does, or unrolled over the ordered lags. Run one-shot, GConvGRU and DCRNN reach 857 and 878 MW, respectively, within the range of the static convolutions but above the 750 MW aggregate, so they do not improve on the best results of the static sweep.

6

TABLE V F RENCH NATIONAL LOAD , 2019 TEST SET. MAPE AND RMSE (MW), SWEEP ORDERED BY RMSE, FORECASTS UNRECONCILED . I NTERVALS ON THE GNN ROWS ARE BLOCK - BOOTSTRAP STANDARD ERRORS ON A SINGLE TRAINING RUN ( CF. TABLE II), THE BASELINES BEING POINT ESTIMATES . T HE REPRODUCIBILITY SECTION AT THE END OF THE PAPER STATES WHICH ROWS ARE BORROWED AND WHICH ARE TRAINED HERE . Model Persistence (1 day) Chronos-Bolt ARIMA-X Chronos-2 XGBoost GAM (regional) GAM (national) LSTM

MAPE (%)

RMSE (MW)

5.78 2.99 2.51 1.78 1.63 1.84 1.67 1.46

4507 2408 1706 1647 1416 1248 1200 1078

Static GNN sweep (GraphToolbox) LEConv 1.14 ± 0.03 GatedGraphConv 1.22 ± 0.03 GATv2Conv 1.22 ± 0.03 ARMAConv 1.24 ± 0.03 ResGatedGraphConv 1.22 ± 0.03 SGConv 1.24 ± 0.03 RGATConv 1.22 ± 0.03 RGCNConv 1.23 ± 0.03 GCNConv 1.20 ± 0.03 GATConv 1.23 ± 0.03

834 ± 35 881 ± 29 881 ± 30 887 ± 32 890 ± 36 898 ± 33 899 ± 34 910 ± 43 925 ± 46 935 ± 47

Expert aggregation Uniform average MLpol (bottom-up) MLpol (top-level)

828 ± 38 754 ± 37 750 ± 36

1.11 ± 0.03 0.98 ± 0.03 0.98 ± 0.03

B. Net-load forecasting The second case study forecasts national net-load, the demand that remains once solar and wind generation are subtracted, which is the quantity a system operator must actually balance. This target is harder than gross load, because it changes sign and inherits the volatility of renewable production. It also exposes a subtlety that the framework handles explicitly. The percentage error is meaningless for series that pass through zero, such as solar generation at night, so the pipeline restricts that metric to strictly positive series and reports root mean squared error and a symmetric percentage error elsewhere. We write sMAPE for the symmetric mean absolute percentage error, which is bounded between 0% and 200% and remains well defined near zero. On the same 2019 test set (Table VI), the 10 convolutions, all at the same validation-selected configuration (1 layer of width 64), occupy a band from 3.26% to 3.69% sMAPE, GCN2Conv the best at 2038 MW, and online aggregation tightens the direct sweep to 1963 MW. The classical references, however, tell a different story than on gross load. Persistence, the zero-shot Chronos-Bolt at 3676 MW, and ARIMA-X at 2933 trail far behind, and Chronos-2 reaches 2344 without matching the fitted supervised models. Gradient boosting at 2031 MW already matches the best direct convolution. The national GAM, one additive model per half-hour fitted on weather and calendar covariates, is stronger than every direct convolution and their aggregate at 1763 MW, and the recurrent baseline reaches 1746 MW. Fitting the same GAM per region and summing its forecasts gives the lowest RMSE, 1723 MW, while the recurrent baseline has

the lowest sMAPE. Net-load is dominated by weather-driven wind and solar generation, which the additive models capture directly through their power-weighted meteorological terms, so the spatial inductive bias that carries gross load does little for a model that predicts net-load in one piece. We thus evaluate a natural decomposition of net-load. Since net-load = load − wind − solar, forecasting each component with its own GNN and recombining the forecasts, using the same configuration and training budget per model, improves 9 of the 10 convolutions and lowers RMSE by 162 MW averaged over all 10, GatedGraphConv reaching 2.97% and 1927 MW and the decomposed aggregate 2.87% and 1884 MW. The gain is systematic and insufficient: a paired test separates the best decomposed forecast from the best direct one (p = 0.03), and at the seed means the paired gap is 200 MW against a standard deviation of 58, yet the decomposed aggregate remains 121 MW above the national GAM, 161 MW above the regional GAM, and 138 MW above the recurrent model. Physical structure moves the graph models from the rear of this field into its middle, and no further. The decomposition also costs three models for one forecast, so its gain is not matched on total training compute. The decomposition identifies wind as the main source of residual error. Across the 10 operators, wind RMSE ranges from 1709 to 1786 MW, against 825 to 1055 MW for load, and 392 to 521 MW for solar. The squared wind RMSE accounts for roughly three quarters of the sum of the squared component RMSEs, which approximates its share of the recombined error when the component errors are independent. Its error correlates only moderately with the recombined error across operators (rank correlation 0.43). The recurrent cells behave as on load. Unrolled over the 48 lags, GConvGRU and DCRNN reach 2280 and 2289 MW, inside the band of the direct sweep of Table VI but above its median and ahead of only its weakest operator, while A3TGCN and T-GCN exceed 5% sMAPE, so the static sweep keeps the advantage here as well. C. Interpretability diagnostics The visualization module displays spatial attributions on the geographic domain, with one explanation graph per period alongside the accumulated local effects of the feature groups. Applied to the French load model, this diagnostic reveals substantial instability in the edge attributions, at monthly resolution. Figure 2 draws the 10% most important edges of the attention model in January and in July, and the two panels contain markedly different sets of leading edges. Over the 12 months of the test year, two monthly rankings share 30% of their leading edges, against 10% for random subsets of that size, and adjacent months are no more similar than distant ones. The attention weights account for the instability. Averaged over the test year they depart from the inverse indegree of the receiving node by a median of 1.4%, and their largest and smallest values differ by under 1%, so the model assigns nearly uniform weights to each node’s neighbors. Across edges, two monthly maps correlate at 0.37 and two half-year maps at 0.87: one month fixes about a third of the

7

group, normalized across groups, gives 67.6% to the net-load lags, 19.6% to the temperature lags, 11.9% to the remaining exogenous covariates and below 1% to the calendar indicators. Unlike the unstable ranking of nearly equal attention weights, this analysis exposes the additive terms of the fitted predictor directly; it remains associational and is not a causal attribution.

TABLE VI F RENCH NATIONAL NET- LOAD , 2019 TEST SET. S MAPE AND RMSE (MW) FOR DIRECT AND DECOMPOSED GNN S AND THEIR AGGREGATIONS . D IRECT AND DECOMPOSED MODELS USE THE SAME PER - MODEL CONFIGURATION AND TRAINING BUDGET. I NTERVALS ARE BLOCK - BOOTSTRAP STANDARD ERRORS FROM A SINGLE TRAINING RUN ( CF. TABLE II); BASELINES ARE POINT ESTIMATES . B EST VALUES ARE IN BOLD . Model Persistence (1 day) Chronos-Bolt ARIMA-X Chronos-2 XGBoost GAM (national) LSTM GAM (regional)

sMAPE (%)

RMSE (MW)

VI. D ISCUSSION AND LIMITATIONS

8.89 5.70 5.04 3.43 3.33 2.92 2.83 2.86

5636 3676 2933 2344 2031 1763 1746 1723

The case studies support two methodological observations. The evaluated graph models are shallow, with 2 messagepassing layers for load and 1 for net-load, which keeps their computational cost modest. This design is consistent with over-smoothing analyses that predict a loss of node-specific signal as depth grows [36]; depth should therefore be treated as a hyperparameter selected on validation data rather than assumed in advance. Online expert aggregation improves the best individual operator under the paired predictive-accuracy test, both at the national level and when applied bottom-up. Hierarchical reconciliation is beneficial only when the toplevel forecast is more accurate than the forecasts it constrains. This condition is not met for French load, and the integrated pipeline allows it to be checked without implementing a separate reconciliation workflow. A number of limitations remain. GraphToolbox is at an early release stage, and its public interface may evolve. It targets homogeneous node-level graphs and therefore excludes heterogeneous graphs and point-cloud operators. The interpretability tools inherit the known caveats of attention-based and perturbation-based attribution, and Section V illustrates this limitation on a small dense graph, where the resulting edge rankings are unstable. The significance analysis conditions on a single training run per architecture, so its intervals and tests quantify variability over the test window and not across retraining; the seed study of Section III measures the second source separately, and the two would be better quantified jointly. The package ships with continuous integration and a test suite.

Direct GNN, one model on net-load GCN2Conv 3.26 ± 0.10 FastRGCNConv 3.40 ± 0.10 RGATConv 3.41 ± 0.11 GatedGraphConv 3.43 ± 0.11 ResGatedGraphConv 3.48 ± 0.11 FiLMConv 3.52 ± 0.12 GMMConv 3.42 ± 0.11 PANConv 3.51 ± 0.11 GATConv 3.61 ± 0.13 SuperGATConv 3.69 ± 0.14

2038 ± 62 2078 ± 62 2145 ± 71 2150 ± 64 2188 ± 69 2193 ± 68 2205 ± 74 2233 ± 77 2257 ± 73 2325 ± 94

Direct GNN, expert aggregation MLpol (top-level) 3.14 ± 0.11 Uniform average 3.09 ± 0.10 MLpol (bottom-up) 3.06 ± 0.10

2000 ± 71 1980 ± 67 1963 ± 65

Decomposed GNN, one model per component GatedGraphConv (best) 2.97 ± 0.10 1927 ± 59 MLpol (top-level) 2.90 ± 0.09 1890 ± 60 MLpol (bottom-up) 2.87 ± 0.09 1884 ± 62

January

July

5

5

7

7 6

2

4

3 1

10

0.0852

6

2

4

3 1

10

0

0.0849

8 9

11

0.0851 0.0850

0 8

Edge importance (mean edge_mask)

Explanation Graph (model): GATConv

9

11

0.0848

Fig. 2. Explanation graph of the attention model on French regional load, in January and in July. Only the 10% of edges with the highest mean attention are drawn, and edge color encodes that mean. The color scale runs from 0.0848 to 0.0853, a spread of 0.6%, so the model weights a node’s neighbors almost uniformly.

ranking variation it produces, and its edge ordering should not be read. This result concerns this model on a dense 12-node graph and does not establish a general limitation of attribution methods. The maps are therefore used here as diagnostics rather than as evidence for stable spatial relations. The feature-level side of the module is demonstrated on netload with the intrinsically additive graph model. Its 2386 MW RMSE is not competitive with the forecasting models of Table VI, so the purpose of the fit is diagnostic rather than another benchmark claim. Aggregating its accumulated local effects by the mean per-feature RMS importance within each

VII. C ONCLUSION GraphToolbox unifies graph construction, architectural benchmarking, online aggregation and interpretability behind one configuration-driven interface on PyTorch Geometric. It runs 51 convolution operators unmodified. On French load, aggregation reaches 0.98% MAPE, while the identity-graph ablation isolates a smaller but systematic contribution from message passing. On net-load, the same experimental interface shows instead that direct GNNs trail classical baselines and that physical decomposition narrows, but does not close, this gap. The significance analysis identifies no unique best architecture, while hierarchical reconciliation degrades the forecasts in this case. The package is released under GPL-3.0 with documentation and examples. R EPRODUCIBILITY AND DISCLOSURE The source code, documentation, and the example notebooks underlying the case studies are publicly available at https://github.com/eloicampagne/GraphToolbox and installable

8

via pip install graphtoolbox, with a lock file fixing the interpreter and every library version, the random seeds pinned, and the hardware of the benchmarks stated. A reproducibility guide maps each result to its runner. The repository carries the script that fits every baseline of Table V and Table VI, and the cached forecasts of both sweeps and of the identity-graph ablation, so that the bootstrap, the Diebold–Mariano tests and the Model Confidence Set rerun on stored predictions; the trained weights are distributed as a release asset. Since the tables report one fixed set of checkpoints, retraining a row reproduces it only up to the seed spread of Table II. The French load and net-load data come from RTE open data [30] and Météo-France SYNOP observations [31], January 2015 to December 2019 at half-hourly resolution over the 12 regions. In Table V, the rows persistence, Chronos-Bolt, gradient boosting and the national additive model, come from a prior study on the same data [9]; everything else was produced with GraphToolbox for this paper, the baselines being trained as national-direct forecasters except the regional additive model, which is fitted per region and summed. E. Campagne carried the project end to end, from the design of the software and the experiments to the writing. Y. Amara-Ouali, Y. Goude, and A. Kalogeratos contributed through review and supervision. The authors declare no competing interests. R EFERENCES [1] M. Fasiolo, S. N. Wood, M. Zaffran, R. Nedellec, and Y. Goude, “Fast calibrated additive quantile regression,” J. Amer. Statist. Assoc., vol. 116, no. 535, pp. 1402–1412, 2021. [2] J. de Vilmarest, J. Browell, M. Fasiolo, Y. Goude, and O. Wintenberger, “Adaptive probabilistic forecasting of electricity (net-) load,” IEEE Trans. Power Syst., vol. 39, no. 2, pp. 4154–4163, 2023. [3] S. Guo, Y. Lin, N. Feng, C. Song, and H. Wan, “Attention based spatialtemporal graph convolutional networks for traffic flow forecasting,” in AAAI Conference on Artificial Intelligence, 2019. [4] E. Campagne, Y. Amara-Ouali, Y. Goude, and A. Kalogeratos, “Leveraging graph neural networks to forecast electricity consumption,” in Joint European Conference on Machine Learning and Knowledge Discovery in Databases. Springer, 2024, pp. 309–328. [5] M. Fey and J. E. Lenssen, “Fast graph representation learning with PyTorch Geometric,” in ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019. [6] B. Rozemberczki, P. Scherer, Y. He, G. Panagopoulos, A. Riedel, M. Astefanoaei, O. Kiss, F. Beres, G. López, N. Collignon, and R. Sarkar, “PyTorch Geometric Temporal: Spatiotemporal signal processing with neural machine learning models,” in ACM International Conference on Information and Knowledge Management, 2021, pp. 4564–4573. [7] A. Cini and I. Marisca, “Torch Spatiotemporal,” https://github.com/ TorchSpatiotemporal/tsl, 2022. [8] A. F. Ansari, L. Stella, C. Turkmen, X. Zhang, P. Mercado, H. Shen, O. Shchur, S. S. Rangapuram, S. Pineda Arango, S. Kapoor et al., “Chronos: Learning the language of time series,” Trans. Mach. Learn. Res., 2024. [9] E. Campagne, Y. Amara-Ouali, Y. Goude, I. Zehavi, and A. Kalogeratos, “Graph neural networks for electricity load forecasting,” Preprint arXiv:2507.03690, 2025. [10] S. Salvador and P. Chan, “FastDTW: Toward accurate dynamic time warping in linear time and space,” in KDD Workshop on Mining Temporal and Sequential Data, 2004. [11] P. Humbert, B. Le Bars, L. Oudre, A. Kalogeratos, and N. Vayatis, “Learning Laplacian matrix from graph signals with sparse spectral representation,” J. Mach. Learn. Res., vol. 22, no. 195, pp. 1–47, 2021. [12] G. Li, M. Müller, A. Thabet, and B. Ghanem, “DeepGCNs: Can GCNs go as deep as CNNs?” in IEEE/CVF International Conference on Computer Vision, 2019, pp. 9267–9276.

[13] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in International Conference on Learning Representations, 2017. [14] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in Advances in Neural Information Processing Systems, vol. 30, 2017. [15] J. Gasteiger, A. Bojchevski, and S. Günnemann, “Predict then propagate: Graph neural networks meet personalized PageRank,” in International Conference on Learning Representations, 2019. [16] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph attention networks,” in International Conference on Learning Representations, 2018. [17] S. Brody, U. Alon, and E. Yahav, “How attentive are graph attention networks?” in International Conference on Learning Representations, 2022. [18] Y. Shi, Z. Huang, S. Feng, H. Zhong, W. Wang, and Y. Sun, “Masked label prediction: Unified message passing model for semi-supervised classification,” in International Joint Conference on Artificial Intelligence, 2021, pp. 1548–1554. [19] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in Advances in Neural Information Processing Systems, vol. 29, 2016. [20] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in International Conference on Learning Representations, 2015. [21] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, “Optuna: A nextgeneration hyperparameter optimization framework,” in ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019. [22] P. Gaillard, G. Stoltz, and T. Van Erven, “A second-order bound with excess losses,” in Conference on Learning Theory, 2014. [23] D. W. Apley and J. Zhu, “Visualizing the effects of predictor variables in black box supervised learning models,” J. Roy. Statist. Soc. Ser. B, vol. 82, no. 4, pp. 1059–1086, 2020. [24] Z. Ying, D. Bourgeois, J. You, M. Zitnik, and J. Leskovec, “GNNExplainer: Generating explanations for graph neural networks,” in Advances in Neural Information Processing Systems, vol. 32, 2019. [25] F. X. Diebold and R. S. Mariano, “Comparing predictive accuracy,” J. Bus. Econ. Statist., vol. 13, no. 3, pp. 253–263, 1995. [26] W. K. Newey and K. D. West, “A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix,” Econometrica, vol. 55, no. 3, pp. 703–708, 1987. [27] D. Harvey, S. Leybourne, and P. Newbold, “Testing the equality of prediction mean squared errors,” Int. J. Forecasting, vol. 13, no. 2, pp. 281–291, 1997. [28] H. R. Künsch, “The jackknife and the bootstrap for general stationary observations,” Ann. Statist., vol. 17, no. 3, pp. 1217–1241, 1989. [29] P. R. Hansen, A. Lunde, and J. M. Nason, “The model confidence set,” Econometrica, vol. 79, no. 2, pp. 453–497, 2011. [30] Réseau de Transport d’Électricité (RTE), “éCO2mix: Open data on french electricity generation and consumption,” https://www.rte-france. com/en/eco2mix, 2024, regional half-hourly load and net-load data. [31] Météo-France, “SYNOP observational weather data,” https: //donneespubliques.meteofrance.fr, 2024, ground weather station observations. [32] L. Zhao, Y. Song, C. Zhang, Y. Liu, P. Wang, T. Lin, M. Deng, and H. Li, “T-GCN: A temporal graph convolutional network for traffic prediction,” IEEE Trans. Intell. Transp. Syst., vol. 21, no. 9, pp. 3848–3858, 2020. [33] J. Bai, J. Zhu, Y. Song, L. Zhao, Z. Hou, R. Du, and H. Li, “A3T-GCN: Attention temporal graph convolutional network for traffic forecasting,” ISPRS Int. J. Geo-Inf., vol. 10, no. 7, p. 485, 2021. [34] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in International Conference on Learning Representations, 2018. [35] Y. Seo, M. Defferrard, P. Vandergheynst, and X. Bresson, “Structured sequence modeling with graph convolutional recurrent networks,” in International Conference on Neural Information Processing, 2018. [36] K. Oono and T. Suzuki, “Graph neural networks exponentially lose expressive power for node classification,” in International Conference on Learning Representations, 2020.