PaperlayerView exact revision on arXiv

Source-linked AI summary

Forecasting Time Series with VARMA Recursions on Graphs

Elvin Isufi, Andreas Loukas, Nathanael Perraudin, Geert Leus

arXiv:1810.08581v2eess.SPcs.SYecon.EMeess.SY

TL;DR

Forecasting high-dimensional multivariate series remains challenging, and the paper investigates how graph structure can help model their evolution. It combines approximate time-vertex stationarity with VAR/VARMA models and spectral-domain fitting. Experiments report benefits over classical and state-of-the-art graph-based alternatives, while the authors identify implications of the stationarity assumption and graph structure for accuracy as open questions.

  • Problem

    High-dimensional multivariate forecasting is difficult from both estimation and computational perspectives, while how to exploit graph structure for this task remains incompletely understood.

    Source: §I · p. 1

  • Method

    The paper combines approximate time-vertex stationarity with graph VAR/VARMA models, spectral-domain fitting, low-rank representations, and subset-node tracking.

    Source: §I · p. 2, §VII · pp. 13–14

  • Results

    Across three scenarios, the proposed approaches outperform classical univariate and multivariate ARMA techniques and state-of-the-art graph-based techniques.

    Source: §I · p. 2

  • Takeaways & Limitations

    The graph structure plays a major role in easing the computational complexity of high-dimensional multivariate models, particularly with limited training data.

    Source: §VII · p. 14

    Further work is needed to understand how the stationarity assumption affects prediction and how graph structure affects prediction accuracy.

    Source: §VII · p. 14

Abstract

from arXiv · show

Graph-based techniques emerged as a choice to deal with the dimensionality issues in modeling multivariate time series. However, there is yet no complete understanding of how the underlying structure could be exploited to ease this task. This work provides contributions in this direction by considering the forecasting of a process evolving over a graph. We make use of the (approximate) time-vertex stationarity assumption, i.e., timevarying graph signals whose first and second order statistical moments are invariant over time and correlated to a known graph topology. The latter is combined with VAR and VARMA models to tackle the dimensionality issues present in predicting the temporal evolution of multivariate time series. We find out that by projecting the data to the graph spectral domain: (i) the multivariate model estimation reduces to that of fitting a number of uncorrelated univariate ARMA models and (ii) an optimal low-rank data representation can be exploited so as to further reduce the estimation costs. In the case that the multivariate process can be observed at a subset of nodes, the proposed models extend naturally to Kalman filtering on graphs allowing for optimal tracking. Numerical experiments with both synthetic and real data validate the proposed approach and highlight its benefits over state-of-the-art alternatives.

I. INTRODUCTION

Forecasting high-dimensional multivariate series is challenging, and graph-based methods offer a way to use data connections to constrain prediction models. This work combines approximate time-vertex stationarity with VAR/VARMA models, spectral-domain fitting, and low-rank representations; experiments report benefits over classical and graph-based alternatives.

I. INTRODUCTION

  • Graph-structured forecasting faces estimation and computational challenges, especially when observations are fewer than time series; existing responses include factor models, shrinkage, and low-rank representations.Prior graph predictors often rely on deterministic-process assumptions or couple graph interactions to a fixed time step, potentially ignoring hidden node correlations.Source: §I · p. 1
  • Approximate joint stationarity preserves process statistics across time and graph localization, extending temporal stationarity to graph signals.The introduction notes that perfect stationarity is difficult to observe, while real datasets have shown degrees of graph/time stationarity that may be exploited for prediction.Source: §I · pp. 1–2
  • The proposed VAR and VARMA recursions model graph time series, while spectral-domain separation yields uncorrelated univariate series for fitting and low-rank representation further reduces estimation cost.The paper also introduces sub-graph tracking from subset measurements and connects it to Kalman filtering on graphs.Source: §I · p. 2
  • Across weather, dynamic-mesh, and epidemic scenarios, the proposed approaches outperform classical univariate and multivariate ARMA techniques and state-of-the-art graph-based techniques.The experiments also assess the joint-stationarity hypothesis.Source: §I · p. 2

A. Graph Signal Processing

  • A graph signal maps graph vertices to complex values, and the graph Fourier transform represents it in a spectral domain where graph filtering becomes element-wise multiplication.The graph frequencies are the eigenvalues of the graph Laplacian.Source: A. Graph Signal Processing · p. 2

B. Time-Vertex Signal Processing

The time-vertex framework applies joint spectral transforms to graph signals evolving over time, enabling filtering across both dimensions. Its ARMA recursion combines temporal input-output shifts with graph-Laplacian polynomial shifts and encompasses several established time-vertex models.

B. Time-Vertex Signal Processing

  • The joint Fourier transform diagonalizes time-vertex filtering, so filtering corresponds to element-wise multiplication in the joint frequency representation.The joint spectrum is indexed by graph frequency and temporal frequency.Source: B. Time-Vertex Signal Processing · p. 3
  • The time-vertex ARMA recursion combines current and past inputs with previous outputs, using temporal shifts for time filtering and powers of the graph Laplacian for graph filtering.The coefficients and polynomial orders determine the recursion structure.Source: B. Time-Vertex Signal Processing · p. 3
  • The general recursion encompasses separable time-vertex ARMA, finite-impulse-response time-vertex filtering, and the model of.Its polynomial graph shifts also support a distributed implementation.Source: B. Time-Vertex Signal Processing · p. 3

C. Stationarity on Networks

Graph wide-sense stationarity requires covariance structure aligned with the graph Fourier basis, while joint wide-sense stationarity extends the invariance condition across graph and time. The joint covariance is diagonalized in the joint Fourier basis and characterized by a joint power spectral density.

C. Stationarity on Networks

  • A zero-mean graph process is wide-sense stationary when its covariance is diagonalizable by the graph Fourier basis; its spectral power is the graph power spectral density.This expresses second-order statistical invariance with respect to graph shifts.Source: C. Stationarity on Networks · pp. 3–4
  • Joint wide-sense stationarity requires first- and second-order moments to be preserved jointly across graph and time, with covariance diagonalizable by the joint Fourier basis.The corresponding joint power spectral density is the joint Fourier transform of the autocorrelation function.Source: C. Stationarity on Networks · p. 4
  • For a joint wide-sense stationary process, the covariance is related to a joint time-vertex filter that can drive the process from zero-mean white noise.The filter magnitude squared equals the joint power spectral density.Source: C. Stationarity on Networks · p. 4

III. MODELING OF NETWORK PROCESSES

The paper develops causal graph VARMA models for forecasting jointly stationary network processes after identifying the computational and noncausal drawbacks of a general joint-filter model. Polynomial graph VAR provides a distributed alternative with lower per-iteration cost under the stated sparse-graph setting.

A. Joint Non-Causal Models on Graphs

  • The general joint noncausal model can represent a jointly stationary process, but its computational burden and noncausality make it unsuitable for timely forecasting from process history.The paper introduces causal models to address both issues.Source: A. Joint Non-Causal Models on Graphs · p. 4

B. Joint Causal Models on Graphs

  • The graph-VARMA model is a causal graph-process model that can describe all causal jointly stationary processes for sufficiently large P and Q.Its per-iteration cost is O((P + Q)N^2), comprising P + Q matrix-vector multiplications.Source: B. Joint Causal Models on Graphs · p. 4
  • The authors omit the moving-average part from the polynomial graph model because polynomial graph VARMA was observed to produce unstable predictors or, at best, larger MSE than graph polynomial-VAR.The text distinguishes this observation from the behavior of the graph-Fourier-fitted G-VARMA model.Source: B. Joint Causal Models on Graphs · p. 4
  • For |E| ≪ N^2, graph polynomial-VAR has distributed per-iteration complexity O(Lmax|E|), below graph-VARMA's O((P + Q)N^2) cost.Unlike the restricted model of, its graph polynomial order can vary with time index and is not constrained by Lp ≤ p; the paper says this permits more hidden node interactions and yields better predictors.Source: B. Joint Causal Models on Graphs · p. 5

C. Optimal predictor

The paper derives a one-step predictor from the conditional expectation given past observations, establishes its MSE optimality for fixed parameters, and outlines a structured parameter-fitting strategy.

C. Optimal predictor

  • The one-step predictor is the conditional expectation of the process given its past realizations.For the GP-VAR model, it follows by setting Q = 0 and using the specified graph-filter coefficients.Source: C. Optimal predictor · p. 5

D. Mean Square Error Analysis

  • For fixed parameters of the joint causal model, the one-step predictor achieves the smallest achievable MSE.The analysis examines how the predictor parameters affect MSE; the proposition's displayed expression is not included in the supplied passage.Source: D. Mean Square Error Analysis · p. 5

IV. MODEL FITTING: PARAMETER DESIGN

  • Parameter design fits G-VARMA kernels or GP-VAR scalars, using GFT decoupling and low-rank estimation to trade prediction accuracy against computational cost.The section also considers estimating GP-VAR parameters from graph-process statistics and joint frequency smoothing.Source: §IV · p. 5

A. G-VARMA fitting

G-VARMA fitting uses graph-frequency decoupling to simplify estimation, then restricts fitting to selected frequencies for further savings, with a stated accuracy trade-off.

A. G-VARMA fitting

  • GFT decoupling replaces a nonlinear system with N smaller, uncorrelated fitting problems, each equivalent to fitting a univariate temporal ARMA model.Each frequency problem has T equations and P + Q unknowns, and can be solved using established methods such as Gauss-Newton.Source: A. G-VARMA fitting · pp. 5–6
  • The full G-VARMA approach is better suited to small- or medium-sized graphs because Laplacian eigendecomposition adds estimation cost; low-rank fitting trades accuracy for savings.The low-rank option requires fitting only K temporal ARMA series.Source: A. G-VARMA fitting · p. 6
  • The optimal K-rank approximation retains graph-frequency components indexed by the largest diagonal elements of the covariance in the graph spectral basis.For a JWSS process, rotating by the graph Laplacian eigenvectors yields the best rank-K approximation.Source: A. G-VARMA fitting · p. 6
  • Fitting only K selected frequencies avoids a full Laplacian eigendecomposition and reduces G-VARMA estimation cost by N/K.The paper notes that real data can have K ≪ N; only the eigenvectors for the selected frequencies are needed.Source: A. G-VARMA fitting · p. 6

B. GP-VAR fitting

GP-VAR fitting estimates scalar coefficients from prediction error or process autocorrelations, while spectral smoothing can improve fitting in low-sample regimes by exploiting neighboring-frequency similarity.

B. GP-VAR fitting

  • GP-VAR fitting estimates scalar coefficients ψ_l,p by minimizing prediction MSE or solving a convex problem based on process autocorrelations.The coefficients can alternatively be estimated with a Yule-Walker least-squares approach.Source: B. GP-VAR fitting · pp. 6–7
  • The graph-imposed coefficient structure makes GP-VAR parameter estimation computationally easier than classical VAR and more robust to larger multivariate dimensions.The autocorrelation estimates can incorporate shrinkage or random-matrix methods for low-sample regimes, and estimators for heavy-tailed non-Gaussian data.Source: B. GP-VAR fitting · pp. 6–7
  • A Yule-Walker strategy can estimate GP-VAR coefficients from subsampled measurements while keeping the optimization computationally lighter.The passage describes combining it with the approach in.Source: B. GP-VAR fitting · p. 7

C. Spectral Smoothing

  • Spectral smoothing convolves the joint power spectral density with smoothly decaying windows; adjacent-frequency similarity can reduce estimation error in low-sample settings.The Gaussian window width controls a bias-variance trade-off, and the paper says smoothing will often lower estimation error when abrupt spectral changes reflect noise or finite-sample inaccuracies.Source: C. Spectral Smoothing · p. 7

V. TRACKING ON A SUBSET OF NODES

The tracking method combines a G-VARMA prediction with current measurements at sampled nodes to estimate a chosen set of nodes, with error depending on both sampling and model residuals.

V. TRACKING ON A SUBSET OF NODES

  • The tracker uses a G-VARMA prior and measurements from a subset of nodes to estimate the process on a target set.The target may be all nodes, the unsampled complement, or a combination; the example uses uniform random sampling.Source: §V · pp. 7–8
  • The estimate combines the prior prediction with inverse filtering on observed vertices and spreads the result onto the set of interest.The paper identifies the prediction as the a priori estimate and the updated value as the a posteriori estimate.Source: §V · p. 8
  • Tracking MSE depends on the sampling set and the G-VARMA residual covariance.The authors use uniform random sampling to demonstrate the tracker, although sparse-sensing strategies could improve MSE by 1–2 dB over random uniform sampling at higher complexity.Source: §V · p. 8
  • The Kalman-filter reformulation can incorporate another predictor, including GP-VAR, by supplying its a priori prediction error covariance.The passage states that GP-VAR's covariance can be derived from its parameter-fitting expression.Source: §V · p. 8

VI. NUMERICAL RESULTS

Across four scenarios, the proposed graph predictors perform best in Molene and exploit graph parameterization effectively in NOAA. Results vary by scenario and horizon: graph-based methods are stable on the dog mesh, while epidemic forecasts favor the proposed approaches beyond one step.

A. Considered scenarios and data pre-processing

  • The evaluation tests four scenarios: Molene weather data, NOAA hourly temperatures from 109 U.S. stations in 2010 (8,759 hours), a walking-dog mesh of 251 nodes over 59 time steps, and simulated SI epidemic diffusion over 100 international airports.Source: §VI · p. 8, A. Considered scenarios and data… · p. 9

B. Prediction performance

  • For all scenarios, training, cross-validation, and test data use 35%, 15%, and 50%, respectively, and model fitting minimizes average rNMSE over five-step-ahead prediction.Parameters are selected on validation data and the selected model is refit on all in-sample data before out-of-sample testing.Source: §VI · p. 8, B. Prediction performance · p. 9
  • In Molene, G-VARMA and GP-VAR achieve the best predictive performance; in NOAA, standard VAR training repeatedly fails while graph-based methods are trainable and the proposed methods yield the minimum performance.The NOAA disjoint ARMA performs close to the proposed algorithms, while the graph topology is held fixed for direct comparison.Source: B. Prediction performance · pp. 9–10, Figure 4 · p. 10
  • On the dog mesh, classic VAR predictions are unstable while graph-based approaches are more stable; for epidemic forecasts, proposed methods outperform adaptive methods beyond one step, but univariate ARMA catches up from three steps ahead.For the epidemic scenario, graph-based methods have similar one-step rNMSE, except RGP-VAR performs worse; from three steps ahead, univariate ARMA has similar rNMSE to the proposed methods.Source: B. Prediction performance · pp. 10–11, Figure 5 · p. 11
  • For the dog mesh, three-step prediction is worse than one-step prediction, with larger errors around the most mobile body parts, including the legs and tail.The three-step G-VARMA prediction matches the dog position well on most nodes.Source: B. Prediction performance · p. 11, Figure 6 · p. 12
  • On Molene, graph-based methods are more robust than standard VAR when in-sample data are scarce; VAR could not be trained with 10% in-sample data.The passage attributes robustness to fewer parameters and notes a role for spectral smoothing when training data are limited.Source: Figure 7 · p. 12
  • In the epidemic scenario, G-VARMA follows one-step I-state evolution with error up to 4%, but performance degrades at three steps ahead, partly due to outlier nodes.The passage attributes higher rNMSE in this scenario to abrupt transitions when infected targets recover every 12 days.Source: B. Prediction performance · pp. 11–12

C. In-sample data effect

On Molene, graph-based methods are more robust when in-sample data are scarce, while GP-VAR improves on G-VARMA at 10%. The paper also reports a model-complexity and training-data trade-off for its proposed methods.

C. In-sample data effect

  • The setup kept cross-validation data at 30% of the in-sample data and selected parameters to minimize one-step-ahead rNMSE.Source: C. In-sample data effect · p. 12
  • The reported outcome may improve with other model parameters.Source: C. In-sample data effect · p. 12
  • 0.13 maximum rNMSE loss accompanies a speed-up of 16 on NOAA for low-rank G-VARMA estimation.Fig. 8 plots gained speed-up against prediction rNMSE loss with both axes on a log scale.Source: Figure 8 · p. 12
  • At 10% in-sample data, GP-VAR improves on G-VARMA due to its fewer DoFs.The multivariate VAR suffers most in low-data regimes; graph-based methods are more robust, which the authors associate with fewer parameters and spectral smoothing.Source: C. In-sample data effect · p. 12
  • The proposed methods provide a topology-aware trade-off between model complexity and training-data amount; VAR was unstable in other scenarios.The authors omit those scenario results because VAR, the comparison recursion, did not yield stable predictors.Source: C. In-sample data effect · p. 12

D. Low-rank model effect

On Molene, G-VARMA tracking outperforms the alternatives even with only eight current measurements; low-rank estimation reveals a trade-off between prediction accuracy and computation. The conclusion attributes computational benefits to graph structure while noting that stationarity assumptions and graph effects on accuracy need further study.

E. Tracking performance

  • On Molene, G-VARMA outperforms the other tracking approaches with only N/4 = 8 current samples; Wiener inpainting is suboptimal because it omits historical data.Fig. 9 plots tracking rNMSE against observation-set size.Source: Figure 9 · p. 13, E. Tracking performance · p. 13

D. Low-rank model effect

  • Ignoring increasing portions of the process JPSD trades prediction accuracy for computational speed-up; NOAA is most sensitive and the SI epidemic data set least sensitive.The authors attribute NOAA's greater sensitivity to mismatch between the graph and stationarity assumptions, and point to JPSD sparsity as an opportunity to reduce estimation costs.Source: D. Low-rank model effect · pp. 12–13

E. Tracking performance

  • On NOAA, G-VARMA prediction and tracking have lower rNMSE than KKF, while Wiener inpainting performs worst because it does not use the temporal dimension.Cross-validation slightly improves KKF performance.Source: E. Tracking performance · p. 13

VII. CONCLUSION

  • The conclusion describes graph-spectral coefficient design as reducing VAR and VARMA degrees of freedom, with an optimal low-rank representation further reducing estimation costs.A graph tracker is also introduced for measurements available at a subset of nodes.Source: §VII · pp. 13–14
  • The experiments show benefits over classical multivariate forecasting and state-of-the-art graph-based approaches, while graph structure eases high-dimensional model computation.Source: §VII · p. 14
  • Further work is needed to understand how the stationarity assumption affects prediction and how graph structure affects prediction accuracy.Source: §VII · p. 14

APPENDIX

The appendix completes a one-step prediction MSE derivation and establishes an optimal low-rank representation by retaining the largest graph-spectral components. It also identifies zoomed-in versions of two figures.

APPENDIX

  • For fixed coefficients, the one-step prediction MSE depends only on the unknown innovation at time t and is the smallest achievable MSE.Source: APPENDIX · p. 14
  • The optimal rank-K representation retains the K largest components of g(λ); any other rotation gives a larger approximation error.Source: APPENDIX · p. 14

ZOOMED IN FIGURES

  • The appendix provides zoomed-in versions of Fig. 4c and Fig. 10 as Figs. 11 and 12, respectively.Source: ZOOMED IN FIGURES · p. 14

Related papers

Loading 1810.08581v2…