Scalable estimation of VARMA models

Vector autoregressive moving-average (VARMA) models are powerful tools for multivariate time series, but their estimation has historically been limited to low-dimensional settings. The likelihood is non-convex, and standard estimation methods require a full pass over the entire dataset for every optimization step, making them computationally prohibitive for large series or high dimensions. This paper addresses this bottleneck by developing a scalable estimation framework that decouples the optimization cost from the series length.

The authors propose a framework that transforms the VARMA estimation problem into one that depends on the data only through fixed-size sufficient statistics. By using a partial-autocorrelation reparametrization, they guarantee stationarity and invertibility by construction, allowing for unconstrained optimization. The framework includes two primary point estimators: a regularized least-squares fit and a covariance-marginalized maximum-a-posteriori estimator. To further improve efficiency, they employ a Parseval (Fourier) identity to evaluate the loss at near-linear cost relative to the truncation length, effectively removing the dependence on the total series length $T$ during the optimization process.

The study demonstrates that this framework maintains statistical consistency, recovering the infinite-autoregressive representation of the true process at a near-parametric rate. Empirically, the estimators perform robustly in dimensions ranging from $d=10$ to $d=40$, where traditional conditional maximum likelihood estimation often fails due to non-invertible fits and divergent forecasts. The method consistently matches or outperforms standard baselines, including VAR, Bayesian-VAR, and sparse-VARMA, across diverse datasets such as retail demand, meteorology, and air quality. The framework also extends naturally to seasonal dynamics and exogenous regressors (VARMAX) without increasing the leading computational cost.