Scalable estimation of VARMA models

Scalable estimation of VARMA models Daniel Paulin

daniel.paulin@ntu.edu.sg

College of Computing & Data Science Nanyang Technological University Singapore

Vı́ctor Elvira

victor.elvira@ed.ac.uk

arXiv:2608.06340v1 [stat.ML] 6 Aug 2026

School of Mathematics University of Edinburgh Edinburgh, UK

Abstract Vector autoregressive moving-average (VARMA) models have long been considered impractical beyond moderate dimensions: the likelihood is non-convex, the parametrization is identified only up to equivalence, and every evaluation costs a pass over the entire series. Yet their moving-average term captures with a few parameters what a pure autoregression matches only with many lags. We introduce an estimation framework that removes this computational barrier: each optimization iteration is independent of the series length T . The framework combines a partial-autocorrelation reparametrization that guarantees stationarity and invertibility by construction, penalized least-squares and covariance-marginalized losses that depend on the data only through fixed-size sufficient statistics, and a Parseval (Fourier) evaluation whose cost is near-linear in the truncation length. This yields two point estimators: a regularized least-squares fit and a covariancemarginalized maximum-a-posteriori estimator; in both cases, using Gaussian priors on the reparametrized VARMA coefficients with separate standard deviations for diagonal and non-diagonal coefficients. We further prove that both estimators recover the infiniteautoregressive representation of the true process at a near-parametric rate in fixed dimension, so the truncation does not introduce asymptotic bias: both estimators still converge to the true dynamics at a near-parametric rate. The same sufficient-statistics machinery extends, at the same leading cost, to seasonal dynamics, exogenous regressors (VARMAX), and rolling-window refits. Empirically, the estimators stay close to the oracle forecast error from d = 10 to d = 40 (where the classical conditional-MLE implementation returns non-invertible fits whose forecasts diverge) and match or beat VAR, Bayesian-VAR, component-wise ARMA, and sparse-VARMA baselines on retail-demand, meteorological, and air-quality data. These results bring likelihood-based VARMA estimation, at a periteration cost independent of the series length, to the problem sizes at which practitioners have until now relied on VAR models. Keywords: VARMA, VARMAX, seasonal VARMA, high-dimensional time series, scalable estimation, partial autocorrelation, forecasting

1

Paulin and Elvira

1 Introduction Multivariate time series arise throughout science, and vector autoregressive moving average (VARMA) processes are a natural model class for them (Box et al., 2015; Reinsel, 2003; Lütkepohl, 2005; Tsay, 2014; Hannan and Deistler, 2012; Kilian and Lütkepohl, 2017). Joint modelling captures cross-variable dependence that componentwise analyses miss, and the class is closed under marginalization. As a result, the observed components of a partially observed vector autoregressive (VAR) model follow a VARMA process (see e.g. Akaike, 1974, or Section 3.7 of Tsay, 2014). The estimation of vector autoregressive (VAR) processes is comparatively simple; it reduces to solving large linear systems of equations and makes VARs applicable to highdimensional series. VARMA estimation is however not easy to be applied. Despite much effort, its non-convex likelihood and identifiability problems have kept it from scaling, and the standard packages (Tsay et al., 2022; Seabold and Perktold, 2010) are often unreliable beyond about ten dimensions. In the dense regime, existing methods face important challenges. The likelihood-based methods (state-space Kalman-filter ML, the Whittle likelihood, Bayesian VARMA) are statistically efficient but scale poorly in the dimension or the series length, while the scalable ones (Hannan–Rissanen two-stage least squares, subspace identification) replace the likelihood by a cheaper convex criterion, remaining consistent but statistically less efficient. We review them in Section 2.1. To the best of our knowledge, there is still missing a one-stage, dense, likelihood-motivated estimator whose per-iteration cost does not grow with T and that stays tractable as the number of series grows. We address that gap in this paper. Faced with these difficulties, practitioners have long favoured VARs over VARMAs even where a moving-average term would be more parsimonious (Lütkepohl, 2006, Section 1.1). In this paper, we develop a dense VARMA framework with computational efficiency and sound theoretical guarantees that can scale up to high dimension. We minimize a least-squares-type prediction loss, written in reparametrized coordinates (without loss of generality) in which stationarity and invertibility hold by construction. The loss depends on the data only through fixed-size sufficient statistics, so that each optimization iteration is independent of the series length T . As the experiments show, the resulting estimators are competitive in accuracy with the maximum likelihood estimator (MLE) on small models and remain tractable at dimensions where MLE is not. Although we develop the method for the VARMA model, the same sufficient-statistics framework accommodates, without changing its leading computational cost, two extensions that matter in applications: a multiplicative seasonal VARMA for periodic data (Section 5.5), and a VARMAX model that conditions the dynamics on exogenous regressors. The exogenous terms enter the very same pre-computed sufficient statistics, so VARMAX inherits the same T -independent gradient cost. This underpins a central part of our empirical study: one of our three real-data benchmarks forecasts retail demand (units) conditional on the observed price path (Section 7.2), a genuine VARMAX problem in which the movingaverage dynamics and the exogenous price together prove decisive. Table 1 summarises the per-gradient computational cost of these approaches. Both existing MLE-based methods scale linearly in T per gradient evaluation, since each optimizer step must touch the entire training history through either the recursive innovations or the 2

Scalable estimation of VARMA models

Table 1: Per-gradient computational complexity of VARMA(p, q) fitting methods. Our methods replace the per-gradient O(T ) cost of MLE-based approaches by a one-time prepass that scales with T , after which gradient evaluations are independent of T . The Fourier (Parseval) evaluation of S 2 further reduces the dependence on the truncation length m from m2 to m log m, making large-m regimes tractable. The MLE bounds are representative dense companion-form costs, not lower bounds for every state-space algorithm. In the per-gradient column, the final O((p2 +q 2 )d3 ) term for our methods is the T - and mindependent PAR-reparameterisation overhead; the preceding, data-dependent terms are those the Fourier evaluation reduces. Method Existing approaches MLE, recursive cond. likelihood (Tsay, 2014) MLE, state-space + Kalman filter (Durbin and Koopman, 2012; Metaxoglou and Smith, 2007) This paper: truncated VARMA(p, q, m) Sufficient-statistics evaluation (Section 4.1) + Fourier (Parseval) evaluation of S 2 (Section 6)

One-time pre-pass

Per-gradient evaluation

none

O T d3 (p + q)2

none

O T d3 max(p, q + 1)3

O T m d2

O T log m d2

O m2 d3 + (p2 +q 2 )d3

O m log m d2 + m d3 + (p2 +q 2 )d3

Kalman filter. Our truncated VARMA approach absorbs the T -dependence into a one-time pre-pass that produces a fixed-size sufficient-statistic summary, after which every gradient evaluation is independent of T . The Fourier-based variant additionally avoids materialising the m2 d2 second-moment block matrix and reduces the per-gradient cost to nearly linear in m, making large truncation lengths computationally tractable. The figures in Table 1 are the dominant, data-dependent part of each gradient; the PAR reparametrisation adds a further O((p2 +q 2 ) d3 ) per evaluation, and the multiplicative seasonal VARMA of Section 5.5 evaluates at the same leading T -independent cost, with its full dependence on the orders (p, q, P, Q), the period s, and (m, d) given precisely in (61). Contributions. (i) A scalable estimation framework for dense VARMA whose periteration optimization cost is independent of the series length T : the fixed-size sufficient statistics that summarize a VAR regression are shown to drive the nonlinear VARMA problem as well. (ii) Statistical guarantees: both estimators recover the infinite-VAR representation at a near-parametric rate in fixed dimension, so scalability does not sacrifice consistency. (iii) Practical extensions including (a) seasonal VARMA, exogenous regressors (VARMAX), and (b) incremental rolling-window refits, both at the same leading cost, validated on retail-demand, meteorological, and air-quality benchmarks. The rest of the paper is organized as follows. Section 2.1 introduces VARMA processes and reviews the two standard MLE-based fitting strategies and their computational bottleneck (Section 2.2), the basic least-squares formulation (Section 2.3), the identifiabil3

Paulin and Elvira

ity issues that we do not address (Section 2.4), and the existing software (Section 2.5). Section 3 defines the truncated VARMA(p, q, m) model and the stationarity/invertibility constraints it must obey. Sections 4 and 5 develop the two pointwise loss functions we use to fit the truncated model, the regularised least-squares (RLS) loss and the Σ-marginalised MAP loss, written directly in the unconstrained PAR coordinates of Section 2.1 so that the stationarity and invertibility constraints are enforced by construction, and show that both depend on the data only through pre-computed sufficient statistics; Section 5.6 proves that, for fixed dimension, both estimators recover the infinite-VAR representation of the true process at a near-parametric rate as T and m grow jointly. Both objectives are min(p+q)d2 Θ imised by L-BFGS on the unconstrained PAR parameters (M Φ , with no 1:p , M 1:q ) ∈ R 2 2 projection step required. Section 6 replaces the materialised m d second-moment block by a Fourier/Parseval representation, bringing the per-gradient evaluation cost down to O(m log m d2 + m d3 ). Section 7 presents our experimental results. Appendix E gives the LG-SSM equivalence used in Section 2.2, and Appendix F shows how the sufficient statistics, and hence the entire fit, can be updated incrementally as the training window slides forward, so that refits cost O(∆ m d2 ) per shift.

2 Background 2.1 VARMA models For t ∈ Z, let Zt =

p X

Φj Z t−j + at −

j=1

q X

Θj at−j

(1)

j=1

be a stationary d-dimensional VARMA model with mean 0. Here Z t ∈ Rd , Φj ∈ Rd×d for 1 ≤ j ≤ p, Θj ∈ Rd×d for 1 ≤ j ≤ q, at ∈ Rd is a vector white-noise process with mean 0 and covariance matrixP Σ ∈ Rd×d (independent over time). P Let Φ(z) := I − pj=1 Φj z j and Θ(z) := I − qj=1 Θj z j be the matrix-valued AR and MA lag polynomials (so that Φ(0) = Θ(0) = I), in terms of which the model (1) reads Φ(B)Z t = Θ(B)at with B the backshift operator (BZ t = Z t−1 ); here z ∈ C. The condition required for the stationarity of this model can be established in the spectral domain as follows. Assumption 1 (Stationarity and invertibility) We assume that det Φ(z) ̸= 0 and det Θ(z) ̸= 0 if |z| ≤ 1, z ∈ C,

(2)

i.e., all complex roots of det Φ(z) = 0 and det Θ(z) = 0 are outside the unit circle. Companion matrices and spectral radius reformulation. Assumption 1 admits an equivalent characterisation in terms of the spectral radii of the AR and MA block companion matrices. Define CΦ ∈ Rpd×pd , CΘ ∈ Rqd×qd as     0 I 0 ... 0 0 I 0 ... 0 0 0 0 I ... 0  0 I ... 0        .. .. .. . . ..  , C .. .. . . ..  . CΦ :=  ... := (3)   Θ . .  . .  . . . .   .  0 0 0 0 ... I  0 0 ... I  Φp Φp−1 . . . . . . Φ1 Θq Θq−1 . . . . . . Θ1 4

Scalable estimation of VARMA models

The non-zero eigenvalues of CΦ (resp. CΘ ) are the reciprocals of the roots of det Φ(z) (resp. det Θ(z)), so the root-outside-the-unit-circle conditions of Assumption 1 are equivalent to spr(CΦ ) < 1 and spr(CΘ ) < 1, where spr(·) denotes the spectral radius. Section 3.2 returns to the companion matrices to derive the geometric decay of the infinite-VAR coefficients. Partial-autocorrelation reparametrisation.

The stable set

{(Φ1:p , Θ1:q ) : spr(CΦ ) < 1, spr(CΘ ) < 1} is open, non-convex, and bounded by a complex algebraic variety, making constrained optimisation over it cumbersome. We will exploit the well-known fact that this set admits a smooth bijection from unconstrained Euclidean space, via the multivariate Levinson recursion in partial-autocorrelation form (Ansley and Kohn, 1986; Heaps, 2023). Reparametrisations of this kind have been used for estimation before: Roy et al. (2019) compute exactlikelihood estimates of causal invertible VARMA models by unconstrained optimization in closely related coordinates, and Heaps (2023) enforces stationarity of Bayesian VARs through a prior in these coordinates. Our contribution lies not in the parametrisation but in the loss evaluated on it: the truncated sufficient-statistic criteria of Sections 4–5, whose per-iteration cost is independent of T . The construction is identical for the AR and MA sides, so we first describe a generic transformation acting on a stack of r unconstrained lag matrices, then express the full VARMA reparametrisation as a factorised product. Generic transformation τr . For any r ∈ N, we define a smooth bijection τr : Rr×d×d −→



A1:r ∈ Rr×d×d : spr(CA ) < 1 ,

(4)

where CA ∈ Rrd×rd is the block companion matrix of A1:r defined analogously to (3). Given an unconstrained stack M 1:r = (M 1 , . . . , M r ) ∈ Rr×d×d , the map τr produces A1:r = τr (M 1:r ) in two stages. Stage 1 (per-lag contraction). Write chol(Σ) for the Cholesky factor of a positivedefinite matrix Σ, the unique lower-triangular matrix with positive diagonal entries satisfying chol(Σ) chol(Σ)⊤ = Σ. We state both stages in the Cholesky form in which they are implemented. For each s = 1, . . . , r, let LA,s := chol(I + M ⊤ s M s ) and define P s := M s L−⊤ A,s .

(5)

⊤ −1 2 2 Then P ⊤ s P s = I − (LA,s LA,s ) , whose eigenvalues are σi (M s ) /(1 + σi (M s ) ) < 1, so −1 σmax (P s ) < 1; conversely, given any P with σmax (P ) < 1, setting LA = chol(I − P ⊤ P ) and M = P L⊤ A recovers the unique preimage. Hence (5) is a smooth diffeomorphism from Rd×d onto the open ball of strict contractions {P ∈ Rd×d : σmax (P ) < 1}. Stage 2 (Whittle–Ansley–Kohn recursion, innovation-anchored). The recursion is anchored at a fixed innovation covariance, taken to be I d , and is carried out on Cholesky factors Lfs , Lbs of the forward and backward partial-prediction covariances Σfs = Lfs (Lfs )⊤ , Σbs = Lbs (Lbs )⊤ . The anchor is first propagated backward: set Lfr = I d and, for s = r, r − 1, . . . , 1,
−1 Lfs−1 = Lfs chol I − P s P ⊤ . (6) s

5

Paulin and Elvira

The terminal factor determines the stationary variance Γ0 = Lf0 (Lf0 )⊤ consistent with innovation covariance I d , and the backward chain starts from Lb0 = Lf0 . The coefficient pass then runs forward: for s = 1, . . . , r, A(s) = Lfs−1 P s (Lbs−1 )−1 , s e (s) = Lb P ⊤ (Lf )−1 , A s s−1 s s−1 (s)

(s−1)

Aj = Aj

(s−1)

e − A(s) s As−j ,

j = 1, . . . , s − 1,

(7)

e (s) = A e (s−1) − A e (s) A(s−1) , A j j s s−j
Lbs = Lbs−1 chol I − P ⊤ s Ps , the forward factors Lfs−1 being those produced by the backward pass (6). All factors are ⊤ well defined, since I − P s P ⊤ s ≻ 0 and I − P s P s ≻ 0 when σmax (P s ) < 1, and a product of lower-triangular matrices with positive diagonals is again one, so every Lfs and Lbs is a genuine Cholesky factor. The two passes are mutually consistent by construction: the f f ⊤ implied forward variance updates Σfs = Lfs−1 (I − P s P ⊤ s )(Ls−1 ) terminate at Σr = I d , (r)

(r)

the innovation anchor. The output is τr (M 1:r ) := (A1 , . . . , Ar ); the intermediate stack e (r) encodes the time-reversed (backward-AR) coefficients and is not returned. Cholesky A 1:r factorization is smooth (indeed real-analytic) on the positive-definite cone, so τr is smooth, including at M = 0, where every factor equals I d and τr is locally the identity. The construction is the matrix partial-autocorrelation parametrisation of Ansley and Kohn (1986) and Heaps (2023), the multivariate counterpart of the scalar partial-autocorrelation diffeomorphism of Barndorff-Nielsen and Schou (1973); it has been used for VARMA coefficients before by Roy et al. (2019). The variant that uses principal symmetric square roots in place of Cholesky factors differs from (7) only by a two-sided orthogonal change of the standardized coordinates of each P s and has the same image; we work with the Cholesky form throughout because it is the implemented one: it costs a triangular factorization and triangular solves per stage, and its automatic differentiation is well behaved at coincident eigenvalues (as at M = 0), where eigendecomposition-based square-root implementations are not, although the principal square-root map itself is smooth there. Anchoring at the innovation covariance, rather than at a fixed marginal variance, is what makes the image the full stable set. Proposition 1 (PAR bijection) For every r ≥ 1 and d ≥ 1, the map τr defined by (5), (6) and (7) is a smooth bijection from Rr×d×d onto the open stable set {A1:r ∈ Rr×d×d : spr(CA ) < 1}, with smooth inverse. The proof, in Appendix A, verifies that the backward pass (6) is the unique solution of the anchoring equations, that every contraction stack produces a Schur-stable coefficient stack (a matrix form of the classical Schur–Levinson stability argument, via the J-unitarity of the Cholesky-normalized lattice section), that any stable coefficient stack arises from a unique contraction stack (through the stationary autoregression with innovation covariance I d ), that the forward pass reproduces the coefficients, and that both directions are smooth. The anchoring is what carries the full stable set: with a fixed marginal variance Γ0 = I d 6

Scalable estimation of VARMA models

instead, a stationary VAR(1) with coefficient A1 has innovation covariance I d − A1 A⊤ 1 ≻ 0, so representability forces σ (A ) < 1, a strict subset of the stable set (for instance max 1   1/2 2 A1 = 0 1/2 is stable with σmax > 1). Joint VARMA reparametrisation. Applying τr separately to the AR side (with r = p) and the MA side (with r = q) yields the factorised reparametrisation  2 T : R(p+q) d −→ (Φ1:p , Θ1:q ) : spr(CΦ ) < 1, spr(CΘ ) < 1 ,

Θ Φ Θ T MΦ (8) 1:p , M 1:q := τp (M 1:p ), τq (M 1:q ) , p×d×d q×d×d where M Φ and M Θ are unconstrained matrices. Because τp and 1:p ∈ R 1:q ∈ R τq are independent diffeomorphisms (Proposition 1), T is itself a smooth bijection between 2 the unconstrained space R(p+q)d and the open stable set on the right-hand side of (8). We exploit this in Sections 4.3 and 5 to write the RLS and Σ-marginalised MAP losses (p+q)d2 Θ as unconstrained functions of (M Φ , eliminating the need for projection 1:p , M 1:q ) ∈ R steps or augmented-Lagrangian penalties during L-BFGS optimisation, as in Roy et al. (2019). As explained in (Box et al., 2015, Section 14.4.1), under Assumption 1, we have

Z t :=

∞ X

Πj Z t−j + at ,

(9)

j=1

where the matrices Πj satisfy the recursions Πj = Θ1 Πj−1 + Θ2 Πj−2 + . . . + Θq Πj−q + Φj , Πj = Θ1 Πj−1 + Θ2 Πj−2 + . . . + Θq Πj−q

j = 1, 2, . . . , p,

for j > p,

(10) (11)

where Π0 = −I, Πj = 0 for j < 0. 2.2 Maximum-likelihood estimation and its computational cost The natural estimator of the VARMA(p, q) parameters (Φ1:p , Θ1:q , Σ) is the maximumiid

likelihood estimator (MLE) under the Gaussian assumption at ∼ N (0, Σ). Two standard ways of evaluating the Gaussian likelihood and its gradient are used in the literature, both with per-iteration cost that scales linearly in the sample size T . Recursive evaluation of the conditional likelihood. The first approach evaluates the conditional likelihood p(Z 1:T | Φ1:p , Θ1:q , Σ) recursively in t by computing the onestep-ahead innovations Ẑ t+1|1:t from the VARMA recursion; see Section 14.4.5 of Box et al. (2015) or Section 5.3.1 of Reinsel (2003) for the standard formulas. Tsay (2014) implements
a conditional-likelihood gradient descent in which each iteration costs O d3 T (p + q)2 flops; the dependence on T comes from the inner forward pass, and the (p + q)2 factor from the recursive structure of the gradient. State-space evaluation via Kalman filter. The second approach exploits the fact that every stationary invertible VARMA(p, q) admits an equivalent finite-dimensional linear Gaussian state-space (LG-SSM) representation; we give an explicit construction in Appendix E. The Gaussian likelihood and its gradient with respect to (Φ1:p , Θ1:q , Σ) can then 7

Paulin and Elvira

be evaluated by running a Kalman filter on the equivalent LG-SSM. Durbin and Koopman (2012) provides a comprehensive treatment of state-space likelihood methods, and Metaxoglou and Smith (2007) apply this machinery specifically to VARMA estimation. A standard companion-form LG-SSM realisation of a VARMA(p, q) has state dimension dx = d max(p, q + 1)

(12)

(the minimal McMillan degree can be smaller under pole–zero cancellations or reduced-rank coefficient structure, but this companion realisation is the one used in practice), giving a representative per-likelihood-evaluation cost dominated by the O(d3x ) Kalman covariance update at each time step:

O T d3x = O T d3 max(p, q + 1)3 . (13) Common bottleneck. Both approaches share two characteristics that make them expensive in high dimensions. First, every gradient evaluation must touch all T observations through the recursion (or the filter), so cost scales linearly in T . Second, the cubic dependence on d (or on dx ) forbids scaling much beyond a few tens of dimensions even with modest T . In practice we observe that for d ≳ 10 MLE-based fits become impractical, taking hours per fit and exhibiting frequent convergence failures from the non-convex likelihood surface. 2.3 Least squares estimators When the goal is forecasting rather than parameter identification, a common alternative is to fit the VARMA model by minimising the one-step-ahead squared prediction error. With Ẑ t+1|1:t = E[Z t+1 | Z 1:t = z 1:t ], where the expectation is with respect to the VARMA(p, q) predictive distribution with parameters (Φ1:p , Θ1:q , Σ) treated as a deterministic tuple. We let L∞ (Φ1:p , Θ1:q , Σ) :=

T −1 X

∥Ẑ t+1|1:t − z t+1 ∥2 .

(14)

t=0

For Σ = I and Gaussian vector white-noise process at , L∞ coincides with the conditional Gaussian negative log-likelihood up to a constant (see Section 14.5 of Box et al. (2015)). Direct evaluation of L∞ is, however, no cheaper than the maximum-likelihood approach: the recursive structure of Ẑ t+1|1:t forces every gradient evaluation to use all T observations, with the same O(T ) scaling as the recursive MLE. 2.4 Identifiability of VARMA parameters Even at the population level, the VARMA(p, q) parameters (Φ1:p , Θ1:q , Σ) are not in general uniquely determined by the observed process: distinct (Φ1:p , Θ1:q , Σ) tuples can produce the same transfer function Φ(L)−1 Θ(L), and thus the same Wold representation and the same predictive distribution of Z t . The classical fix is to restrict attention to a canonical form (e.g. the Echelon form) in which the parametrisation is uniquely identified by a tuple of so-called Kronecker indices, but locating those indices is itself a combinatorial search over an 8

Scalable estimation of VARMA models

O((p + q)d )-size set, and therefore scales poorly with d. Düker et al. (2025) review this and related identifiability questions; Wilms et al. (2023) propose a convex two-stage procedure that first fits a truncated VAR to estimate residuals and then estimates (Φ1:p , Θ1:q , Σ) via penalised least squares with a sparsity-inducing regulariser, side-stepping the Kroneckerindex search at the cost of relying on the ℓ1 structure of the chosen penalty. In this paper we do not address VARMA identifiability. Our focus is instead on building forecasts for high-dimensional time series, for which any tuple (Φ1:p , Θ1:q , Σ) in an equivalence class of observationally indistinguishable parametrisations gives identical predictions of Z t , and is therefore equally good. 2.5 Existing software for VARMA fitting VARMA is implemented in several time-series packages, including statsmodels (Seabold and Perktold, 2010) for Python and MTS (Tsay, 2014) for R, both of which use one of the two MLE approaches above. Empirically these packages are slow even for moderate dimensions (d ≳ 10) and moderate sample sizes (T ≳ 104 ); a single fit can take hours and is prone to convergence failures, as discussed at length in Lütkepohl (2006). Section 7.1 documents these failures on simulated data. The remainder of this paper develops an alternative approach that scales well with both d and T .

3 Truncated VARMA models The classical VARMA(p, q) likelihood is expensive to evaluate (Section 2.2). We exploit instead the infinite-VAR representation (9) of a stationary invertible VARMA(p, q) and truncate it at a finite lag m, producing what we call a truncated VARMA(p, q, m) model : a VAR(m) whose coefficients are constrained to lie on a lower-dimensional structured subset of Rm×d×d defined by the recursions (10)-(11). The model is well-defined provided that the spectral radii of the AR and MA companion matrices stay strictly below 1 (Section 3.2), and choosing m large enough relative to the MA companion’s spectral radius makes the truncation arbitrarily accurate (Section 3.3). The associated loss functions and pre-computed sufficient statistics, derived in Sections 4–6, keep the parameter count and per-iteration cost independent of T . 3.1 The truncated VARMA(p, q, m) model A stationary invertible VARMA(p, q) admits the infinite VAR representation (9), Zt =

∞ X

Πj Z t−j + at ,

j=1

with coefficients Π1 , . . . , Πp contained in the (Π1:p , Θ1:q , Σ) parametrization, and (Πj )j>p generated by the recursion (11) from (Π1:p , Θ1:q ). It is natural to consider a truncated version of this infinite sum. Approximating an invertible process by a finite autoregression of slowly growing order is the classical longautoregression device of Hannan and Kavalieris (1984) and Lewis and Reinsel (1985). We adopt the same approximation but do not fit the m d2 coefficients of an unrestricted 9

Paulin and Elvira

VAR(m): the truncation enters only the loss, whose free parameters remain the (p + q) d2 VARMA coefficients. Definition 2 (Truncated VARMA(p, q, m)) Let p, q, m be positive integers with m ≥ p + q. The truncated VARMA(p, q, m) model parameterised by (Π1:p , Θ1:q , Σ) is Zt =

m X

Πj (Π1:p , Θ1:q ) Z t−j + at ,

(15)

j=1

where at ∈ Rd is a vector white-noise process with mean 0 and covariance matrix Σ ∈ Rd×d , and the lag matrices Πj ∈ Rd×d are obtained from the recursion (11) truncated at m terms. Equivalently, (Π1 , . . . , Πm ) lies in the image
of the polynomial map (Π1:p , Θ1:q ) 7→ Π1:m , m×d×d 2 a subset of R with at most (p + q)d degrees of freedom. Under Assumption 1 the spectral radius spr(CΘ ) of the MA companion in (3) is strictly j below 1, and Eq. (17) together P∞ with Gelfand’s formula (18) imply ∥Πj ∥ = O(ρ ) for any ρ ∈ (spr(CΘ ), 1). The tail j=m+1 Πj Z t−j representing the truncation error is therefore exponentially small in m. The truncated VARMA(p, q, m) admits a VAR(m) representation, and therefore has Markov order at most m, but is parameterized by only (p + q) d2 free coordinates rather than the m d2 of an unconstrained VAR(m). The constraint encodes the assumption that the VAR coefficients arise from an underlying VARMA(p, q) structure, which provides both interpretability and statistical efficiency, since the effective parameter count is independent of m. 3.2 Spectral radius constraints For the truncated VARMA(p, q, m) of Definition 2 to be a meaningful approximation of the underlying VARMA(p,q), the spectral radii of the AR and MA companion matrices must remain strictly below 1. This subsection makes Assumption 1 explicit as a constraint on (Φ1:p , Θ1:q ) and links it to the geometric decay of the infinite-VAR coefficients Πj . Throughout the rest of the paper we enforce these constraints by construction via the unconstrained PAR reparametrisation of Section 2.1, so no projection step is needed during the L-BFGS optimisation of the loss functions of Sections 4–5. Decay of Πj via the companion matrices. Recall from Section 2.1 that Assumption 1 is equivalent to spr(CΦ ) < 1 and spr(CΘ ) < 1 for the block companion matrices defined in Eq. (3). Inverting the recursion (10) for j = 1, . . . , p also yields the reverse identity min(j,q)

Φj = Πj −

X

Θk Πj−k ,

j = 1, . . . , p,

(16)

k=1

with Π0 = −I, which expresses the AR coefficients Φ1:p entering CΦ in terms of (Π1:p , Θ1:q ). This recovery costs O(pq d2 ) flops. The invertibility constraint spr(CΘ ) < 1 also governs the rate at which the truncation approximates the infinite-VAR representation (9): setting V Π := (Πp−q+1 , . . . , Πp )′ (with Π0 = −I and Πj = 0 for j < 0), the recursion (11) gives, for j ≥ p + 1,
′ j−p Πj = E ′q CΘ V Π, E q := 0, . . . , 0, I ∈ Rqd×d . (17) 10

Scalable estimation of VARMA models

Gelfand’s formula (Lax, 2002) states that for any submultiplicative matrix norm ∥ · ∥, spr(C) = lim ∥C k ∥1/k , k→∞

∥C k ∥1/k ≥ spr(C) for every k ≥ 1.

(18)

Applied to (17), it implies that ∥Πj ∥ decays geometrically: for every ρ > spr(CΘ ) there is a constant Cρ with ∥Πj ∥ ≤ Cρ ρj , so ∥Πj ∥ → 0 as j → ∞ (the precise, non-asymptotic constant is given in Appendix B). The truncation error of the VARMA(p, q, m) relative to the infinite-VAR representation therefore decays geometrically in m once both spectral radii are bounded away from 1; the next subsection makes this quantitative. 3.3 Truncation bias The truncated VARMA(p, q, m) generates a different stationary process than the original VARMA(p, q), because the truncation drops a non-zero (though geometrically small) tail of the infinite-VAR sum. We quantify the mismatch through the stationary autocovariances. ′ Let Γ∞ l := E[Z t Z t−l ] denote the autocovariance of the original VARMA(p, q) at lag l, and (m)

Γl the autocovariance of the truncated VARMA(p, q, m) generated by Eq. (15) with the same (Φ1:p , Θ1:q , Σ). Proposition 3 (Non-asymptotic geometric truncation bias) Suppose Assumption 1 holds, let ρ ∈ (spr(CΘ ), 1), write A(z) := Θ(z)−1 Φ(z) for the infinite-VAR operator, and set
µmin := min λmin A(z)∗ A(z) , Gρ := max I − A(z) op , |z|=1 |z|=1/ρ

√ (19) β• := inf σmin A(z) , β◦ := min σmin A(z) = µmin . |z|≤1

|z|=1

Under Assumption 1 all four are finite, with 0 < β• ≤ β◦ (Gρ measures the non-normality √ of the moving-average inversion; β◦ = µmin is the on-circle conditioning of the spectral density of Basu and Michailidis, 2015). Define n Gρ ρ m+1 β• o m0 := min m ≥ p + q : ≤ . 1−ρ 2

(20)

Then for every m ≥ m0 the truncated VARMA(p, q, m) is itself stationary and (m)

Γ∞ l − Γl

op

≤ Cρ ρm

for every l ≥ 0,

Cρ :=

6 ∥Σ∥op Gρ ρ 3/2

.

(21)

(1 − ρ) µmin

The constants Gρ , µmin , β• , β◦ , m0 are computable functions of (Φ, Θ, ρ) alone, and Cρ depends additionally on Σ; in particular the order-zero covariance Γ∞ 0 is recovered up to an explicit O(ρm ) error, uniformly in lag l. The bound (21) follows by comparing the spectral densities of the two stationary processes: they share the same P innovation covariance Σ and their autoregressive operators differ only by the dropped tail j>m Πj , whose operator norm decays geometrically by a Cauchy estimate on the analytic operator Θ(z)−1 Φ(z) (with an explicit, non-asymptotic constant), the passage to the autocovariance gap being governed by the on-circle conditioning µmin of the spectral density. The complete, fully non-asymptotic proof is given in Appendix B. 11

Paulin and Elvira

In practice the constants in (21) are unimportant: for any ρ < 1, choosing m = O(log(1/ε)/ log(1/ρ)) makes the bias smaller than any prescribed tolerance ε. The truncated VARMA(p, q, m) is therefore a faithful surrogate for the original VARMA whenever m is several times the time-scale set by the slowest mode of CΘ .

4 Least Squares loss for truncated VARMA models We now develop the loss functions used to fit the truncated VARMA(p, q, m) model of Definition 2 from data. Section 4.1 introduces the basic least-squares objective Lm and shows that it admits a sufficient-statistics decomposition that makes each evaluation cost O(m2 d3 ), independent of the number of observations T . For a pure VAR, the reduction of least-squares criteria to sample autocovariances is classical (see, e.g., Lütkepohl, 2005); the point here is that the same fixed-size summary, sufficient for the conditional Gaussian truncated-VAR regression criterion, continues to drive the nonlinear truncated-VARMA problem, in which Π1:m depends on the free parameters through the map (22). Section 4.2 then explains how to evaluate Lm in practice from these sufficient statistics. Section 4.3 adds Bayesian regularisation to obtain the regularised least-squares (RLS) loss. We work throughout in the unconstrained PAR parametrisation of Section 2.1: the free (p+q) d2 Θ parameters are (M Φ , with the VARMA coefficients obtained as 1:p , M 1:q ) ∈ R Θ (Φ1:p , Θ1:q ) = T (M Φ 1:p , M 1:q ), min(j,q)

Πj = Φj +

X

Θk Πj−k

Π0 = −I,

(j = 1, . . . , p),

(22)

k=1

followed by the recursion (11) to extend Π1:p to Π1:m . Because T is a smooth bijection 2 of R(p+q)d onto the open strictly-stable set (Section 2.1), the spectral-radius constraints of Section 3.2 are enforced by construction and no projection step is required during L-BFGS optimisation. The objectives Lm , LRLS m , and the marginalised MAP loss of Section 5 are all (p+q)d2 Θ minimised as unconstrained functions of (M Φ via the composite map 1:p , M 1:q ) ∈ R of (22). 4.1 Efficient Least Squares loss function of VAR Least squares loss functions have been widely used in time series forecasting, and their consistency for estimating VAR process parameters under fourth moment noise assumptions (i.e. not necessarily Gaussian noise) has been established in Mann and Wald (1943). In the setting of VAR processes with Gaussian noise, the least squares estimator is known to be equivalent to the MLE (see Sections 3.3 - 3.4 of Lütkepohl (2005)). We consider the one-step-ahead squared prediction error of our truncated VARMA Lm (Π1:p , Θ1:q ) =

T −1 X

m X

t=m

j=1

2

Πj z t+1−j − z t+1

,

(23)

where Π1:m are the first m coefficients of the infinite VAR representation (9) of the VARMA process, these can be computed recursively based on the parameters Π1:p , Θ1:q by (11) for 12

Scalable estimation of VARMA models

j ≥ p + 1. Note that the parameters Π1:p , Θ1:q allow us to recover the original parameter Φ1:p by using (10) for 1 ≤ j ≤ p. This loss function attempts to find the parameters of a VARMA process whose mtruncated infinite VAR representation has the smallest 1-step prediction error (in MSE). Note that unlike ℓ∞ , ℓm does not depend on the noise covariance parameter Σ (as this does not affect the coefficients of the infinite VAR representation). 4.2 Evaluation of the loss function First, we introduce some necessary notation. Let     zt Π1     z t:t−m+1 =  ...  , Π1:m =  ...  , z t−m+1

Πm

 Π′1   Π1:m =  ...  . Π′m 

Here z t:t−m+1 is a column vector of dimension md, while Π1:m and Π1:m are block matrices with md rows and d columns (m blocks of size d × d), with Π corresponding to block-wise transpose. Using these, and expanding the square with the cyclic property of the trace, the loss can be written as ′

′

Lm (Π1:m ) = S 0 − 2 tr(Π1:m S 1 ) + tr(Π1:m S 2 Π1:m )

(24)

for S 0 :=

T −1 X

∥z t+1 ∥2 ,

S 1 :=

t=m

T −1 X

z t:t−m+1 z ′t+1 ,

t=m

S 2 :=

T −1 X

z t:t−m+1 z ′t:t−m+1 .

t=m

After we have computed S 0 , S 1 and S 2 , further evaluation of the loss function no longer requires access to the observations, i.e. the cost no longer depends on the number of observations T . This is in contrast with the original log-likelihood function, where we need to access all observations for each likelihood evaluation. Note that Lm is a function of the truncated AR coefficients Π1:m alone (it does not take Θ as a separate argument). For a VARMA(p, q) the Π1:m are the truncation of its AR(∞) expansion, produced from (Φ1:p , Θ1:q ) by the recursion (22); in the unconstrained PAR coΘ ordinates this composes as (Φ1:p , Θ1:q ) = T (M Φ by Π1:m = Π1:m (Φ, Θ), 1:p , M 1:q ) followed
Θ so the objective actually minimised is Lm Π1:m (T (M Φ , M )) , which we abbreviate 1:p 1:q Φ Θ Lm (T (M , M )) in the sequel (and likewise R(T (·)) for the residual scatter). 4.3 Regularized Least Squares The least-squares loss Lm of Section 4.1 has no explicit regularization on the coefficient matrices Π1:p and Θ1:q . When the number of free parameters pd2 + qd2 is of the same order as the effective sample size N := T − m, the unregularized optimum can become highly unstable. In this subsection we describe a lightweight Bayesian regularization that preserves the sufficient-statistics computational structure of Lm while introducing data-adaptive shrinkage on the coefficient matrices. We call the resulting estimator the Regularized Least Squares (RLS) estimator. 13

Paulin and Elvira

Gaussian likelihood with isotropic noise. We adopt the simplest Gaussian model in which the innovation covariance is held fixed at the identity, Σ = I d . Under this assumption the conditional Gaussian negative log-likelihood of the truncated VAR(m) working model (up to a constant that does not depend on the parameters) reduces to
− log p z m+1:T | Π1:p , Θ1:q , Σ = I d = 21 Lm (Π1:m ) + const,

(25)

i.e., half the quadratic criterion (24). In particular, the likelihood depends on the observations only through the sufficient statistics S 0 , S 1 , S 2 of Section 4.1, so each evaluation is O(m2 d3 ), independent of T . Gaussian shrinkage prior on the PAR parameters. Following the approach of Heaps (2023), we place a centred Gaussian shrinkage prior directly on the unconstrained PAR Θ parameters (M Φ 1:p , M 1:q ), partitioning the entries into four groups (diagonal vs off-diagonal, AR vs MA side) and assigning each group a single scalar variance. For every j = 1, . . . , p, k = 1, . . . , q, and 1 ≤ a, b ≤ d with a ̸= b:
iid 2 2 (M Φ j )aa | σΦ,diag ∼ N 0, σΦ,diag ,
iid 2 2 (M Θ k )aa | σΘ,diag ∼ N 0, σΘ,diag ,

iid 2 2 (M Φ j )ab | σΦ,off ∼ N 0, σΦ,off ,
iid 2 2 (M Θ k )ab | σΘ,off ∼ N 0, σΘ,off .

2 2 2 2 control the amount of shrinkage , σΘ,off , σΘ,diag , σΦ,off The four scalar variances σΦ,diag applied to each group. Because T is locally the identity at M = 0 (the contraction map of Eq. (5) satisfies dP s /dM s = I there), the diagonal/off-diagonal split on M closely matches the corresponding split on the original (Φ, Θ): shrinking diagonal M -entries less aggressively allows the within-channel partial autocorrelations to remain near their datadriven values, while off-diagonal entries (cross-channel partial effects) are penalised more heavily. Separating AR from MA variances accommodates the different roles the two parts play in the Wold representation. Pushed forward by T , this Gaussian prior on M induces a smooth, fully supported prior on the strictly stable VARMA coefficients (Φ1:p , Θ1:q ) that increasingly penalizes approaches to the stability boundary in the unconstrained coordinates (Heaps, 2023). The result is a soft, stationarity-encouraging prior with no truncation discontinuity. The four group variances are treated as fixed hyperparameters. Rather than hand-tune them, we select them once per data set by cross-validation: a dyadic grid search chooses the variances minimizing an out-of-sample predictive log-likelihood on a held-out tail of the first training window. We defer the details of this selection procedure to the numerical 2 2 2 2 experiments (Section 7), and for now simply treat σΦ,diag , σΦ,off , σΘ,diag , σΘ,off as given constants.

A separate, fixed ridge on the first-lag own-channel diagonals. The within-channel Θ first-order dynamics (the diagonal entries (M Φ 1 )aa and (M 1 )aa of the lag-one AR and MA matrices) carry the dominant persistence of each series and are typically large and well Θ identified. By contrast, the higher-lag own-channel diagonals (M Φ j )aa and (M k )aa with j, k ≥ 2 decay quickly and are mostly noise. Folding all lags into a single diagonal variance is therefore undesirable: the many small higher-lag entries pull the shared variance down and over-shrink the lag-one persistence, which is exactly the signal we wish to preserve. 14

Scalable estimation of VARMA models

We therefore exclude the lag-one own-channel diagonals from the cross-validated diagonal 2 groups and give them their own centred Gaussian prior with a separate fixed variance σdiag,1 ,
iid 2 (M Θ 1 )aa ∼ N 0, σdiag,1 ,

iid 2 (M Φ 1 )aa ∼ N 0, σdiag,1 ,

a = 1, . . . , d,

(26)

2 chosen large (σdiag,1 = 4 in our experiments) so that these coefficients are only lightly ridged and stay close to their data-driven values, while the off-diagonal and higher-lag diagonal 2 groups can be shrunk more aggressively. Keeping σdiag,1 separate from the cross-validated groups is deliberate: a single variance shared with the near-zero higher-lag diagonals would couple the lag-one persistence to them, defeating the purpose. With this split, the crossvalidated diagonal groups below run over the remaining lags only (j = 2, . . . , p and k = 2, . . . , q), so their entry counts become (p − 1) d and (q − 1) d (and vanish entirely for first-order models p = q = 1); the off-diagonal groups are unchanged.

The regularized least-squares loss. Writing ndiag := d and noff := d(d − 1), and collecting the four cross-validated group sums of squares (the diagonal sums running over lags ≥ 2) together with the two lag-one diagonal sums held out for the fixed ridge, diag := SΦ

diag := SΘ

p X d X 2 (M Φ j )aa ,

off SΦ

p X X 2 (M Φ := j )ab , j=1 a̸=b

q X d X

q X X

k=2 a=1

off SΘ :=

d X 2 (M Φ 1 )aa , a=1

j=2 a=1 2 (M Θ k )aa ,

diag,1 := SΦ

2 (M Θ k )ab ,

diag,1 := SΘ

d X 2 (M Θ 1 )aa , a=1

k=1 a̸=b

the negative log of the Gaussian shrinkage prior is, up to additive constants, a sum of quadratic group-ridge penalties 21 Sg /σg2 . Adding these to the half-quadratic data term 1 2 Lm gives the regularized least-squares (RLS) loss diag,1 diag,1
SΦ + SΘ Φ Θ Φ Θ 1 LRLS (M , M ) = L T (M , M ) + m 1:p 1:q 1:p 1:q 2 m 2 2 σdiag,1 diag diag off off SΦ SΘ SΦ SΘ + + + + . 2 2 2 2 2 σΦ,diag 2 σΦ,off 2 σΘ,diag 2 σΘ,off

(27)

Θ Φ Θ The data-fit term Lm T (M Φ 1:p , M 1:q ) uses Eq. (22) to obtain Π1:m from (M 1:p , M 1:q ) and depends on the data only through the sufficient statistics S 0 , S 1 , S 2 ; the five ridge penalties depend on no data. The per-group variances are the fixed, cross-validated hyperparameters 2 introduced above (with the lag-one own-channel diagonals held at σdiag,1 ). The objective 2

(p+q)d Θ is a smooth function of the unconstrained (M Φ , minimized by L-BFGS 1:p , M 1:q ) ∈ R without any spectral-radius constraint.

Properties. Θ • Unconstrained parameter space. The state vector for L-BFGS is (M Φ 1:p , M 1:q ) ∈ 2 R(p+q)d , which is unconstrained (no spectral-radius constraints, no projections, no augmented-Lagrangian penalties). The group variances are fixed hyperparameters

15

Paulin and Elvira

chosen by cross-validation, not free parameters, so the optimization is over the coefficients alone, avoiding both the extra dimensions and the parameterization sensitivity of a joint MAP that retains log σ 2 as variables.
Θ • Sufficient-statistics evaluation. The likelihood term 21 Lm T (M Φ 1:p , M 1:q ) depends on the data only through S 0 , S 1 , S 2 (cf. Section 4.1); the ridge-prior terms depend on no data. The PAR map T adds an O((p2 +q 2 ) d3 ) overhead per evaluation (a Durbin–Levinson recursion, O(k 2 d3 ) per factor of order k). Each value-and-gradient evaluation of LRLS therefore costs O(m2 d3 + (p2 +q 2 ) d3 ), independent of T . m • Relation to ridge regression. The penalty in (27) is a group-wise Tikhonov (ridge) regularizer on the PAR coordinates. Specialising to Θ1:q ≡ 0 and a regime where T is well-approximated by the identity (small M Φ , i.e. near-white-noise) with a com1 2 2 mon variance σ 2 across groups, it reduces to 2σ1 2 ∥M Φ 1:p ∥F + const ≈ 2σ 2 ∥Φ1:p ∥F , i.e. Tikhonov-regularized (ridge) regression on the truncated VAR representation. Compared to the Σ-marginalized MAP estimator of Section 5, the RLS estimator is the cheaper of the two: it avoids the O(d3 ) Cholesky and log-determinant of Ψ0 + R at every iteration, and its likelihood is simply the quadratic Lm already computed from the sufficient statistics. The price paid is that Σ is assumed isotropic, an assumption that is accurate only after the observations z t have been approximately whitened, e.g., by a rolling standardization step as in our applications.

5 MAP estimation with marginalized covariance The least-squares loss Lm from the previous section treats all residual directions equally and does not model the noise covariance Σ. In this section, we extend the approach to a full Bayesian framework in which Σ is given an inverse-Wishart prior and then marginalized out analytically. The result is a marginalized MAP loss that (i) still depends on the data only through precomputed sufficient statistics, (ii) adapts to the noise structure without requiring explicit estimation of Σ, and (iii) stays numerically well-behaved (bounded below even when the residual scatter is rank-deficient) thanks to the positive-definite inverseWishart scale Ψ0 . Coefficient shrinkage is supplied separately, by the PAR-coordinate priors of Section 4.3, not by the covariance marginalization. 5.1 Gaussian likelihood and residual scatter matrix Let N := T − m and define the residual vectors rt := z t+1 −

m X

Πj z t+1−j ,

t = m, . . . , T − 1.

j=1

Conditionally on Σ, the conditional Gaussian negative log-likelihood of the truncated VAR(m) working model, given the first m observations (up to a constant), is − log p(z m+1:T | Π1:p , Θ1:q , Σ) = 16

N 1 log det(Σ) + tr Σ−1 R , 2 2

(28)

Scalable estimation of VARMA models

where R ∈ Rd×d is the residual scatter matrix R :=

T −1 X

rt r′t .

(29)

t=m

Crucially, R can be expressed in terms of the precomputed sufficient statistics from Section 4.1. Let T −1 X (2) S 0 := z t+1 z ′t+1 ∈ Rd×d t=m

be the matrix-valued analogue of the scalar statistic S 0 . Then expanding the definition of rt gives (2)

R = S0 −

m X j=1

Πj S 1 (j) −

m X

S 1 (j)′ Π′j +

j=1

m X m X

Πi S 2 (ij)Π′j .

(30)

i=1 j=1 (2)

Thus R, and hence the likelihood (28), depends on the observations only through S 0 , S 1 , and S 2 . 5.2 Inverse-Wishart prior on Σ and marginalization We place a conjugate inverse-Wishart prior on the noise covariance: Σ ∼ IW(Ψ0 , ν0 ), with scale matrix Ψ0 ∈ Rd×d and degrees of freedom ν0 . The prior is proper precisely when Ψ0 ≻ 0 and ν0 > d − 1; the improper (Jeffreys-like) limit Ψ0 = 0, ν0 = 0 is also admissible and is discussed below. In our implementation Ψ0 = ψ0 I d for a scalar ψ0 ≥ 0, and throughout this paper we use the default ψ0 = 1, ν0 = d + 4, i.e. Σ ∼ IW(I d , d + 4), which is proper (I d ≻ 0 and d + 4 > d − 1). The choice ν0 = d + 4 is the smallest integer degrees-of-freedom value for which the inverse-Wishart prior on Σ has finite second moments (the inverse-Wishart variance exists only for ν0 > d + 3), so the prior is only mildly informative while still letting the data dominate. The density of the inverse-Wishart distribution is proportional to   |Σ|−(ν0 +d+1)/2 exp − 21 tr(Ψ0 Σ−1 ) . Combining with the Gaussian likelihood (28) and integrating over Σ, we obtain the marginal likelihood Z p(z m+1:T | Π1:p , Θ1:q ) = p(z m+1:T | Π1:p , Θ1:q , Σ) p(Σ) dΣ ∝ |Ψ0 + R|−(ν0 +N )/2 . This is a standard result from the conjugacy of the inverse-Wishart with the multivariate Gaussian (see, e.g., Chapter 3 of Gelman et al., 2013). The negative log marginal likelihood is therefore − log p(z m+1:T | Π1:p , Θ1:q ) =

ν0 + N log det(Ψ0 + R) + const. 2 17

(31)

Paulin and Elvira

Remark 4 (Default prior and improper limit) Throughout, the default is the proper prior Σ ∼ IW(I d , d + 4), which adds the unit ridge I d to R inside the log-determinant in (31). Setting ν0 = 0, Ψ0 = 0 recovers the improper (Jeffreys-like) prior p(Σ) ∝ |Σ|−(d+1)/2 , under which the marginal likelihood reduces to N2 log det R + const; we do not use this limit, as the R-only criterion can degenerate when R is singular. 5.3 Global coefficient regularization We place the same Gaussian shrinkage prior on the unconstrained PAR parameters as for RLS (Section 4.3): the fixed lag-one own-channel diagonal ridge (26), and the four crossvalidated groups with the diagonal groups running over lags ≥ 2, inducing via push-forward through T a smooth, stability-respecting prior on the VARMA coefficients (Heaps, 2023). In terms of the group sums of squares of Section 4.3, the negative log-prior is, up to constants, RegΦ =

diag diag,1 off SΦ SΦ SΦ + + , 2 2 2 2 σΦ,diag 2 σΦ,off 2 σdiag,1

(32)

RegΘ =

diag diag,1 off SΘ SΘ SΘ + + , 2 2 2 2 σΘ,diag 2 σΘ,off 2 σdiag,1

(33)

identical to the coefficient penalties in the RLS loss (27). This global (lag-pooled) parameterisation uses one variance per group rather than one per lag, which empirically improves stability when p and q are large relative to N , at the cost of preventing lag-specific shrinkage. 5.4 The Σ-marginalized MAP loss Combining the marginalized likelihood (31) with the M-prior regularization (32)–(33), the full marginalized MAP loss is Φ Θ Lmarg m (M 1:p , M 1:q ) =

ν0 + N Θ log det Ψ0 + R(T (M Φ 1:p , M 1:q )) + RegΦ + RegΘ , 2

(34)

where R is the residual scatter matrix (30) evaluated at (Π1:m , Θ1:q ) obtained from Θ (M Φ 1:p , M 1:q ) via the composite map of (22), and N = T − m. The objective is min(p+q)d2 Θ imised as an unconstrained function of (M Φ ; the PAR map T enforces 1:p , M 1:q ) ∈ R spr(CΦ ) < 1 and spr(CΘ ) < 1 by construction, so no projection or augmented-Lagrangian step is needed. This loss has several desirable properties: 1. Covariance-free. The noise covariance Σ has been marginalized out analytically and does not appear as a parameter. This reduces the parameter count by d(d + 1)/2 compared to approaches that optimize over Σ directly. 2. Sufficient-statistics form. As with RLS, the loss depends on the data only through (2) S 0 , S 1 , S 2 ; after the one-time pre-pass, each evaluation costs O(m2 d3 + (p2 +q 2 )d3 ), independent of T . 3. Numerical stability from Ψ0 . The positive-definite scale Ψ0 keeps the data term bounded below: since Ψ0 +R ⪰ Ψ0 ≻ 0, we have log det(Ψ0 +R) ≥ log det Ψ0 > −∞, 18

Scalable estimation of VARMA models

so even a rank-deficient residual scatter R (as when the truncated model interpolates the in-sample residuals) cannot drive the criterion to −∞. The inverse-Wishart prior therefore regularizes the covariance estimate, shrinking it toward Ψ0 ; shrinkage of the coefficient matrices Π1:p , Θ1:q comes separately, from the PAR-coordinate priors of Section 4.3. 4. Unconstrained optimisation. Stationarity and invertibility are guaranteed by T , with no constraints on the state vector. Remark 5 (Recovering Σ̂) Although Σ is not optimized over, a point estimate can be recovered after fitting. The posterior mode of Σ under the inverse-Wishart prior, given Θ Φ Θ
Φ the MAP estimates (M̂ 1:p , M̂ 1:q ) and the corresponding R̂ = R T (M̂ 1:p , M̂ 1:q ) , is Σ̂ = (Ψ0 + R̂)/(ν0 + N + d + 1). Optimisation by L-BFGS. The RLS loss (27) and the marginalised MAP loss (34) are (p+q)d2 Θ both smooth functions of (M Φ , composed of matrix multiplications, 1:p , M 1:q ) ∈ R FFTs, a Cholesky log-determinant, and the fixed-variance ridge penalties; a value-andgradient pair is obtained in one reverse-mode autodiff pass (Bradbury et al., 2018), storing no recursion intermediates. The PAR map T enforces spr(CΦ ) < 1 and spr(CΘ ) < 1 by 2 construction, so both losses are minimised as unconstrained functions on R(p+q)d . Because the parameter dimension (p + q)d2 grows quadratically with d and easily reaches the thousands, we use the limited-memory quasi-Newton method L-BFGS (Liu and Nocedal, 1989), which builds an implicit positive-definite Hessian approximation from the last few value-gradient pairs and produces each search direction in time linear in the parameter dimension. L-BFGS also warm-starts naturally across rolling-window refits, the previous window’s optimum serving as the initial point when the window shifts by small steps (Appendix F); the experiments of Section 7 instead refit each window from a pilot estimate (Section 7). We optimize this nonconvex objective numerically and do not claim convergence to the global optimum. The recovery guarantee of Section 5.6 is stated for a global empirical minimizer over a compact ball (and, more generally, for any point whose objective is near the minimum, Theorem 6 and Remark 8); it therefore bounds the statistical error of a well-optimized fit but does not certify that the L-BFGS iterate reaches it. 5.5 Multiplicative seasonal VARMA Many series, hourly or daily measurements in particular, carry a strong periodic component: even after a deterministic seasonal mean is removed, the anomalies retain a season-recurring autocorrelation that an ordinary VAR or dense VARMA can only mimic with many extra lags. The multivariate analogue of the multiplicative seasonal ARMA of Box et al. (2015, Ch. 9), the seasonal VARMA(p, q) × (P, Q)s , captures this parsimoniously by multiplying a regular and a seasonal lag polynomial on each side: 

I−

P X j=1

(s) Φj B sj



I−

p X

i





Φi B Z t = I −

i=1

Q X k=1

(s) Θk B sk



I−

q X

 Θl B l at ,

(35)

l=1

where B is the backshift operator (BZ t = Z t−1 ) and s is the seasonal period. Expanding the two products gives an ordinary VARMA whose AR and MA polynomials have orders 19

Paulin and Elvira

p + sP and q + sQ, so (35) is a restricted VARMA: of the p + sP AR lags only those at the regular lags 1:p, the seasonal lags s, 2s, . . . , sP , and the interaction lags sr + j (1 ≤ r ≤ P , 1 ≤ j ≤ p; possibly overlapping the regular and seasonal lags) carry nonzero coefficients, and these are determined by just (p + P )d2 free parameters (and likewise (q + Q)d2 on the MA side) rather than the (p + sP + q + sQ)d2 of the unrestricted expansion. A single seasonal factor at lag s thus encodes a full period of dynamics that an unrestricted VAR would need ≥ s lags and s d2 parameters to represent. The model (35) plugs directly into the truncated-VARMA machinery of Sections 3–5: the four factor polynomials are convolved once into the expanded (Φ1:p+sP , Θ1:q+sQ ), the infinite-VAR coefficients Πj that solve the Box–Jenkins relation (10)–(11) are formed directly from these factors by a single FFT inversion, truncated at j = m, and fed to the same sufficient-statistic loss. Stationarity and invertibility are again enforced coordinatefree: each of the four factors is parameterized in its own unconstrained PAR coordinates (s) (s) through the map τ of Section 2.1, so writing (M Φ , M Φ , M Θ , M Θ ) for the four unconstrained stacks, the maps τp , τP , τq , τQ make each factor’s block-companion matrix Schur stable, so each factor polynomial has no determinantal zeros in the closed unit disk. The determinant of the expanded AR (or MA) operator is the product of its factors’ determinants, evaluating the seasonal factor at z s (with |z s | ≤ 1 whenever |z| ≤ 1), hence also nonvanishing there, so (35) is stationary and invertible by construction, and the Σ-marginalized MAP loss (34) (or, with Σ = I, the RLS criterion of Section 4.3) is minimized as an unconstrained function of the (p + P + q + Q)d2 factor parameters. We apply this seasonal VARMA in the meteorological and air-quality benchmarks of Sections 7.3 and 7.4. Relation to the standard multiplicative seasonal VARMA. Equation (35) is the multiplicative seasonal VARMA of Lütkepohl (2005); Tsay (2014) (see also Reinsel, 2003), written in the convention in which the seasonal factor left-multiplies the regular factor: in the scalar case the ordering of the factors is immaterial, but matrix coefficients need not commute, so for d > 1 the convention fixes the restricted family, and we use the stated one throughout. Our treatment changes how the model is parameterized and fit, not the model class. The PAR map acts per factor, and the four factor polynomials are then multiplied  in coefficient space to form the expanded (Φ1:p+sP , Θ1:q+sQ ); since det Φ(z) Φ(s) (z s ) = det Φ(z) det Φ(s) (z s ) (likewise on the MA side), the expanded operator is stable and invertible exactly when each factor is, so parameterizing each factor over its own stable region via τ traverses precisely this restricted family. The departures from the textbook treatment are in estimation only: the per-factor PAR-shrinkage and inverseWishart priors (Sections 4.3–5.3), unconstrained optimization in the PAR coordinates, and the truncated infinite-VAR loss, whose approximation improves as m → ∞. Computational cost. The per-gradient cost of fitting (35) (independent of T and, at leading order, identical to that of a plain VARMA of the same dimension) is analyzed in detail in Appendix G, building on the Fourier evaluation developed in Section 6. Seasonal VARMAX. The exogenous-regressor (VARMAX) extension composes with the seasonal factorization, giving a seasonal VARMAX : the side of (35) is augPright-hand sx k acting on observed mented with a finite exogenous lag polynomial Ξ(B) = Ξ B k=0 k 20

Scalable estimation of VARMA models

covariates xt ∈ Rdx (a known price or temperature path, say): 

I−

P X j=1

(s)

Φj B sj



I−

p X

Q q     X X (s) Φi B i Z t = Ξ(B) xt + I − Θk B sk I − Θl B l at .

i=1

k=1

l=1

P Solving for the innovation gives the truncated residual at = Z t − m j=1 Πj Z t−j − Pm x Π x , in which Π are exactly those of the seasonal VARMA and Πx is 1:m j t−j j=0  (s) s −1 the truncation of Θ (z ) Θ(z) Ξ(z), computed by the same FFT inversion. The sufficient-statistic structure, the T -independent gradients (working dimension d + dx ), the incremental rolling updates, and the built-in stationarity and invertibility all carry over, the exogenous block Ξ being unconstrained; Appendix H gives the details. The case P = Q = 0 is the plain VARMAX of the retail benchmark (Section 7.2); the seasonal case suits periodic series driven by known covariates, forecast conditional on the covariate path. 5.6 Statistical guarantees We close the development of the two point estimators with a statistical guarantee. Its target is the sequence of infinite-VAR coefficients Π1:∞ of Eq. (9), equivalently, the one-step predictor, or the transfer function Θ(z)−1 Φ(z). Unlike the pair (Φ1:p , Θ1:q ), this functional is identifiable (Section 2.4): observationally equivalent VARMA parametrisations share the same Π1:∞ , and it is Π1:m alone that enters our losses and forecasts. Accordingly, in the spirit of the classical transfer-function consistency theory for Gaussian VARMA maximum likelihood (Dunsmuir and Hannan, 1976; Hannan and Deistler, 2012), we prove recovery of Π1:∞ without any identifiability conditions, and (as is standard in that literature) over a compact parameter set, for exact or approximate empirical minimizers. Assumption 2 (i) The observations z 1:T are a sample of a stationary d-dimensional Gaussian VARMA(p̄, q̄) process (1) whose parameters (Φ01:p̄ , Θ01:q̄ ) satisfy Assumption 1, with i.i.d. innovations at ∼ N (0, Σ), Σ ≻ 0. (ii) The fitted orders satisfy p ≥ p̄ and q ≥ q̄. Writing npar := (p + q)d2 and ϑ for npar Θ the stacked PAR coordinates (M Φ , let ϑ0 := T −1 (Φ01:p , Θ01:q ) be 1:p , M 1:q ) ∈ R the PAR coordinates of the true coefficients padded with zero matrices up to orders (p, q) (padding leaves the lag polynomials, hence Assumption 1, unchanged and appends only zero eigenvalues to the block companion spectrum, and T is a bijection onto the stable set by Proposition 1 and the padding observation of Appendix A, so ϑ0 is well defined). Assume ϑ0 ∈ KB := {ϑ ∈ Rnpar : ∥ϑ∥2 ≤ B} for a fixed radius B. Minimizing over KBpneeds no projection step in practice: the smooth radial squashing σB (ϑ) := Bϑ/ B 2 + ∥ϑ∥22 , the scalar analogue of the contraction M 7→ M (I + M ⊤ M )−1/2 inside T itself, is a diffeomorphism of Rnpar onto the open ball, so unconstrained optimization in the squashed coordinates produces KB -valued points whose objective comes arbitrarily close to the minimum over KB ; these are the approximate minimizers that Theorem 6 admits through its tolerance η, although a minimizer lying on the boundary sphere itself is not attained at any finite squashed coordinate. 21

Paulin and Elvira

(iii) ϑ̂T is a measurable KB -valued estimator targeting the RLS loss (27) or the Σmarginalized MAP loss (34) (with Ψ0 ≻ 0 and ν0 > d − 1 in the latter case, so the inverse-Wishart prior is proper). The hyperparameters (Ψ0 , ν0 ) are fixed; the group ridge variances may be data-dependent, selected in any measurable way (for instance by validation on a held-out block of the same trajectory, as in Section 7), provided 2 each selected variance is at least a fixed floor σmin > 0. Writing L for the objective at the selected variances, the estimator is required to satisfy, for a (possibly random) η ≥ 0, L(ϑ̂T ) ≤ inf L(ϑ) + η, ϑ∈KB

and we call the smallest such η its optimization error (η = 0 for an exact minimizer over KB ). Theorem 6 (Recovery of Π1:∞ for approximate empirical minimizers) Let Assumption 2 hold, write Πj (ϑ̂T ) for the infinite-VAR coefficients (10)–(11) of the fitted parameters and Π0j := Πj (ϑ0 ) for those of the truth, and set Λδ := npar log T + log(1/δ). There exist constants ρ∗ ∈ (0, 1), C ≥ 1, c > 0, m1 ∈ N and T0 , depending only on (Φ0 , Θ0 , Σ), the orders (p, q, d), the radius B and the fixed hyperparameters (Ψ0 , ν0 ) and 2 the variance floor σmin (but not on T , m or δ) such that for every T ≥ T0 , every δ ∈ (0, 1/2) with Λδ ≤ c T , and every truncation lag m with m1 ≤ m ≤ T /2, any estimator satisfying Assumption 2(iii), for either loss, obeys with probability at least 1 − δ, ∞ Λ X η 2 δ Πj (ϑ̂T ) − Π0j F ≤ C + ρ∗2m + . (36) T T j=1

No lower bound on m beyond the constant m1 is imposed: the truncation enters only through the reported term ρ2m ∗ , and each corollary below imposes of it what it needs; the choice m ≥ log T / log(1/ρ∗ ) makes that term at most T −2 . Corollary 7 (Consistency and rate in summed operator norm) Under Assumption 2: P P (i) (Consistency.) If m = m(T ) → ∞ with m ≤ T /2 and η/T → 0, then j≥1 ∥Πj (ϑ̂T )− P

Π0j ∥op −→ 0. (ii) (Rate.) If in addition m ≥ log T / log(1/ρ∗ ) and there is an event ET with P(ET ) ≥ 1 − δ ′ on which η ≤ C ′ npar log T for a constant C ′ (in particular, deterministically with δ ′ = 0), then with probability at least 1 − δ − δ ′ , r q   X Λδ log T npar 0 Πj (ϑ̂T ) − Πj op ≤ C = O log T for fixed δ. T T j≥1

22

Scalable estimation of VARMA models

The proof is given in Appendix C. The coupling between m and T is weak: the truncation may be any lag in [m1 , T /2] (the only remaining requirement, Λδ ≤ cT , is a sample-size condition that does not involve m), and the rate regime of part (ii) already holds at m ≍ log T , the same logarithmic growth of the fitted order as in the long-autoregression theory of Hannan and Kavalieris (1984) and Lewis and Reinsel (1985). This reflects the geometric truncation bias of Proposition 3: a bias decaying only polynomially in m would force a polynomial balance between m and T . In practice m is set from the time-scale of the slowest moving-average mode of the fitted model, which is well inside the feasible set and decays at a fixed rate; the bound above is stated through the uniform rate ρ∗ over the whole set, and is correspondingly conservative. The T −1/2 dependence of the rate cannot be improved: already in a one-parameter Gaussian VAR(1) √ submodel, a two-point (Le Cam) argument shows that no estimator does better than c/ T for this loss (Remark 13); we do not pursue a matching npar -dependent lower bound for the summed operator-norm loss. Both estimators are handled by one decomposition of the empirical excess loss, into a self-normalized martingale and a Gaussian quadratic form; the Σ-marginalized JMAP loss localizes to a Σ−1 -weighted version of the same quadratic form, so the two share the rate and differ only in constants, as the experiments of Section 7 bear out, RLS and JMAP separating only through the conditioning of Σ. Gaussianity is used only through the innovations’ moment-generating function and concentration of Gaussian quadratic forms (we do not pursue the classical fourth-moment extensions). Finally, the theorem is a fixed-order, fixeddimension result: its constants depend on (p, q, d) and B, and regimes in which d, p or q grow with T are not treated; the high-dimensional claims of this paper are computational (cf. Basu and Michailidis, 2015 for sparse-VAR theory in the growing-dimension regime). Remark 8 (Optimization error and compactness) Because Theorem 6 is stated for approximate empirical minimizers, it already covers inexact optimization: the guarantee degrades only through the additive η/T , so an exact minimizer (η = 0) attains the rate, and by Corollary 7 so does any ϑ̂ whose objective is within O(npar log T ) of the minimum with high probability. Assumption 2(iii) also covers the cross-validated ridge variances of the experimental pipeline, since the proof touches the selected variances only through the penalty bound PB and the selection respects a fixed floor (the validation search box of Section 7 is bounded below). The implementation minimizes the same objectives by unconstrained LBFGS, which returns a local stationary point of a nonconvex objective; we do not certify that its optimization error η is small, so the theorem informs but does not by itself guarantee the computed estimate. For the statistical guarantee we analyze the compactly constrained version of the estimator, minimizing over the fixed ball KB of Assumption 2(ii), which contains a representative ϑ0 of the true transfer function. Constraining to a compact set is standard for nonconvex, stability-constrained time-series parametrizations (Dunsmuir and Hannan, 1976; Hannan and Deistler, 2012): on KB the spectral radii of all fitted companion matrices are bounded away from one, ensuring the uniform geometric-decay and Lipschitz constants of the implied infinite-VAR coefficients on which the proof rests (Lemma 11 in Appendix C), and, more generally, uniform mixing of the associated processes. The implementation minimizes the same objectives by L-BFGS that is unconstrained in the PAR coordinates, stability and invertibility being supplied by the map itself (Proposition 1), with the ridge penalties keeping the iterates in a bounded region in practice; the ball is a proof 23

Paulin and Elvira

device only and is not enforced, and the radius B may be taken as large as desired, at the price of larger constants in Theorem 6. These constants are not uniform in B, and they deteriorate as the fitted models approach the stability boundary: a larger B admits companion matrices with spectral radius closer to one, so the uniform decay rate ρ∗ of Lemma 11 tends to one while the geometry constants Ḡ, L̄ grow, and the implicit constants degrade for processes whose roots sit near the unit circle. The guarantee extends to the multiplicative seasonal model of Section 5.5 to the seasonal model with only notational changes, because the estimators enter the analysis solely through the expanded infinite-VAR coefficients and the smooth, per-factor PAR parametrization. Corollary 9 (Seasonal VARMA) Let the observations be a sample of a stationary Gaussian multiplicative seasonal VARMA (35) with period s and factor orders (p̄, q̄, P̄ , Q̄), each of whose four factor polynomials is zero-free on the closed unit disk, with i.i.d. innovations N (0, Σ), Σ ≻ 0. Fit the seasonal model with the same period, factor orders p ≥ p̄, q ≥ q̄, P ≥ P̄ , Q ≥ Q̄, and truncation m, by minimizing the RLS or Σ-marginalized MAP loss (the latter with Ψ0 ≻ 0 and ν0 > d − 1) over a ball KB in the four unconstrained PAR stacks (s) (s) (M Φ , M Φ , M Θ , M Θ ) ∈ Rnpar , where now npar := (p + q + P + Q) d2 and KB contains the inverse-PAR coordinates of the true four factor polynomials after zero-padding their coefficient arrays to the fitted orders. Then the conclusions of Theorem 6 and Corollary 7 hold for the expanded infinite-VAR coefficients Π1:∞ , with constants depending in addition on the period s. Moreover, ρ∗ may be taken of the form ρ̃ 1/s with ρ̃ < 1, so the rate regime of Corollary 7(ii) requires m ≳ s log T : the truncation must resolve the seasonal time-scale, consistent with the settings m ≫ s used in Sections 7.3 and 7.4. The short proof (stability and analyticity of the expanded operators on a disk of radius (1/ρ̄)1/s , after which the argument of Appendix C applies verbatim) is given in Appendix C.6.

6 Efficient likelihood evaluation via Fourier methods Both losses developed in the previous two sections, the regularised least-squares loss LRLS m of Section 4.3 and the Σ-marginalised MAP loss Lmarg of Section 5, are built around the m truncated least-squares criterion Lm of Section 4.1. The RLS loss is 21 Lm plus four dataindependent ridge penalties; the MAP loss is ν0 +N log det(Ψ0 + R) plus the same ridge 2 penalties, where R is the residual scatter matrix (30) that produces Lm = tr(R). In both cases the entire data-dependence of the loss reduces to the same block-quadratic form in the truncated VAR coefficients, Q(Π1:m ) :=

m X m X

Πi S 2 (i, j) Π′j ∈ Rd×d ,

(37)

i=1 j=1

P −1 where the m × m block S 2 has (i, j) block S 2 (i, j) = Tt=m z t+1−i z ′t+1−j (cf. Eq. (24)). Materialised as a 4-tensor S 2 ∈ Rm×m×d×d , this object would require m2 d2 floats of storage and Θ(m2 d3 ) flops per evaluation of Q, both of which are prohibitive in the large-m regime 24

Scalable estimation of VARMA models

needed for the truncation (21) to be accurate, since L-BFGS makes many evaluations per fit and we refit on a rolling basis. We avoid materialising S 2 entirely by exploiting its block-Toeplitz-plus-correction structure. Each block depends almost exclusively on the lag τ = i − j: b i−j − ∆head − ∆tail , S 2 (i, j) = C i,j i,j

(38)

PT −|τ |

z s+max(−τ,0) z ′s+max(τ,0) is the (untruncated) empirical autocovariance P −τ b −τ = C b′ . z s z ′s+τ , every index lies in [1, T ], and C of Z at lag τ ; for τ ≥ 0 this reads Ts=1 τ tail remove the terms of C b i−j whose summation index falls outside The corrections ∆head , ∆ i,j i,j the (narrower) range of S 2 (i, j). By the block symmetry S 2 (i, j) = S 2 (j, i)′ it suffices to record them for i ≥ j (so τ = i − j ≥ 0); blocks with i < j then P follow by transposition, −i which sidesteps any negative-lag indexing. For i ≥ j, S 2 (i, j) = Ts=m−i+1 z s z ′s+τ sums b τ , and the difference is over [m − i + 1, T − i] rather than the full range [1, T − τ ] of C exactly b τ := where C

s=1

∆head = i,j

m−i X

z s z ′s+τ ,

s=1

∆tail i,j =

TX −i+j

z s z ′s+τ

(τ = i − j ≥ 0),

(39)

s=T −i+1

whose indices all lie in [1, T ] (the ranges [1, m−i], [m−i+1, T −i], [T −i+1, T −τ ] partition [1, T − τ ]) and which involve only the boundary segments Y head := (z 1 , . . . , z m−1 )′ ∈ R(m−1)×d and Y tail := (z T −m+1 , . . . , z T )′ ∈ Rm×d . The decomposition (38) replaces the m2 d2 tensor S 2 by three lighter objects: b 0:m−1 ∈ Rm×d×d (md2 floats); • the autocovariance sequence C • the boundary blocks Y head ∈ R(m−1)×d and Y tail ∈ Rm×d (∼ 2md floats); total O(md2 ) storage, a factor-m reduction. These objects are computed once per training window in a single FFT-based pre-pass over the data, costing O(T log m d2 ) (Appendix G), and are stored alongside S 0 , S 1 as part of the sufficient statistics of Section 4.1. Toeplitz part via Parseval. Frequency-domain evaluation of Gaussian time-series criteria goes back to the Whittle likelihood (Whittle, 1953); the identity below applies the same Parseval mechanism not to the data likelihood, whose frequency sum grows with T , but to the quadratic form in the truncated coefficients, whose length is set by m. Substituting (38) into (37) splits Q into a block-Toeplitz quadratic plus two boundary Grams: Q(Π1:m ) =

m X

b i−j Π′ − (Π ∗ Y head )⊤ (Π ∗ Y head ) − (Π ⋆ Y tail )⊤ (Π ⋆ Y tail ), Πi C j | {z } {z } | i,j=1 head Gram tail Gram | {z } Toeplitz part Qtoep

(40) where ∗ denotes linear convolution and ⋆ a cross-correlation, both of length-m filters Π with the boundary segments. The Toeplitz part is the Gram of the truncated residuals:P writing P d for the length-m filtering of Z by Π, one has Qtoep (Π) = ′ r t := m Π z ∈ R k t−k k=1 t rt rt , 25

Paulin and Elvira

which makes its d × d shape and positive-semidefiniteness manifest. By the convolution ef z ef = Π ef , so Parseval’s theorem applied to this matrix Gram gives theorem r toep

Q

Nf −1 1 X e e e∗ (Π) = Πf C f Πf , Nf

(41)

f =0

e f is the DFT e f = Pm Πk e−2πif k/Nf is the DFT of the coefficient sequence, C where Π k=1 b of the Hermitian-extended autocovariance C at Nf = 2m frequencies (a real, positivesemidefinite spectral-density estimate), and ∗ denotes conjugate transpose. The orientation e f on the left and Π e ∗ on the right) is the one consistent with the block convention (Π f b i−j Π′ of (40). Eq. (41) is computed in three steps: Πi C j e ← FFTN (Π): O(m log m · d2 ); 1. Π f ff ← Π ef C e f for each f : O(m d3 ); 2. M PNf −1 f e ∗
3. Qtoep ← N1f f =0 ℜ M f Πf : O(m d3 ), a direct frequency-domain contraction (the real part discards the imaginary terms, which cancel by the Hermitian symmetry of the summand). e depends on the data only, not on Π, and is cached once per window. The spectral density C The zero-padding to Nf = 2m makes the length-Nf circular convolution reproduce the linear filtering exactly, with no aliasing. Boundary parts via short convolutions. The head and tail Grams in (40) are quadratic forms in length-m convolutions of Π with the boundary blocks of size O(m × d). Each convolution is a single FFT-based filter at length Nf = 2m, costing O(m log m · d2 ); the Gram matrix is then an O(md2 ) outer-product reduction. Both boundary corrections are therefore cheaper than the Toeplitz part. Total cost. Summing the three contributions in (40), the per-evaluation cost of Q (and hence of R in (30) and of Lm in (24)) is O(m log m · d2 ) + | {z } FFTs

O(m d3 ) | {z }

.

(42)

per-frequency d×d matmuls and final Π⊤ H contraction

Both the storage and the compute scale only linearly in m (up to logarithmic factors), so the Fourier construction makes long-memory regimes tractable. The savings persist down the autodiff graph: reverse-mode differentiation through Eq. (41) is again a length-Nf FFT pair plus a per-frequency conjugate matmul, so the gradient costs the same order as the forward. Remarks. • The same Parseval trick handles the residual scatter R in Eq. (30) verbatim: the P 3 linear term j Πj S 1 (j) is already O(md ) without help, but the quadratic term P ′ i,j Πi S 2 (i, j)Πj is exactly Q(Π1:m ). 26

Scalable estimation of VARMA models

• The single O(T log m d2 ) pre-pass is amortised across all L-BFGS iterations of a single fit, and across all daily refits in a rolling-window scheme via incremental updates of b τ , Y head and Y tail . the cumulative sums underlying C The multiplicative-seasonal factors of Section 5.5 and the exogenous regressors of the VARMAX extension leave this evaluation unchanged at leading order: a seasonal VARMA evaluates at the same O(m log m d2 + m d3 ) per-gradient cost as a plain VARMA of equal dimension, and the VARMAX at the same cost with d replaced by d + dx , independent of T in all cases. The full per-term accounting (61), including the dependence on the orders (p, q, P, Q) and the period s, is given in Appendix G.

7 Numerical results We evaluate the two pointwise estimators of Sections 4.3–5: the regularized least-squares estimator (RLS) and the Σ-marginalized MAP estimator (JMAP). Both use the PARshrinkage prior of Sections 4.3–5.3 with per-group ridge variances chosen by cross-validation 2 = 4 of Eq. (26). We report forecast (Section 7) and the fixed lag-one diagonal ridge σdiag,1 mean-squared error (MSE) by forecast horizon. Measured timings for both are reported in Appendix D. Common cross-validation protocol. On all four benchmarks we tune every method’s shrinkage and order hyperparameters by a single, common protocol, so the model classes compete on an equal footing. We hold out a validation tail of the first training window, one fifth of its length in every experiment, fit each candidate on the observations preceding it, and score its out-of-sample predictive log-likelihood over the validation origins; the minimizing hyperparameters are then frozen for the test. For the PAR-shrinkage VARMA(X) the prior variances are selected on a dyadic grid. The own-channel diagonal groups are held at σdiag = 2, with the lag-one own-channel 2 = 4 of Eq. (26), and one shared value is swept over the diagonals at the fixed ridge σdiag,1 1 0 −1 eight points σoff ∈ {2 , 2 , 2 , . . . , 2−6 }; this value governs the off-diagonal Φ and Θ groups together with the exogenous groups where a VARMAX is fit. The eight fits share no state, so the selection runs one process per grid point; an earlier version of these experiments used a derivative-free Nelder–Mead search over the four group variances, whose simplex is sequential and needed about sixty objective evaluations per candidate order. Each fit, in the variance search and the final refit alike, runs L-BFGS for 2000 iterations at limited-memory depth 100 on the sufficient-statistics loss of Sections 4.3–5, initialized from a Hannan– Rissanen pilot (Hannan and Rissanen, 1982) fit on the training observations and mapped to the PAR coordinates; the selected configuration is refit on the full training window under the same budget. The discrete order (p, q[, s]) is swept by an outer grid, the variance search run within each candidate order. The ridge and Bayesian VAR baselines jointly select the lag order p and their prior variances, as detailed below; the componentwise baselines select their orders by the same criterion, per product on the retail panel and shared across the series on the hourly panels. VAR baselines. The ridge VAR(X) fits each equation by generalized ridge least squares 2 , cross-series coefficients by on the lag design, penalizing own-lag coefficients by 1/σdiag 27

Paulin and Elvira

2 , and, where price enters, the exogenous block by its own variance; multi-step fore1/σoff casts are iterated from the fitted coefficients, conditioning on the known future price where applicable. The Bayesian VAR(X) places the corresponding diagonal/off-diagonal Gaussian prior on the coefficients and treats Σ as a plug-in, the residual covariance refined by a few posterior-mean/residual-covariance fixed-point iterations; conditional on that Σ the coefficient posterior is Gaussian in closed form. Forecasts are posterior-predictive means, averaged over 120 posterior coefficient draws at every horizon: on the retail panel the hstep forecast is iterated per draw and the forecasts averaged, and for the pure VAR the same average is computed by averaging the per-draw iterated coefficient maps, which is equivalent because the forecast is linear in those maps. At h = 1 the forecast is also linear in the coefficients, so the average reduces to a posterior-mean plug-in; at longer horizons it is the genuine model average. On the retail panel the diagonal variance is pinned and the cross-series and exogenous variances are tied to the swept dyadic value; on the hourly panels the diagonal and off-diagonal variances are swept on their own grids. Orders and prior variances are selected by the same predictive-log-likelihood criterion as every other model.

Feasible set. Stability and invertibility hold at every iterate by construction, through the PAR map of Section 2.1, so the optimization is unconstrained and needs no projection step. Two devices from the theory are implemented on top of the map: the coefficients are rescaled as Φj 7→ c j Φj (and likewise for Θ) with the fixed cap c = 0.995 (substituting λ = cµ in the characteristic polynomial shows this scales the companion spectrum by exactly c), and the ball KB of Assumption 2(ii) is enforced through the squashing σB of that assumption with radius B = 100. Both enter through the same map for the two MAP fits, so the estimators explore one parameter set. On the synthetic and retail benchmarks the fits are interior and the devices are inert: the fitted PAR norms stay below 23, where σB is a rescaling by 0.97 or more, and the cap is slack. On the hourly panels they are not inert, and this is their designed role: the diurnal cycle is close to deterministic, the fitted seasonal factors push toward the unit root, and the PAR coordinates grow without bound as the companion radius approaches one, so several JMAP fits reach norms near B with the capped radius attaining c. There the ball and the cap supply the persistence bound that the compactness assumption calls for, which is exactly the job the theory gives them; the forecasts are insensitive to the residual drift in these coordinates, which runs along weakly identified AR–MA cancellations. The recovery guarantee is unaffected, being stated for any KB -valued near-minimizer, which the squashed optimizer produces by construction. The optimization error remains uncertified (Remark 8), and it, rather than the feasible-set devices, is what separates the computed estimate from the analysed one. Significance tests. Forecast comparisons are Diebold–Mariano tests on the per-origin squared-error loss, with a Newey–West long-run variance and the Harvey et al. (1997) small-sample correction against a tn−1 reference. Every comparison uses the JMAP fit of the panel’s leading model class as a fixed reference rather than the per-horizon winner, so a column reports one estimator against all others; the reference row itself carries no stars. Stars give the two-sided level of that comparison, ∗ p < 0.1, ∗∗ p < 0.05, ∗∗∗ p < 0.01. On the hourly panels n runs to 3.3 × 104 (Beijing) and 5.3 × 104 (Singapore) origins, so the test 28

Scalable estimation of VARMA models

resolves differences of a tenth of a percent and a star there records detectability rather than practical size, for which the MSE values themselves should be read. 7.1 Synthetic data We simulate stationary and invertible VARMA(2, 2) processes directly in the unconstrained Θ PAR coordinates of Section 2.1: each unconstrained matrix M Φ j , M k (j, k = 1, 2) has independent Gaussian entries with standard deviation sdiag on the diagonal and soff off the diagonal, and is mapped through T to a strictly stable (Φ1:2 , Θ1:2 ). The innovations are Gaussian, at ∼ N (0, Σ), with a random, strongly ill-conditioned covariance Σ = Q diag(λ1 , . . . , λd ) Q⊤ , where Q is a Haar-random orthogonal matrix and the eigenvaliid

ues are log-uniform, log λi ∼ Unif(−10, 0), giving a condition number up to e10 ≈ 2 × 104 (median ≈ 5 × 103 over replications); each series is afterwards standardized to unit training variance. This anisotropic, cross-correlated noise is deliberate: it is what makes the covariance model matter, and hence what separates the Σ-aware JMAP from the isotropic-noise RLS. We consider two regimes: a dense regime (sdiag , soff ) = (0.5, 0.5) and a sparse regime (0.5, 0.2) with weaker cross-channel coupling. For each regime we sweep the dimension d ∈ {10, 20, 40}, and let the training length grow with it: T ∈ {200, 400, 800} at d = 10, {400, 800, 1600} at d = 20, and {800, 1600, 3200} at d = 40. A VARMA(2, 2) carries 4d2 coefficients against a scalar sample of size N d, so the information ratio N/(4d) is held near 4.5 on the shortest rung of every dimension and doubles along each ladder; a grid fixed in T instead confounds dimension with sample size, and its d = 40, T = 200 corner supplies fewer observations than parameters. Each cell is averaged over 8 independent replications. All estimators are fit with truncation m = 20; the per-group ridge variances are selected by the cross-validation of Section 7, and JMAP places the inverse-Wishart prior IW(I d , d + 4) on the innovation covariance. Both fits run L-BFGS for 2000 iterations. Because the datagenerating process is known, we compute the exact minimum-MSE VARMA forecast and report each estimator’s error relative to this oracle, MSE/MSEoracle ≥ 1, so that 1 is optimal. Tables 2 and 3 report the oracle-relative error at h = 1, 2, and Figure 1 plots the h = 1 value against T . The Σ-aware JMAP leads everywhere. It improves on the isotropic-noise RLS in all eighteen cells and in all 144 cell-seed pairs, by 20.7% on average at h = 1 (paired t = −17.8 across the pooled seeds), the gain coming from modelling the cross-channel covariance RLS ignores. The advantage is a finite-sample one: it shrinks monotonically along every ladder, from 30.9% to 12.6% at d = 10 and from 27.5% to 9.9% at d = 40, as both estimators converge on the oracle. The two estimators also select different shrinkage: along the dense d = 10 ladder RLS relaxes from σoff = 0.25 to 1.5 while JMAP stays at 0.25; at the larger dimensions both relax with T . The oracle-relative error still grows with d and decays toward 1 with T , and the sparse regime sits closer to its oracle than the dense one, its weak cross-channel block carrying little recoverable signal. The sparse-VARMA package bigtime trails throughout, by 16 to 41% in the dense regime and 11 to 24% in the sparse one, the gap widening with d and barely closing as T grows. Two standard R VARMA packages serve as external reference points on the same data. bigtime (Wilms et al., 2023), a sparse VARMA fit by penalised least squares, is reliable and appears as a column of Tables 2–3: it trails the Σ-aware JMAP in all eighteen cells of the design, and beats the isotropic RLS in nine of them, concentrated at the shortest training 29

Paulin and Elvira

Table 2: Synthetic data: oracle-relative forecast MSE at h = 1 (lower is better; 1 is the oracle). The training length grows with the dimension, holding the information ratio N/(4d) near 4.5 on the shortest rung of each d. Best in each row in bold. bigtime (Wilms et al., 2023) is the sparse-VARMA package baseline (cross-validated, same 8 seeds). dense (0.5, 0.5)

sparse (0.5, 0.2)

d

T

RLS

JMAP

bigtime

RLS

JMAP

bigtime

10

200 400 800

1.45 1.32 1.15

1.11 1.05 1.02

1.33 1.23 1.19

1.41 1.28 1.16

1.10 1.04 1.02

1.26 1.17 1.13

20

400 800 1600

1.48 1.32 1.15

1.09 1.04 1.03

1.42 1.33 1.30

1.37 1.31 1.16

1.07 1.03 1.02

1.29 1.22 1.19

40

800 1600 3200

1.41 1.22 1.14

1.11 1.06 1.04

1.55 1.47 1.43

1.34 1.16 1.12

1.07 1.05 1.03

1.33 1.26 1.24

Table 3: Synthetic data: oracle-relative forecast MSE at h = 2 (columns as in Table 2). dense (0.5, 0.5)

sparse (0.5, 0.2)

d

T

RLS

JMAP

bigtime

RLS

JMAP

bigtime

10

200 400 800

1.32 1.21 1.11

1.08 1.04 1.02

1.19 1.12 1.10

1.30 1.20 1.11

1.07 1.03 1.02

1.16 1.10 1.07

20

400 800 1600

1.36 1.25 1.11

1.07 1.04 1.03

1.23 1.17 1.15

1.29 1.22 1.12

1.06 1.02 1.02

1.18 1.13 1.10

40

800 1600 3200

1.32 1.17 1.11

1.09 1.06 1.03

1.26 1.22 1.19

1.25 1.13 1.09

1.05 1.04 1.03

1.17 1.13 1.12

length of each dimension and in the sparse regime, where its sparsity assumption is closer to the truth. Its gap to JMAP widens with the dimension, from 0.17 at d = 10, T = 800 to 0.39 at d = 40, T = 3200 in the dense regime. The other, MTS (Tsay, 2014), implements the classical Gaussian conditional MLE with no constraint on the fitted operators. We fit VARMA(2, 2) with it on the d = 10, T = 800 series of both regimes, the same sixteen series the estimators above are fit on. Two of the sixteen stop with a singular Hessian. Of the fourteen that return a fit, twelve have a non-invertible moving-average operator, with companion spectral radius between 1.016 and 1.444; every fitted autoregressive operator is stable, with radius at most 0.907, so it is invertibility alone that fails. A non-invertible Θ makes the innovation recursion that generates the forecasts grow geometrically in its own past output, and the forecasts follow: across those twelve runs the oracle-relative error at h = 1 ranges from 9.8×1023 to numerical overflow, against 1.02 for JMAP on the same series. The two invertible fits forecast at 1.10 30

Scalable estimation of VARMA models

200

400 T

800

1.3 1.25 1.2 1.15 1.1 1.05 1

RLS JMAP bigtime

200

400 T

800

MSE / MSEoracle

MSE / MSEoracle

sparse (0.5, 0.2), d = 10

1.3 1.25 1.2 1.15 1.1 1.05 1

dense (0.5, 0.5), d = 20 RLS JMAP bigtime

400

800 T

MSE / MSEoracle

1.4 1.35 1.3 1.25 1.2 1.15 1.1 1.05 1

1600

sparse (0.5, 0.2), d = 20 RLS JMAP bigtime

400

800 T

1600

MSE / MSEoracle

RLS JMAP bigtime

MSE / MSEoracle

MSE / MSEoracle

dense (0.5, 0.5), d = 10

1.35 1.3 1.25 1.2 1.15 1.1 1.05 1

1.35 1.3 1.25 1.2 1.15 1.1 1.05 1

dense (0.5, 0.5), d = 40 RLS JMAP bigtime

800

1600 T

3200

sparse (0.5, 0.2), d = 40

1.25 1.2 1.15

RLS JMAP bigtime

1.1 1.05 1

800

1600 T

3200

Figure 1: Oracle-relative forecast MSE, horizon-averaged over h = 1, 2, 3, versus training length T (both axes log), for the dense (0.5, 0.5) regime (top) and the sparse (0.5, 0.2) regime (bottom) at d = 10, 20, 40. The training length grows with the dimension, so each panel has its own ladder. The Σ-aware JMAP approaches the oracle (dotted line at 1) from above; RLS and bigtime trail it at every (d, T ), and bigtime’s gap widens with d. Error bars are ±1 standard error over the eight seeds.

and 1.15, still behind JMAP. We therefore exclude MTS from the tables. This is exactly the failure the PAR reparametrisation removes, the map T rendering every iterate strictly stable and invertible by construction. 7.2 Dominick’s retail scanner data: demand with price as an exogenous driver Our second real-data benchmark turns the retail panel into a VARMAX problem: forecasting category-level demand (units) while conditioning on the contemporaneous per-unit price, which embeds the store’s promotional activity. From the Dominick’s Finer Foods scanner database (Kilts Center, University of Chicago Booth) we forecast a d = 12 basket comprising every food and drink category of the archive with under 10% of weeks missing (cereal, cookies, bottled juice, cheese, canned soup, snack crackers, crackers, front-end candies, frozen juice, canned tuna, soft drinks, and frozen entrees) over the 399-week record; every excluded food or drink category is missing at least 23.6%, and the members share the shopping-basket and promotional co-movement a multivariate model can exploit. The endogenous series are log units, the exogenous series the log volume-weighted price per unit (total dollars over total units), so a promotion enters the model as a transient price cut. We model the first difference of the logs, which renders the unit-root-like levels stationary, and truncate the implied autoregressive representation at m = 20 lags. Forecasts use a rolling 31

Paulin and Elvira

Table 4: Dominick’s (d = 12 grocery basket, demand with price exogenous, first-differenced growth, truncation m = 20): standardized-scale forecast MSE on the 96-origin rolling test, by horizon. Every model is fit once on the first 299-week window and frozen for all origins; order and shrinkage are chosen by the cross-validation of Section 7 (selected order in parentheses). Because a small fraction of category-weeks are missing, each entry is the per-origin mean squared (standardized-growth) error over the categories observed at that origin, averaged equally over the 96 origins. Best per column in bold; stars give the Diebold– Mariano significance of the difference from the JMAP fit (Harvey et al. (1997)-corrected, Newey–West HAC long-run variance, t95 reference): ∗ p < 0.1, ∗∗ p < 0.05, ∗∗∗ p < 0.01. model

h=1

h=2

h=3

PAR-shrinkage VARMAX(p, q, s) with price exogenous, order and ridge variances by CV: RLS (1, 1, 1) 0.4442∗∗∗ 0.6938 0.7213 JMAP (1, 1, 1) 0.5051 0.6583 0.6789 PAR-shrinkage VARMA(p, q) on units only (price-free ablation), order and ridge by CV: RLS (1, 1) 0.7516∗∗∗ 1.0678∗∗∗ 1.0857∗∗∗ ∗∗∗ ∗∗∗ JMAP (1, 1) 0.7619 1.0719 1.0912∗∗∗ Ridge VARX(p, p), order and shrinkage by CV: VARX (p=1) 0.6541∗

0.8757∗

0.9268∗

Bayesian BVARX(p, p), order and prior by CV, model-averaged: BVARX (p=1) 0.6343∗ 0.8257∗

0.8650∗

Componentwise ARMAX (per-product SARIMAX), orders (p, q) and own-price lag r by CV: comp. ARMAX 0.5999∗ 0.7919∗ 0.7938 Hierarchical-lasso sparse VARX (bigtime), lag and penalty by CV: sparse VARX (p=4) 0.5361 0.7511

0.7855

Naive baselines: training-set mean stay-in-place (RW)

1.0848∗∗∗ 1.0848∗∗∗

1.0735∗∗∗ 1.0735∗∗∗

1.0743∗∗∗ 1.0742∗∗∗

test of 96 weekly origins at h = 1, 2, 3: every model is fit once on the first 299-week window and frozen, with the price for each forecast week supplied. This is therefore conditional forecasting given the observed (realized) price path, not an unconditional demand forecast: conditioning isolates the demand–price dynamics, and a planning application would substitute the retailer’s posted price for the realized one. Every method is tuned by the common protocol of Section 7, jointly over structure and shrinkage (selected orders in Table 4); the JMAP variant marginalizes Σ under the default IW(I d , d + 4) prior. The truncation is m = 20 and the validation tail is one fifth of the 299-week training window. The predictive log-likelihood used for selection here conditions on the known future price and is evaluated only on the channels observed at each origin, so the exogenous VARMAX models and the price-free ablations are scored on one criterion. Two contrasts are built in. To quantify the exogenous price we add the price-free VARMA counterparts (units alone); to separate the moving-average term from the cross-series coupling we add a ridge VARX and Bayesian BVARX (no MA), a per32

Scalable estimation of VARMA models

product componentwise ARMAX (SARIMAX, no cross-series coupling), and the hierarchicallasso sparse VARX of bigtime, the closest match the package offers, as it implements VARX but not VARMAX. Naive random-walk and training-set-mean baselines anchor the scale. Price dominates: every price-conditioned model roughly halves the naive h=1 error, and dropping price (the units-only VARMA) inflates the best h=1 MSE by two thirds. The moving-average term is also decisive, the full VARMAX beating the moving-average-free ridge and Bayesian VARX by 25–37% across horizons. This is over-differencing (firstdifferencing a near-integrated series induces a near-unit moving-average root), so the baselines that recover it otherwise, the longer sparse autoregression of bigtime (0.536 at h = 1, lag 4) and the componentwise ARMAX (0.600), fall between the VARMAX and the short VARX; as on the seasonal panels, the gain comes from the moving-average term, not the cross-series coupling. Among the VARMAX estimators, both cross-validated to order (1, 1, 1), the isotropic-noise RLS attains the lowest h=1 error, 12% below JMAP (the only significant difference between the two, p < 0.001), while JMAP is ahead by 5–6% at h = 2, 3; over horizons the two average within 1%. The remaining models trail with the significance shown in Table 4 (DM on the per-origin loss against the JMAP fit as the fixed reference, Newey–West HAC, Harvey et al. (1997) correction, n = 96); the sparse VARX is the closest and is not separated from JMAP. 7.3 Singapore meteorological data: a seasonal VARMA Our third real-data benchmark isolates a structure absent from the previous two: strong periodicity. Hourly weather is dominated by a diurnal cycle, and even after its deterministic part is removed the anomalies retain a daily-recurring autocorrelation. This is the regime for which the multiplicative seasonal VARMA is designed, and it lets us add a fourth model to the comparison, a seasonal VARMA whose seasonal factor an ordinary VAR or dense VARMA can only mimic with many extra lags. Model and prior. We fit the multiplicative seasonal VARMA of Section 5.5, Eq. (35), with diurnal period s = 24 and the four factor orders tied to a common p = P = q = Q chosen by cross-validation; stationarity and invertibility hold by construction through the per-factor PAR maps. All four factor-stacks share the same PAR-shrinkage prior, fixed lagone diagonal ridge, and per-group ridge variances as the regular VARMA, and (for JMAP) the same IW(I d , d + 4) prior on Σ of Sections 4.3–5. JMAP minimizes the Σ-marginalized MAP loss, RLS the isotropic-noise criterion. Data and protocol. We forecast d = 11 hourly meteorological series for Singapore (2 m temperature, relative humidity, dew point, precipitation, mean-sea-level pressure, total cloud cover, the zonal and meridional 10 m wind components, wind gusts, shortwave radiation, and wind speed) from the Open-Meteo reanalysis archive (ECMWF-IFS) over 2017– 2025. Each series is deseasonalized by regressing out, per channel by ordinary least squares on the training window only, an additive calendar design—an intercept, 24 hour-of-day indicators, a piecewise-linear annual cycle (a tent basis interpolating between monthly knots), and a piecewise-linear weekly cycle—and standardized by its training-window standard deviation, so the targets are standardized hourly anomalies on which a zero (climatologicalanomaly) forecast scores MSE ≈ 1. We use a rolling scheme of six annual test windows 33

Paulin and Elvira

Table 5: Singapore weather (d = 11 hourly anomalies): forecast MSE per origin (the mean over the 11 series), averaged equally over the test origins of the six annual windows (2020– 2025), at horizons h = 1, . . . , 6 hours. Order and ridge/prior hyperparameters are chosen once by the common cross-validation protocol on the first training window and held fixed across windows (selected order in parentheses). Best in each column in bold. model

h=1

h=2

h=3

h=4

h=5

h=6

Seasonal PAR-shrinkage VARMA(p, p, p, p)24 , p and ridge variances by CV: RLS (p=1) 0.2768∗∗∗ 0.4621∗∗∗ 0.5762∗∗∗ 0.6535∗∗∗ 0.7075∗∗∗ JMAP (p=4) 0.2747 0.4535 0.5605 0.6306 0.6777

0.7473∗∗∗ 0.7111

Dense VARMA(p, p) (no seasonal factor), order and ridge variances by CV: JMAP 0.2797∗∗∗ 0.4614∗∗∗ 0.5690∗∗∗ 0.6383∗∗∗ 0.6842∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ RLS 0.2800 0.4609 0.5683 0.6376 0.6837∗∗∗

0.7161∗∗∗ 0.7160∗∗∗

Seasonal VAR(p, 0, P, 0)24 (no MA), order by CV: RLS 0.2783∗∗∗ 0.4640∗∗∗ 0.5791∗∗∗ ∗∗∗ ∗∗∗ JMAP 0.2786 0.4649 0.5800∗∗∗

0.6571∗∗∗ 0.6576∗∗∗

0.7114∗∗∗ 0.7112∗∗∗

0.7509∗∗∗ 0.7501∗∗∗

Ridge / Bayesian VAR(p), order and prior by CV: ridge VAR 0.2796∗∗∗ 0.4609∗∗∗ 0.5678∗∗∗ ∗∗∗ ∗∗∗ Bayesian VAR 0.2780 0.4584 0.5649∗∗∗

0.6364∗∗∗ 0.6330∗∗

0.6815∗∗∗ 0.6780

0.7128 0.7092

Componentwise univariate baselines, order by CV: AR 0.2904∗∗∗ 0.4703∗∗∗ 0.5779∗∗∗ ARMA 0.2885∗∗∗ 0.4643∗∗∗ 0.5679∗∗∗

0.6478∗∗∗ 0.6349∗∗∗

0.6944∗∗∗ 0.6796

0.7267∗∗∗ 0.7109

Sparse VARMA package bigtime (Wilms et al., 2023), lag and penalty by CV: sparse VAR 0.3032∗∗∗ 0.5082∗∗∗ 0.6462∗∗∗ 0.7458∗∗∗ 0.8190∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ sparse VARMA 0.2954 0.5045 0.6454 0.7453 0.8179∗∗∗

0.8743∗∗∗ 0.8720∗∗∗

Naive baselines: training-set mean stay-in-place (RW)

1.1441∗∗∗ 1.2148∗∗∗

1.1441∗∗∗ 0.3388∗∗∗

1.1441∗∗∗ 0.6138∗∗∗

1.1441∗∗∗ 0.8161∗∗∗

1.1441∗∗∗ 0.9791∗∗∗

1.1441∗∗∗ 1.1094∗∗∗

(2020–2025); for each, every model is refit on the preceding three years (≈ 26,000 hourly observations) and scored at horizons h = 1, . . . , 6 hours. The seasonal period is s = 24 (diurnal) and the truncation m = 1024. Every hyperparameter is chosen once by the common protocol of Section 7, on a validation tail (the last 219 days, one fifth) of the first, 2017–2019 training window, then frozen and refit on each test year. The orders swept are p = P = q = Q ∈ {1, . . . , 4} for the seasonal VARMA and the seasonal-VAR ablation, p ∈ {12, 24, 36} for the nonseasonal dense-VARMA ablation, and p ∈ {24, 48, 96, 168} (AR) and p = q ∈ {24, 48} (ARMA), shared across the series, for the componentwise baselines; the ridge/Bayesian VAR (p ∈ {24, 48, 96}) and the sparse bigtime VAR/VARMA are tuned as on the other panels. The truncation is m = 1024 and the seasonal period s = 24. Table 5 reports the per-origin forecast MSE, averaged over the six test years. Table 5 exhibits a third regime, distinct from the previous two panels: the predictability is dominated by periodicity, and the seasonal VARMA captures it most economically. The seasonal PAR-shrinkage JMAP leads at the first five horizons, ahead of the seasonal RLS 34

Scalable estimation of VARMA models

by 0.76% at h = 1 with the gap widening to 5.1% by h = 6; at h = 6 the Bayesian VAR edges ahead of it by 0.3%. The separations are small in absolute terms but consistent across horizons. The advantage is one of parameter economy: cross-validation selects a seasonal order p = P = q = Q = 4 (just 4 · 4 d2 = 1936 coefficients, encoding an expanded VARMA of AR and MA order 4 + 24 · 4 = 100), whereas without the seasonal factor the dense VARMA climbs to p = 12–36 and the VAR to p = 48, three or more times as many coefficients, to reach across a full diurnal period. A single seasonal factor at lag 24 represents, with a third or less of the parameters, the daily-recurring autocorrelation that the flat models can only approximate with two dozen or more lags, and it also forecasts slightly more accurately, improving on the dense VARMA and the cross-validated VAR by 1.8% at h = 1. The sparse bigtime VAR/VARMA are the weakest multivariate models, behind even the componentwise baselines; sparsity is the wrong bias here, as the denser ozone panel of Section 7.4 examines in detail. Both the seasonal autoregression and the seasonal moving average contribute, and (unlike the differenced retail panel of Section 7.2, where a non-seasonal moving average carried the gain) here the gain comes from the seasonal structure, as on the Beijing panel below. Dropping the seasonal moving-average factor (the seasonal VAR) costs 1.3% at h = 1 and, because that factor reaches a full day ahead, a growing ≈5.5% by h = 6, whereas the ordinary moving average is inert (the dense VARMA tracks the ridge VAR to within half a percent at every horizon). The componentwise baselines, which discard all cross-channel coupling, are weakest throughout, so the coupling helps. A Diebold–Mariano test (Diebold and Mariano, 1995; Harvey et al., 1997) on the pooled hourly path (n = 52,608) makes every between-family gap significant at p < 0.01; the two seasonal PAR-shrinkage estimators are a practical tie, separable only by the large sample. 7.4 Beijing air-quality data: structure versus sparsity Our final benchmark revisits the periodic regime of Section 7.3 on a larger, denser panel, and uses it to ask a sharper question: when the predictable structure is spread diffusely across many channels and lags, does sparse estimation help or hurt? We forecast hourly ozone (O3 ) at the d = 12 monitoring stations of the Beijing multi-site air-quality network (Zhang et al., 2017) over 2013–2019 (stopping before the COVID-19 disruption). Ozone is photochemically driven, with a strong, regular diurnal cycle and pronounced cross-station coupling, so (as in the meteorological panel) its anomalies retain a daily-recurring autocorrelation that the multiplicative seasonal VARMA of Section 5.5 is built to capture. What is new here is the cross-station signal, which is high-dimensional, making this a natural test of the sparse VAR/VARMA estimators of the bigtime package (Wilms et al., 2023) against our dense, shrinkage-based estimators. Data and protocol. Each station’s hourly O3 series is preprocessed causally and independently within each window: a missing station-hour is filled by the contemporaneous cross-sectional mean of the other stations; the series are log-transformed; deseasonalized by regressing out, per station on the training window only, an additive calendar design (intercept, hour-of-day indicators, piecewise-linear annual and weekly cycles) together with a causal linear time trend that removes the well-documented multi-year decline in Beijing pollution; and standardized by the training-window standard deviation. We use four rolling 35

Paulin and Elvira

Table 6: Beijing ozone (d = 12 hourly station anomalies): forecast MSE per origin (the mean over the 12 stations), averaged equally over the test origins of the four rolling windows (test periods 2016–2019), at horizons h = 1, . . . , 6 hours. Order and ridge/prior hyperparameters are chosen by the common cross-validation protocol of Section 7 once on the first training window and held fixed across the four windows. Best in each column in bold. model

h=1

h=5

h=6

Seasonal PAR-shrinkage VARMA(p, p, p, p)24 , p and ridge variances by CV: RLS (p=2) 0.1903 0.3410 0.4584 0.5529 0.6302 JMAP (p=3) 0.1901 0.3412 0.4590 0.5539 0.6311

0.6917 0.6925

Dense VARMA(p, p) (no seasonal factor), order and ridge variances by CV: JMAP 0.1953∗∗∗ 0.3572∗∗∗ 0.4871∗∗∗ 0.5935∗∗∗ 0.6807∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ RLS 0.1951 0.3499 0.4699 0.5660 0.6447∗∗∗

0.7501∗∗∗ 0.7069∗∗∗

Seasonal VAR(p, 0, P, 0)24 (no MA), order by CV: RLS 0.1953∗∗∗ 0.3549∗∗∗ 0.4812∗∗∗ ∗∗∗ ∗∗∗ JMAP 0.1953 0.3550 0.4816∗∗∗

0.5834∗∗∗ 0.5840∗∗∗

0.6669∗∗∗ 0.6676∗∗∗

0.7334∗∗∗ 0.7342∗∗∗

Ridge / Bayesian VAR(p), order and prior by CV: ridge VAR 0.1951∗∗∗ 0.3543∗∗∗ 0.4790∗∗∗ ∗∗∗ ∗∗∗ Bayesian VAR 0.1955 0.3547 0.4794∗∗∗

0.5799∗∗∗ 0.5803∗∗∗

0.6626∗∗∗ 0.6629∗∗∗

0.7281∗∗∗ 0.7284∗∗∗

Componentwise univariate baselines, order by CV: AR 0.2197∗∗∗ 0.4130∗∗∗ 0.5548∗∗∗ ARMA 0.2163∗∗∗ 0.4014∗∗∗ 0.5349∗∗∗

0.6602∗∗∗ 0.6336∗∗∗

0.7384∗∗∗ 0.7067∗∗∗

0.7956∗∗∗ 0.7605∗∗∗

Sparse VARMA package bigtime (Wilms et al., 2023), lag and penalty by CV: sparse VAR 0.2577∗∗∗ 0.4027∗∗∗ 0.5304∗∗∗ 0.6361∗∗∗ 0.7211∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ ∗∗∗ sparse VARMA 0.2508 0.4061 0.5333 0.6385 0.7228∗∗∗

0.7876∗∗∗ 0.7886∗∗∗

Naive baselines: training-set mean stay-in-place (RW)

1.0354∗∗∗ 1.1551∗∗∗

1.0354∗∗∗ 0.2357∗∗∗

h=2

h=3

1.0354∗∗∗ 0.4764∗∗∗

1.0354∗∗∗ 0.6836∗∗∗

h=4

1.0354∗∗∗ 0.8644∗∗∗

1.0354∗∗∗ 1.0214∗∗∗

windows, each three training years and one test year (the last ten months, as the record ends December 2019), aligned to the March start of the data; the preprocessing is redone per window, while order and shrinkage hyperparameters are chosen once on the first window by the common protocol of Section 7 and frozen thereafter, as on the other panels. The seasonal period is s = 24, the truncation m = 1024, the covariance prior IW(I d , d + 4). To isolate the seasonal moving-average term we add a seasonal VAR ablation (q = Q = 0) whose order is chosen over the balanced grid (p, P ) = (k, k), by the same protocol as the seasonal VARMA. For bigtime, whose internal cross-validation is prohibitively slow here, we instead fit its coefficient path at each candidate lag (a single fit returns the whole penalty path) and select lag and penalty on the validation tail; its sparse-VAR and sparse-VARMA modes are scored exactly as our own models. The candidate lags are p ∈ {24, 48, 96} for the sparse VAR and (p, q) ∈ {(24, 2), (48, 2)} for the sparse VARMA, each fit on the cross-validation window and on the full training window, and the VAR(∞) inversion of the sparse VARMA is truncated at 200 lags before the forecast recursion. Table 6 reports the per-origin forecast MSE, averaged over the four windows. 36

Scalable estimation of VARMA models

The seasonal VARMA is again the most accurate model. The two seasonal PARshrinkage estimators lead at every horizon, tied to within 0.2% (JMAP marginally lower at h = 1, RLS thereafter). The remaining models cluster just above: at h = 1 the seasonal VARMA beats every non-seasonal alternative (dense VARMA, seasonal VAR, ridge VAR, all nearly coincident) by 2.6–2.8%, and the margin widens with horizon (to ∼5% over the VAR family and ∼6% over the seasonal VAR by h = 6), because the seasonal factor reaches a full day ahead by construction. At h = 1 neither the seasonal VAR (drop the MA) nor the dense VARMA (drop the seasonal factor) improves on the plain ridge VAR, so neither ingredient alone helps: it is the seasonal moving-average factor, acting with the seasonal autoregression, that delivers the gain. As on the meteorological panel, the moving average is inert only in its non-seasonal form (the dense VARMA-RLS tracks the ridge VAR, sparse VAR ≈ sparse VARMA), while the seasonal MA is not. The gain again comes cheaply: the seasonal VARMA reaches this accuracy with on the order of 8d2 coefficients, whereas a pure autoregression must span a full diurnal cycle to approximate what the seasonal factor encodes in a single low-order term. Averaged over horizons and windows, the seasonal VARMA leads, the dense VARMA-RLS and the VAR family follow nearly coincident, and the componentwise and sparse models trail; the componentwise baselines, which discard cross-station coupling, are again among the weakest, confirming the coupling is informative. The panel also speaks to sparsity. The bigtime sparse VAR and sparse VARMA are the least accurate multivariate models, and at h = 1 they fall behind the stay-in-place baseline as well (∼32–36% behind the seasonal VARMA at h = 1, and worse than every dense model at every horizon), and their own cross-validation explains why: on every window it drives the penalty to the least-sparse end of the path and selects the shortest lag offered, fighting the sparsity the method imposes and never reaching the diurnal lags the dense VAR exploits. Ozone’s spatiotemporal dependence is diffuse, many small cross-station and multi-lag couplings rather than a few large ones, so an ℓ1 /hierarchical penalty that zeroes coefficients discards usable signal, whereas the dense ℓ2 shrinkage of our estimators (and of the ridge VAR) keeps every coefficient, merely pulled toward zero. Sparsity is the wrong inductive bias here. Two qualifications are in order. First, the comparison mixes two differences: bigtime fits the moving-average part by penalized two-stage least squares rather than by a likelihood-based criterion, so part of its gap to the MAP estimators may reflect statistical efficiency rather than the sparsity penalty itself; the synthetic experiments of Section 7.1 point the same way, with bigtime trailing the Σ-aware estimators even in the weakly-coupled regime once Σ is estimable. Second, the conclusion is about these data, not about sparse methods in general: for data generated by a VARMA with sparse coefficients and a limited sample, the ℓ1 approach would likely recover the coefficients, and forecast, more accurately than our dense estimators. With the meteorological panel, this benchmark indicates that on real, periodically-structured, densely-coupled data the most accurate and most parameter-economical forecaster is a structured dense seasonal VARMA rather than a sparse one, stably across four independent pre-COVID years. A Diebold–Mariano test (Diebold and Mariano, 1995; Harvey et al., 1997) on the pooled path (n = 32,720) confirms every gap: the best seasonal VARMA leads its closest non-seasonal competitor (the dense VARMA-RLS) by 2.1–2.6% across horizons, the seasonal VAR by up to 6.0%, the componentwise baselines by 10–21%, and the sparse models by 14–36%, all at p < 0.01 (the two seasonal PAR-shrinkage estimators a practical tie). 37

Paulin and Elvira

Conclusion We have presented a scalable framework for estimating dense, truncated VARMA(p, q, m) models. It combines a partial-autocorrelation reparametrization, under which stationarity and invertibility hold by construction so that fitting is unconstrained; regularized leastsquares and Σ-marginalized MAP losses that see the data only through fixed-size sufficient statistics, making every gradient evaluation independent of the series length T ; and a Fourier/Parseval evaluation that brings the per-gradient cost down to near-linear in the truncation length m. Within this framework we fit two estimators (an isotropic-noise regularized least-squares point estimate and the Σ-marginalized MAP) and extend both, at similar computational cost, to multiplicative seasonal dynamics, to exogenous regressors (VARMAX), and to incremental rolling refits. Empirically the estimators stay close to the oracle on data where the classical conditional-MLE implementation returns non-invertible fits, while remaining tractable well beyond the reach of existing packages, and on retaildemand, meteorological, and air-quality benchmarks they are competitive with or better than VAR, Bayesian-VAR, component-wise ARMA, and sparse-VARMA baselines. In these examples the moving-average and seasonal moving-average terms account for much of the advantage, capturing with a single low-order factor the over-differencing or diurnal structure that a pure autoregression can only approximate with many extra lags; and where the predictable signal is diffuse rather than sparse, the dense structured VARMA is both more accurate and more parameter-economical than its sparse counterpart. A natural further step is full Bayesian inference, averaging forecasts over the Σ-marginalized posterior rather than reporting its mode. Sampling this posterior efficiently is a challenging open problem (the marginalization leaves a heavy-tailed, high-dimensional density), and we leave it for future research. Code availability. The implementation, together with drivers that reproduce every experiment in the paper, is available at https://github.com/paulindani/VARMA. The code is Python/JAX throughout, with an optional R dependency for the bigtime baseline; each experiment runs from a single script, and a master driver reruns the full study.

Acknowledgments and Disclosure of Funding DP is supported by a Nanyang Technological University Start-up Grant, project number: 024968-00001.

References Yasin Abbasi-Yadkori, Dávid Pál, and Csaba Szepesvári. Improved algorithms for linear stochastic bandits. In Advances in Neural Information Processing Systems 24, pages 2312–2320, 2011. Hirotugu Akaike. Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes. Annals of the Institute of Statistical Mathematics, 26(1):363–387, Dec 1974. ISSN 1572-9052. doi: 10.1007/BF02479833. 38

Scalable estimation of VARMA models

Craig F. Ansley and Robert Kohn. A note on reparameterizing a vector autoregressive moving average model to enforce stationarity. Journal of Statistical Computation and Simulation, 24(2):99–106, 1986. O. Barndorff-Nielsen and G. Schou. On the parametrization of autoregressive models by partial autocorrelations. Journal of Multivariate Analysis, 3(4):408–419, 1973. Sumanta Basu and George Michailidis. Regularized estimation in sparse high-dimensional time series models. The Annals of Statistics, 43(4):1535–1567, 2015. George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung. Time series analysis: forecasting and control. John Wiley & Sons, 2015. James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs. Software available at http://github.com/jax-ml/jax (version 0.10.0 used in this paper), 2018. Victor H. de la Peña, Tze Leung Lai, and Qi-Man Shao. Self-Normalized Processes: Limit Theory and Statistical Applications. Springer, Berlin, 2009. Francis X. Diebold and Roberto S. Mariano. Comparing predictive accuracy. Journal of Business & Economic Statistics, 13(3):253–263, 1995. Marie-Christine Düker, David S Matteson, Ruey S Tsay, and Ines Wilms. Vector autoregressive moving average models: A review. Wiley Interdisciplinary Reviews: Computational Statistics, 17(1):e70009, 2025. William Dunsmuir and Edward J. Hannan. Vector linear time series models. Advances in Applied Probability, 8(2):339–364, 1976. James Durbin and Siem Jan Koopman. Time Series Analysis by State Space Methods. Oxford Statistical Science Series. Oxford University Press, 2nd edition, 2012. Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. Bayesian Data Analysis. Chapman and Hall/CRC, Boca Raton, FL, 3rd edition, 2013. E. J. Hannan and L. Kavalieris. Multivariate linear time series models. Advances in Applied Probability, 16(3):492–561, 1984. E. J. Hannan and Jorma Rissanen. Recursive estimation of mixed autoregressive-moving average order. Biometrika, 69(1):81–94, 1982. Edward James Hannan and Manfred Deistler. The statistical theory of linear systems. SIAM, 2012. David Harvey, Stephen Leybourne, and Paul Newbold. Testing the equality of prediction mean squared errors. International Journal of Forecasting, 13(2):281–291, 1997. 39

Paulin and Elvira

Sarah E. Heaps. Enforcing stationarity through the prior in vector autoregressions. Journal of Computational and Graphical Statistics, 32(1):74–83, 2023. Lutz Kilian and Helmut Lütkepohl. Structural vector autoregressive analysis. Cambridge University Press, 2017. Béatrice Laurent and Pascal Massart. Adaptive estimation of a quadratic functional by model selection. Annals of Statistics, 28(5):1302–1338, 2000. Peter D. Lax. Functional analysis. John Wiley & Sons, 2002. Richard Lewis and Gregory C. Reinsel. Prediction of multivariate time series by autoregressive model fitting. Journal of Multivariate Analysis, 16(3):393–411, 1985. Dong C. Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical Programming, 45:503–528, 1989. Helmut Lütkepohl. New introduction to multiple time series analysis. Springer Science & Business Media, 2005. Helmut Lütkepohl. Forecasting with VARMA models. In Graham Elliott, Clive W. J. Granger, and Allan Timmermann, editors, Handbook of Economic Forecasting, volume 1, pages 287–325. Elsevier, 2006. Henry B Mann and Abraham Wald. On the statistical treatment of linear stochastic difference equations. Econometrica, 11(3/4):173–220, 1943. Konstantinos Metaxoglou and Aaron Smith. Maximum likelihood estimation of varma models using a state-space em algorithm. Journal of Time Series Analysis, 28(5):666– 685, 2007. Martin Morf, Augusto Vieira, and Thomas Kailath. Covariance characterization by partial autocorrelation matrices. The Annals of Statistics, 6(3):643–648, 1978. Gregory C Reinsel. Elements of multivariate time series analysis. Springer Science & Business Media, 2nd edition, 2003. Anindya Roy, Tucker S. McElroy, and Peter Linton. Constrained estimation of causal invertible VARMA. Statistica Sinica, 29(1):455–478, 2019. Skipper Seabold and Josef Perktold. Statsmodels: Econometric and statistical modeling with Python. In Proceedings of the 9th Python in Science Conference, pages 92–96, Austin, TX, 2010. Marc N. Spijker. On a conjecture by LeVeque and Trefethen related to the Kreiss matrix theorem. BIT Numerical Mathematics, 31(3):551–555, 1991. Lloyd N. Trefethen and Mark Embree. Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press, 2005. 40

Scalable estimation of VARMA models

Ruey S Tsay. Multivariate time series analysis: with R and financial applications. John Wiley & Sons, 2014. Ruey S. Tsay, David Wood, and Jon Lachmann. MTS R package, 2022. R package version 1.2.1 (2022); manual at https://cran.r-project.org/web/packages/MTS/MTS.pdf. Peter Whittle. The analysis of multiple stationary time series. Journal of the Royal Statistical Society: Series B, 15(1):125–139, 1953. Peter Whittle. On the fitting of multivariate autoregressions, and the approximate canonical factorization of a spectral density matrix. Biometrika, 50(1-2):129–134, 1963. Ines Wilms, Sumanta Basu, Jacob Bien, and David S Matteson. Sparse identification and estimation of large-scale vector autoregressive moving averages. Journal of the American Statistical Association, 118(541):571–582, 2023. Shuyi Zhang, Bin Guo, Anlan Dong, Jing He, Ziping Xu, and Song Xi Chen. Cautionary tales on air-quality improvement in Beijing. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2205):20170457, 2017.

Appendix A. Proof of Proposition 1 Throughout this appendix, chol(·) is the lower-triangular Cholesky factor with positive diagonal. We use two standing facts: (F1) Σ 7→ chol(Σ) and triangular inversion are realanalytic on the positive-definite cone (forward substitution solves LL⊤ = Σ entrywise, dividing only by the positive diagonal entries); (F2) a product of lower-triangular matrices with positive diagonals is again one, and is therefore the Cholesky factor of its own Gram matrix, by uniqueness of the factorization. Step 1: Stage 1 is a diffeomorphism onto the contraction ball. For LA = chol(I + −⊤ −⊤ ⊤ ⊤ = = L−1 we have P ⊤ P = L−1 M ⊤ M ) and P = M L−⊤ A (LA LA − I)LA A M M LA A ⊤ ⊤ ⊤ ⊤ −1 I − (LA LA ) . The eigenvalues of LA LA agree with those of LA LA = I + M M , namely 1 + σi (M )2 , so the eigenvalues of P ⊤ P are σi (M )2 /(1 + σi (M )2 ) < 1 and σmax (P ) < 1. Conversely, for σmax (P ) < 1 the matrix I − P ⊤ P is positive definite; chol(I −
−⊤setting−1R = ⊤ −⊤ ⊤ −1 ⊤ ⊤ −⊤ P P ) and M = P R gives I + M M = R RR + P P R = R R , whose −⊤ ⊤ Cholesky factor is R−1 by (F2), so LA = R−1 and M L−⊤ = P R R = P : the preimage A −1 exists, and it is unique because LA is determined by P through LA = chol(I − P ⊤ P ) and then M = P L⊤ A . Both directions are compositions of matrix products, triangular inversions and Cholesky factorizations of positive-definite matrices, hence smooth by (F1). Step 2: the backward pass is well defined, consistent and unique. Fix contractions f f f P 1:r and write B s := I − P s P ⊤ s ≻ 0. The recursion (6) defines Lr = I, Lr−1 , . . . , L0 , each lower triangular with positive diagonal by (F2), and each invertible, so Σfs := Lfs (Lfs )⊤ ≻ 0. The forward variance updates are consistent with this chain: by (6), Lfs−1 chol(B s ) = Lfs , so Lfs−1 B s (Lfs−1 )⊤ = Σfs , and in particular the terminal forward covariance equals I d : the anchor holds exactly. Uniqueness: if L is lower triangular with positive diagonal and L B s L⊤ = Σfs , then L chol(B s ) = chol(Σfs ) by (F2), so L = Lfs−1 ; the backward chain 41

Paulin and Elvira

is therefore the unique solution of the anchoring equations, and it depends smoothly on P 1:r by (F1). Finally, associate to the recursion the matrix sequence Γ0 := Lf0 (Lf0 )⊤ P (s) and Γs := sj=1 Aj Γs−j for s = 1, . . . , r (with Γ−k := Γ⊤ k ). With these definitions the coefficient and covariance updates of (7) are precisely the block Levinson recursion, so at (s) e (s) solve the forward and backward Yule–Walker systems every order s the stacks A1:s and A 1:s of Γ0:s with residual covariances Σfs and Σbs (Whittle, 1963; Morf et al., 1978); in particular, at s = r the output stack satisfies the order-r Yule–Walker equations of Γ0:r with residual Σfr = I d . Step 3: every contraction stack maps to a stable stack. Associate to the recursion P (s) the matrix polynomials A(s) (z) := I − sj=1 Aj z j and the reversed backward polynomials e (s) (1/z), for which the coefficient updates of (7) read B (s) (z) := z s A e (s) A(s−1) (z). B (s) (z) = z B (s−1) (z) − A s

(s−1) (z), A(s) (z) = A(s−1) (z) − z A(s) s B

Normalize by the Cholesky factors, αs (z) := (Lfs )−1 A(s) (z) and βs (z) := (Lbs )−1 B (s) (z). (s) e (s) = Lb P ⊤ (Lf )−1 and the factor updates Substituting As = Lf P s (Lb )−1 , A Lfs

s−1 s s−1 s s−1 s−1 f b b ⊤ = Ls−1 C f,s , Ls = Ls−1 C b,s with C f,s := chol(I − P s P ⊤ s ), C b,s := chol(I − P s P s )

gives
αs (z) = C −1 f,s αs−1 (z) − z P s βs−1 (z) ,  The block map Θs :=

C −1 f,s ⊤ −C −1 b,s P s

−C −1 f,s P s



C −1 b,s ⊤ −1

⊤ βs (z) = C −1 b,s z βs−1 (z) − P s αs−1 (z) . is J-unitary, Θ⊤ s J Θs = J for J = diag(I, −I):

the block identities reduce to (I − P P ) − P (I − P ⊤ P )−1 P ⊤ = I and (I − P P ⊤ )−1 P = P (I −P ⊤ P )−1 , both immediate from the Neumann series. Consequently, on |z| = 1 (where |z|2 = 1; here ∗ denotes the conjugate transpose of the matrix-valued polynomials, while ⊤ on the real coefficient blocks is the ordinary transpose), αs (z)∗ αs (z) − βs (z)∗ βs (z) = αs−1 (z)∗ αs−1 (z) − βs−1 (z)∗ βs−1 (z), and since α0 = β0 = (Lf0 )−1 is constant, αs (z)∗ αs (z) = βs (z)∗ βs (z) on the unit circle for every s. We now show by induction that det αs (z) ̸= 0 for all |z| ≤ 1. This holds at s = 0. Assume it at s−1 and set W (z) := βs−1 (z) αs−1 (z)−1 ; by the induction hypothesis det αs−1 is nonvanishing on the compact closed disk, hence on an open neighborhood of it, so W is analytic there, and W (z)∗ W (z) = I on |z| = 1 by the equal-modulus identity. For any unit vectors u, v, the scalar u∗ W (z)v is analytic on a neighborhood of the closed disk and bounded by 1 on its boundary, so the maximum principle gives ∥W (z)∥ ≤ 1 on the closed disk. Writing
αs (z) = C −1 ∥z P s W (z)∥ ≤ σmax (P s ) < 1 (|z| ≤ 1), f,s I − z P s W (z) αs−1 (z), all three factors are invertible on the closed disk, and P the induction advances. At s = r the
anchor gives Lfr = I, so A(r) (z) = αr (z) and det I − j Aj z j ̸= 0 for |z| ≤ 1: the output satisfies the root condition of Assumption 1, equivalently spr(CA ) < 1. This is the matrix form of the classical Schur–Levinson stability argument (Whittle, 1963; Morf et al., 1978). 42

Scalable estimation of VARMA models

Together with Step 2, the sequence Γ0:r generated there satisfies the order-r Yule–Walker equations of a stable stack with residual I d ; since for a stable stack this linear system (with the symmetry constraints Γ−k = Γ⊤ k ) is a rewriting of the discrete Lyapunov equation of the companion form and so has a unique solution (Lütkepohl, 2005, Section 2.1), Γ0:r coincides with the stationary autocovariances of the VAR(r) with coefficients τr (M 1:r ) and innovation covariance I d . Step 4: correspondence with a stationary autoregression. Let A1:r be stable and let z t be the stationary VAR(r) process with coefficients A1:r and innovation covariance I d . Its autocovariances Γk exist and are unique: the stacked covariance of (z t , . . . , z t−r+1 ) solves a discrete Lyapunov equation with the stable companion CA as transition P matrix k Q(C ⊤ )k is positive and state-noise covariance Q := diag(I d , 0), and the solution k≥0 CA A definite because the pair (CA , Q) is reachable: Q has rank d only, but a shock entering the first block propagates through every lag block, so the first r propagation blocks span all rd state coordinates. For 0 ≤ s ≤ r let Σfs and Σbs denote the covariances of the errors of the best linear prediction of z t from its s most recent past values and of z t−s from the same block of values, respectively, and let ∆s be the associated partial cross-covariance. The multivariate Levinson–Durbin recursion computes these quantities from Γ0:r , and its coefficient updates coincide with the middle lines of (7) upon setting P s := (Lfs−1 )−1 ∆s (Lbs−1 )−⊤ ,

Lfs−1 = chol(Σfs−1 ),

Lbs−1 = chol(Σbs−1 ),

(s)

b b −1 −1 since then As = Lfs−1 P s (Ls−1 )  = ∆s (Σs−1 ) , the usual reflection-coefficient update.

The joint covariance

Σfs−1 ∆s

b ∆⊤ s Σs−1

of the two prediction errors is positive definite (it is

the Schur complement of the covariance of the conditioning block (z t−1 , . . . , z t−s+1 ) in the positive-definite covariance of the full vector (z t , . . . , z t−s )), and after congruence by
diag (Lfs−1 )−1 , (Lbs−1 )−1 its positive-definiteness is, by the Schur complement, equivalent f ⊤ to σmax (P s ) < 1. The recursion’s variance downdates are Σfs = Lfs−1 (I − P s P ⊤ s )(Ls−1 ) b ⊤ and Σbs = Lbs−1 (I − P ⊤ s P s )(Ls−1 ) , matching Step 2 and the last line of (7). For a VAR(r) the best linear predictor from r lags is the autoregression itself, so Σfr = I d and the order-r predictor coefficients are A1:r ; by stationarity, Σb0 = Σf0 = Γ0 . The contraction stack extracted from A1:r therefore satisfies the anchoring equations of Step 2, whose solution is unique, and running the forward pass (7) on it returns the order-r predictor coefficients, that is, A1:r itself. This gives surjectivity. Injectivity and smoothness now follow together. The stable coefficients and innovation covariance I d determine the Lyapunov solution Γ0:r uniquely, and by the last part of Step 3 this is the sequence generated in Step 2 from the input contractions; the Levinson extraction from Γ0:r is deterministic and, by the anchoring uniqueness of Step 2, returns those contractions, so with Step 1 the preimage is unique, and the inverse is the composition A 7→ Γ(A) 7→ P (Γ) 7→ M (P ). Each constituent map is smooth on the open stable set: the first is the solution of a linear system with invertible, polynomially A-dependent coefficient matrix; the second stays on positive-definite covariances and applies (F1) at every order; the third is the closed form of Step 1. The forward direction is smooth by Steps 1–2, so τr is a smooth diffeomorphism onto the stable set. □ We record one consequence used in Assumption 2(ii) and in Appendix C.6: padding a stable stack A1:p̄ with zero matrices up to order p ≥ p̄ leaves the lag polynomial unchanged, 43

Paulin and Elvira

and the order-p companion matrix has characteristic polynomial z (p−p̄)d det(z p̄ I −. . .), so its spectrum is that of the original companion matrix together with (p − p̄)d zero eigenvalues. Stability is preserved, and Proposition 1 applied at order p gives the padded stack unique order-p PAR coordinates.

Appendix B. Proof of Proposition 3 Throughout, ∥ · ∥ denotes the operator (ℓ2 → ℓ2 ) norm (which coincides with the ∥ · ∥op of Proposition 3), σmin (·) the smallest singular value, and M ∗ the conjugate transpose of M . Introduce the infinite-VAR operator of the VARMA(p, q) and its degree-m truncation, A(z) := I −

X

Πj z j = Θ(z)−1 Φ(z),

Am (z) := I −

m X

Πj z j ,

j=1

j≥1

where the second equality is the operator form of the infinite-VAR representation (9), obtained from the VARMA equation (1) (so Π0 = −I and at = A(B)Z t ), and Am is exactly the autoregressive operator of the truncated VARMA(p, q, m) of Definition 2. We use the constants Gρ , µmin , β• , β◦ , m0 of Proposition 3. Under Assumption 1, det A(z) = det Φ(z)/ det Θ(z) ̸= 0 for all |z| ≤ 1 and A is continuous on the compact disk, so β• , β◦ √ are both strictly positive, with β• ≤ β◦ and β◦ = µmin . Step 1 (geometric decay of the dropped coefficients, via a Cauchy estimate). The zeros of det Θ(z) are the reciprocals of the nonzero eigenvalues P of CΘ , jhence have modulus ≥ 1/spr(CΘ ); since Φ is a matrix polynomial, I − A(z) = j≥1 Πj z is analytic on the disk |z| < 1/spr(CΘ ). Fix ρ ∈ (spr(CΘ ), 1), so the circle |z| = 1/ρ lies strictly inside that disk. Cauchy’s estimate for Taylor coefficients then gives, for every j ≥ 1, I
1 Πj = I − A(z) z −j−1 dz, hence ∥Πj ∥ ≤ ρ j max ∥I − A(z)∥ = Gρ ρ j , 2πi |z|=1/ρ |z|=1/ρ with Gρ as in (19), finite because Θ(z) is invertible on |z| = 1/ρ < 1/spr(CΘ ). This is the explicit, non-asymptotic counterpart of Gelfand’s formula (18), whose limit only asserts that such a constant exists (Trefethen and Embree, 2005). Consequently, for every m ≥ 0, X X Gρ ρ m δm := ∥Πj ∥ ≤ Gρ ρj = ρ . (43) 1−ρ j>m

j>m

Step 2 (stationarity of the truncation and a uniform inverse bound). For |z| ≤ 1, ∥A(z) − P j ≤ δ m+1 /(1 − ρ), the threshold m Am (z)∥ = m by (43). Since δm = Gρ ρ 0 j>m Πj z of (20) is exactly the smallest m ≥ p + q with δm ≤ β• /2; hence δm ≤ β• /2 ≤ β◦ /2 for all m ≥ m0 . For such m, Weyl’s inequality gives, on the closed disk, σmin (Am (z)) ≥ β• − δm ≥ β• /2 > 0, so det Am (z) ̸= 0 on |z| ≤ 1 and the truncated VARMA(p, q, m) is stationary, with a well-defined spectral density; and on the unit circle, σmin (Am (e−iω )) ≥ β◦ − δm ≥ β◦ /2, whence ∥Am (e−iω )−1 ∥ ≤ 2/β◦ . Step 3 (spectral-density gap, controlled by the on-circle conditioning). For m ≥ m0 both processes are stationary, driven by the same white noise of covariance Σ, with matrix spectral densities 1 A(e−iω )−1 Σ A(e−iω )−∗ , f (ω) = 2π

1 f m (ω) = 2π Am (e−iω )−1 Σ Am (e−iω )−∗ .

44

Scalable estimation of VARMA models

Abbreviate A := A(e−iω ), Am := Am (e−iω ); on the unit circle ∥A−1 ∥ = ∥A−∗ ∥ ≤ β◦−1 and −1 −1 −1 ∥A−1 − A−1 m ∥ ≤ 2/β◦ by Step 2. The resolvent identity A m = A (Am − A)Am gives −2 ∥A−1 − A−1 m ∥ ≤ 2β◦ δm , and writing −∗ −1 −∗ −∗ A−1 ΣA−∗ − A−1 − A−1 + A−1 − A−∗ m ΣAm = (A m )ΣA m Σ(A m )

yields, uniformly in ω ∈ [−π, π], ∥f (ω) − f m (ω)∥ ≤

 3∥Σ∥ 1  −2 2β◦ δm ∥Σ∥ β◦−1 + 2β◦−1 ∥Σ∥ 2β◦−2 δm = δm . 2π πβ◦3

−1/2

(The role of β◦−1 = µmin as the on-circle conditioning is discussed in the remark below.) Step 4 (autocovariances). The autocovariances are the Fourier coefficients of the spectral R π ilω Rπ (m) = −π eilω f m (ω) dω. Hence, for every l ≥ 0 and densities, Γ∞ l = −π e f (ω) dω and Γl 3/2

every m ≥ m0 , using Step 3, β◦3 = µmin , and (43), Z π 6∥Σ∥ Gρ ρ m 6∥Σ∥ (m) ∞ ∥f (ω) − f m (ω)∥ dω ≤ 3/2 δm ≤ ≤ Γl − Γl ρ = Cρ ρm . 3/2 −π µmin (1 − ρ) µmin This is exactly (21) with the explicit constant stated there, uniform in l, completing the proof. □ Remark 10 (Dimension dependence) The two constants are explicit but modeldependent: Gρ reflects the non-normality of the moving-average inversion Θ−1 Φ, and µmin is the on-circle floor of the squared smallest singular value of the autoregressive operator A(e−iω ); together with ∥Σ∥ it yields the spectral-density upper bound of Basu and Michailidis (2015) (it is not itself an eigenvalue floor of the spectral density); neither carries an explicit dimension factor. A fully a priori, dimension-explicit decay constant is k ∥ ≤ e (qd) K(C /ρ) ρk with K(·) the instead available from the Kreiss matrix theorem, ∥CΘ Θ Kreiss constant (Spijker, 1991), where the linear factor qd (the order of CΘ ) is unavoidable in the worst case; the resolvent constant Gρ above is the sharper, model-adapted alternative.

Appendix C. Proof of Theorem 6 Throughout, ∥ · ∥ is the operator norm, ∥ · ∥F the Frobenius norm, N = T − m, and C, c′ > 0 denote constants depending only on the quantities listed in Theorem 6, whose values may change from occurrence to occurrence; the constants c, m1 , T0 of the theorem are fixed by finitely many such requirements collected along the proof. We abbreviate Λ := Λδ and use the standing conditions T ≥ T0 , m1 ≤ m ≤ T /2 (so N ≥ T /2) and Λ ≤ cT without further comment. Because B, p, q, d, and the true-model quantities of Assumption 2 are fixed, every quantity introduced in Lemmas 11 and 12 (ρ∗ , Ḡ, L̄, Π̄, q̄, and their companions) is a constant independent of T and m throughout this appendix. For ϑ ∈ KB let (Φ1:p (ϑ), Θ1:q (ϑ)) := T (ϑ), let Πj (ϑ) be the induced infinite-VAR coefficients (10)–(11), and, as in Appendix B, X A(z; ϑ) := I − Πj (ϑ)z j = Θ(z; ϑ)−1 Φ(z; ϑ). j≥1

45

Paulin and Elvira

Write Π0j := Πj (ϑ0 ), A0 := A(· ; ϑ0 ) and ∆j (ϑ) := Πj (ϑ) − Π0j . Zero-padding does not change the lag polynomials, so Π01:∞ are exactly the infinite-VAR coefficients of the datagenerating process, at = A0 (B)z t . Let Ft := σ(as : s ≤ t); under Assumption 1 the stationary solution is causal, so z t is Ft -measurable and at+1 is independent of Ft . C.1 Preliminaries Lemma 11 (Uniform geometry on KB ) There exist ρ∗ ∈ (0, 1) and Ḡ, L̄ < ∞, depending ϑ′ ∈ KB : (i) ∥Πj (ϑ)∥ ≤ Ḡρ∗j for every j ≥ 1, so P only on (B, p, q, d), such that for all ϑ,P ′ ′ j≥1 ∥Πj (ϑ)∥ ≤ Π̄ := Ḡρ∗ /(1 − ρ∗ ); (ii) j≥1 ∥Πj (ϑ) − Πj (ϑ )∥F ≤ L̄ ∥ϑ − ϑ ∥2 . Proof T is a smooth diffeomorphism onto the open stable set (Proposition 1) and spr(·) is continuous, so ρ̄ := maxϑ∈KB max{spr(CΦ(ϑ) ), spr(CΘ(ϑ) )} is attained on the compact ball and is < 1; set ρ∗ := (1 + ρ̄)/2. On the compact set {|z| ≤ 1/ρ∗ } × KB the polynomial det Θ(z; ϑ) is continuous and zero-free (its zeros have modulus at least 1/ρ̄ > 1/ρ∗ ), so ∥Θ(z; ϑ)−1 ∥ is uniformly bounded there, (z, ϑ) 7→ A(z; ϑ) = Θ(z; ϑ)−1 Φ(z; ϑ) is C 1 , and Ḡ :=

max

|z|=1/ρ∗ , ϑ∈KB

I − A(z; ϑ) ,

L0 :=

max

|z|=1/ρ∗ , ϑ∈KB

∂ϑ A(z; ϑ)

are finite, ∂ϑ A denoting the differential of ϑ 7→ A(z; ϑ) as a linear map from (Rnpar , ∥ · ∥2 ) to (Cd×d , ∥ · ∥F ). The Cauchy coefficient estimate of Step 1 of Appendix B, on the circle |z| = 1/ρ∗ , gives (i); applied to A(· ; ϑ) − A(· ; ϑ′ ) together with the mean value inequality along the segment [ϑ, ϑ′ ] in the convex ball KB , it gives ∥Πj (ϑ)−Πj (ϑ′ )∥F ≤ ρ∗j L0 ∥ϑ−ϑ′ ∥2 , and summing the geometric series yields (ii) with L̄ := L0 ρ∗ /(1 − ρ∗ ).

Lemma 12 (Second moments of filtered processes) Set c0 :=

λmin (Σ) , (1 + Π̄)2

C0 :=

∥Σ∥ , µmin

with of Proposition 3 forP the true model. Let {B j }P j≥1 satisfy P P µmin the on-circle constant 2 2 2 j ∥B j ∥ < ∞ and st := j≥1 B j z t−j . Then (a) c0 j ∥B j ∥F ≤ E∥st ∥ ≤ C0 j ∥B j ∥F ; and (b) for any n, t1 , the covariance matrix of the stacked Gaussian vector (st1 , . . . , st1 +n−1 )
2 P has operator norm at most C0 j ∥B j ∥ . 1 Proof The spectral density of z is f (ω) = 2π A0 (e−iω )−1 Σ A0 (e−iω )
−∗ (Appendix B, Step 3), with on-circle eigenvalue bounds λmin (f ) ≥ λmin (Σ)/ 2π∥A0 ∥2 ≥ c0 /(2π), since max|z|=1 ∥A0 (z)∥ ≤ 1 + Π̄ by Lemma 11(i), and λmax (f ) ≤ ∥Σ∥/(2πµmin ) = C0 /(2π), the bound of Basu and Michailidis (2015). The filtered process is stationary Gaussian P −ijω . Part (a) is Parseval’s efB e ∗ , where B(ω) e with spectral density f s = B := j Bj e R P π 1 2 e 2 = identity 2π j ∥B j ∥F combined with the two-sided eigenvalue bounds −π ∥B∥F dω Rπ 2 on f inside E∥st ∥ = −π tr f s dω. Part (b): the stacked covariance is a principal submatrix of the bi-infinitePblock-Toeplitz covariance of s, whose operator norm is at most
2 2π ess supω λmax (f s ) ≤ ∥B ∥ C (Basu and Michailidis, 2015, Proposition 2.3). 0 j j

46

Scalable estimation of VARMA models

We use P two2 standard probabilistic facts. First (Laurent and Massart, 2000, PLemma 1): if PQ 2= i λi ξi with fixed λi ≥ 0 and ξi i.i.d. N (0, 1), then, writing Σ1 := i λi , Σ2 := i λi ≤ λmax Σ1 , p p

P Q ≤ Σ1 − 2 Σ2 x ≤ e−x , x > 0. (44) P Q ≥ Σ1 + 2 Σ2 x + 2λmax x ≤ e−x , Second, a scalar self-normalized martingale inequality in the spirit of Abbasi-Yadkori et al. (2011, Theorem 1) (see also de la Peña et al., 2009): if ξt is Ft+1 -measurable with E[eλξt | 2 2 Ft ] ≤ eλ st /2 for all λ ∈ R and Ft -measurable st , then for every δ ′ ∈ (0, 1), with probability at least 1 − δ ′ , q X X √
ξt ≤ 2 (1 + V ) log 1 + V /δ ′ , s2t , (45) V := t

t

the sums running over any fixed index window. This scalar specializationP of the method of mixtures follows by integrating the exponential supermartingale exp(λ t ξt − λ2 V /2) against a standard Gaussian in λ (the unit mixing variance is the source of the two 1 + V factors) and applying Markov’s inequality to the resulting supermartingale; see AbbasiYadkori et al. (2011) and de la Peña et al. (2009). C.2 Decomposition and high-probability events For t = m, . . . , T − 1 define v t (ϑ) :=

m X

Π0j − Πj (ϑ) z t+1−j ,

wt :=

j=1

X

Π0j z t+1−j ,

j>m

P so that the residual rt (ϑ) = z t+1 − j≤m Πj (ϑ)z t+1−j appearing in Lm and R decomposes exactly as rt (ϑ) = at+1 + v t (ϑ) + wt , with v t (ϑ0 ) = 0. Set X X X M (ϑ) := v t (ϑ)⊤ at+1 , C a (ϑ) := v t (ϑ)a⊤ W (ϑ) := v t (ϑ)⊤ wt , t+1 , t

Qv (ϑ) :=

t

X

∥v t (ϑ)∥2 ,

Qw :=

t

X

t

∥wt ∥2 ,

Qa :=

t 2

qm (ϑ) := E∥v t (ϑ)∥ ,

X

∥at+1 ∥2 ,

t 2

qw := E∥wt ∥ .

Expanding the squares and outer products yields the exact identities Lm (ϑ) − Lm (ϑ0 ) = 2M (ϑ) + 2W (ϑ) + Qv (ϑ), R(ϑ) − R(ϑ0 ) = V (ϑ) + C(ϑ) + C(ϑ)⊤ ,

(46) V :=

X t

vtv⊤ t ,

C := C a +

X

v t w⊤ t .

t

(47) √ Note M√ = tr C a (so |M P | ≤ ⊤ d ∥C a ∥√F ), tr V = Qv , and, pathwise by Cauchy–Schwarz, |W | ≤ Qw Qv and ∥ t v t wt ∥F ≤ Qv Qw . By Lemmas 11 and 12, applied with B j = 47

Paulin and Elvira

(Π0j − Πj (ϑ))1{j ≤ m} and B j = Π0j 1{j > m} respectively, we have for all ϑ ∈ KB c0

X

X

∥∆j (ϑ)∥2F ≤ qm (ϑ) ≤ C0

j≤m

j≤m



λmax Cov v m (ϑ), . . . , v T −1 (ϑ)



≤ C0

4dC0 Ḡ2 , 1 − ρ2∗

(48)

≤ 4C0 Π̄2 =: λ̄v ,

(49)

∥∆j (ϑ)∥2F ≤ q̄ :=

X

∥∆j ∥

2

j≤m

qw ≤ Cρ2m ∗ ,


 λmax Cov wm , . . . , wT −1 ≤ Cρ2m ∗ ,

(50)

√ √ P using ∥∆j ∥ ≤ ∥∆j ∥F ≤ d ∥∆j ∥ ≤ 2 d Ḡρ∗j and j≥1 ∥∆j ∥ ≤ 2Π̄ in (49). npar , so that log |N | ≤ Fix a τT -net N ⊂ KB with τT := T −4 and |N | ≤
(3B/τT ) 2 5 npar log T for T ≥ T0 , and set x∗ := log 10 d |N | T /δ ≤ CΛ. Let H := H1 ∩ · · · ∩ H5 , where each event below has probability at least 1 − δ/5 (so P(H) ≥ 1 − δ): H1 : max1≤t≤T ∥z t ∥2 ∨ maxm<t≤T ∥at ∥2 ≤ Z 2 := C (1 + log(T /δ)). [Each ∥z t ∥2 (resp. ∥at ∥2 ) is a Gaussian quadratic form with Σ1 ≤ dC0 (resp. tr Σ); apply the upper tail of (44) with x = log(10T /δ) and a union bound over the at most 2T vectors.] H2 : for all ϑ ∈ N : 12 N qm (ϑ) − CΛ ≤ Qv (ϑ) ≤ 2 N qm (ϑ) + CΛ. [ (44) for the Gaussian quadratic form Qv (ϑ), with √ Σ1 =1 N qm (ϑ) and λmax ≤ λ̄v by (49), at x = x∗ ≤ CΛ; absorb the deviation via 2 ab ≤ 2 a + 2b; union over N .] H3 : Qw ≤ Cρ2m ∗ N and Qa ≤ 2N tr Σ. [Upper tail of (44) with (50), x = x∗ ≤ CΛ ≤ 2CcN .] H4 : for all ϑ ∈ N and 1 ≤ k, l ≤ d: q p

(C a (ϑ))kl ≤ 2 1 + ∥Σ∥Qv (ϑ) log 1 + ∥Σ∥Qv (ϑ) + x∗ .
P [(C a )kl = t (v t )k (at+1 )l ; conditionally on Ft the summand N 0, (v t )2k Σll P is exactly with (v t )k Ft -measurable, so (45) applies with V = Σll t (v t )2k ≤ ∥Σ∥Qv (ϑ); union over the d2 |N | pairs.] √ P ⊤ H5 : t at+1 at+1 − N Σ F ≤ C N Λ. [Entrywise by (44) at x = x∗ ≤ CΛ ≤ CN . The diagonal entries are Gaussian quadratic forms; for k ̸= l, polarization writes the summand as a difference of two such forms, (at+1 )k (at+1 )l = 14 ((at+1 )k + (at+1 )l )2 −  2 ((at+1 )k − (at+1 )l ) , and (44) applies to each. Every entry then deviates by at most √ √ C( N x∗ + x∗ ) ≤ C N Λ, and a union bound over the d2 entries with ∥ · ∥F ≤ d maxk,l |·kl | gives the event. H5 is not unioned over N , so a tail parameter C log(d/δ) would suffice; x∗ is kept for uniform bookkeeping.] ′ From the net to KB . Let ϑ ∈ KB and ϑ′ ∈ N with ∥ϑ − τT . On H1 ∩ H3 , Pϑ ∥2 ≤ ⊤ every functional appearing in H2 –H4 (namely M , C a , Qv , W , t v t wt and N qm ) changes by at most 1 between ϑ and ϑ′ , for T ≥ T0 . One representative case: by Lemma 11(ii), ′ ∥v ≤ 2Π̄Z on H1 , so, explicitly, |Qv (ϑ) − Qv (ϑ′ )| ≤ P t (ϑ) − v t (ϑ )∥ ≤ ′L̄ τ
T Z and ∥v t (ϑ)∥ ′ 2 −3 ≤ T −1 ≤ 1 for T ≥ T , 0 t ∥v t (ϑ)∥ + ∥v t (ϑ )∥ ∥v t (ϑ) − v t (ϑ )∥ ≤ CT Z Π̄L̄ τT ≤ CΛ T 2 using Z ≤ CΛ from H1 . The remaining differences are bounded the same way (W and

48

Scalable estimation of VARMA models

√ Qw ≤ C T on H3 ). For H4 both sides of the inequality move: the left side changes by at most 1 as above, while on the right Qv changes q p by at most 1 and the map x 7→ 2(1 + ∥Σ∥x)(log 1 + ∥Σ∥x + x∗ ) is increasing in x ≥ 0,

so, writing F for that map, |(C a (ϑ))kl | ≤ |(C a (ϑ′ ))kl |+1 ≤ F Qv (ϑ′ ) +1 ≤ F Qv (ϑ)+1 +1 on H4 , and the transition from H4 to (U2) absorbs both the shift by 1 and the additive √ constant. Finally, N |qm (ϑ) − qm (ϑ′ )| ≤ N · 2 q̄C0 L̄ τT ≤ CT −3 ≤ 1 by Lemma 12(a). Consequently, on H, for all ϑ ∈ KB : ⊤ t v t w t via Cauchy–Schwarz against

P

(U1)

1 2 N qm (ϑ) − CΛ

(U2)

∥C a (ϑ)∥F ≤ C

(U3)

Qw ≤ Cρ2m ∗ N;

q

√

≤ Qv (ϑ) ≤ 2N qm (ϑ) + CΛ ≤ 2N q̄ + CΛ ≤ CN ; √
1 + Qv (ϑ) Λ + 1, |M (ϑ)| ≤ d ∥C a (ϑ)∥F ; p (U4) ∥C(ϑ)∥F ≤ ∥C a (ϑ)∥F + Qv (ϑ) Qw .

For (U2) we used that, on H2 , ∥Σ∥Qv ≤ CT at net points, so the logarithm in H4 is at most C log T + x∗ ≤ CΛ, and ∥C a ∥F ≤ d maxkl |(C a )kl |. Finally, the penalty in (27) 2 ) on K , where σ 2 and (34) satisfies 0 ≤ pen(ϑ) ≤ PB := B 2 /(2σmin B min is the variance floor of Assumption 2(iii), since the group sums of squares add up to ∥ϑ∥22 ≤ B 2 . This bound is the only place the group variances enter the proof, and it holds pathwise for any datadependent selection respecting the floor: the events H1 , . . . , H5 involve the data alone, so the conclusion covers cross-validated variances without any modification of the argument. C.3 Proof for the RLS estimator RLS RLS 0 ) + η (as ϑ0 ∈ K ), and LRLS = Since ϑ̂ ∈ KB has LRLS B m m (ϑ̂) ≤ minKB Lm + η ≤ Lm (ϑ
1 0 ) ≤ 2 pen(ϑ0 ) − pen(ϑ̂) + 2η ≤ 2P + 2η. Abbreviating L + pen, we get L ( ϑ̂) − L (ϑ m m m B 2 √ Q̂ := Qv (ϑ̂) and using (46), |W | ≤ Qw Qv , (U2) and (U3), on H:  q √  q Q̂ ≤ 2|M (ϑ̂)| + 2|W (ϑ̂)| + 2PB + 2η ≤ 2 d C (1 + Q̂)Λ + 1 + 2 Cρ2m ∗ N Q̂ + 2PB + 2η. √ The two square roots absorb into the left side after rescaling by
their constants: uv ≤
1 2 2 2m 2 (u + v) with (u, v) = Q̂/8, 8C Λ and (u, v) = Q̂/8, 8C ρ∗ N gives q √ q ′ 2m 1 2 d C (1 + Q̂)Λ ≤ 81 Q̂ + C ′ Λ + C ′ , 2 Cρ2m ∗ N Q̂ ≤ 8 Q̂ + C ρ∗ N, q q the first after 1 + Q̂ ≤ 1 + Q̂, so the Q̂-contributions total at most 14 Q̂. Rearranging,
(51) Q̂ ≤ C Λ + ρ2m ∗ N +η ,

after absorbing PB + 1 ≤ CΛ. By (U1), qm (ϑ̂) ≤ 2(Q̂ + CΛ)/N ≤ C(Λ/N + ρ2m ∗ + η/N ), so by (48) and N ≥ T /2, Λ X η 2m 2 ∥∆j (ϑ̂)∥F ≤ qm (ϑ̂)/c0 ≤ C + ρ∗ + , T T j≤m X X
2 ∥∆j (ϑ̂)∥2F ≤ d 2Ḡρ∗j ≤ Cρ2m ∗ , j>m

the tail using ∥M ∥F ≤ displays gives (36).

√

j>m

d ∥M ∥ and Lemma 11(i) for both ϑ̂ and ϑ0 . Summing the two □ 49

Paulin and Elvira

C.4 Proof for the JMAP estimator Write S 0 := Ψ0 + R(ϑ0 ) and Ŝ := Ψ0 + R(ϑ̂); both are positive definite since Ψ0 ≻ 0. The marg 0 estimate has Lmarg + η ≤ Lmarg m (ϑ̂) ≤ minKB Lm m (ϑ ) + η; with 0 ≤ pen ≤ PB , dividing by (ν0 + N )/2 ≥ N/2 gives log det Ŝ − log det S 0 ≤

2PB + 2η . N

(52)

The argument parallels Appendix C.3: the basic inequality (52) caps how much the fitted log-determinant may exceed the reference one; localizing around S 0 splits the difference into a positive-semidefinite signal term, whose trace is proportional to the prediction-error mass Q̂ := Qv (ϑ̂), and a cross term that the events of Appendix C.2 render negligible in operator norm; the resulting scalar inequality as P before. P in Q̂ is then closed exactly 0 0 ⊤ ⊤ Spectral bounds for S . Since t ) and t w t w t ⪰ 0, t (at+1 + w t )(at+1 + w√ PR(ϑ ) = ⊤ we have, on H3 ∩ H5 , using ∥ t P at+1 at+1 − N Σ∥ √ ≤ ∥ · ∥F ≤ C m N Λ from H5 and the ⊤∥ ≤ a w bilinear Cauchy–Schwarz bound ∥ q Qa Qw ≤ CN ρ∗ (for unit vectors u, v: t t+1 t p
P P P √ ⊤ ⊤ 2 ⊤ 2 u⊤ Qa Qw , with Qa , Qw bounded on H3 ), t at+1 w t v ≤ t (u at+1 ) t (w t v) ≤ and noting that the cross term appears together with its transpose, hence with a factor 2, p √ √ λmin (S 0 ) ≥ N λmin (Σ) − C N Λ − 2 Qa Qw ≥ N λmin (Σ) − C N Λ − CN ρm ∗ ≥ N cΣ , √ λmax (S 0 ) ≤ ∥Ψ0 ∥ + N ∥Σ∥ + C N Λ + CN ρm ∗ + Qw ≤ N CΣ , with p cΣ := λmin √ (Σ)/2 and CΣ := 3∥Σ∥, for T ≥ T0 , m ≥ m1 and c small enough that 1 C Λ/N ≤ C 2c ≤ λmin (Σ)/4 and Cρm ∗ ≤ λmin (Σ)/4. Log-determinant localization. By (47), Ŝ = S 0 + V (ϑ̂) + C(ϑ̂) + C(ϑ̂)⊤ , so factoring out S 0 1/2 on both sides, log det Ŝ − log det S 0 = log det(I + E V + E C ), E V := S 0 −1/2 V (ϑ̂)S 0 −1/2 ,

E C := S 0 −1/2 C(ϑ̂) + C(ϑ̂)⊤ S 0 −1/2 .

The two normalized perturbations play different roles. The signal term E V ⪰ 0 can only increase the log-determinant, and its trace measures the prediction-error mass: since tr(S 0 −1 V ) ∈ [tr V /λmax (S 0 ), tr V /λmin (S 0 )] for V ⪰ 0 and tr V = Q̂, Q̂ Q̂ ≤ tr E V ≤ ≤ C, N CΣ N cΣ the last bound because Q̂ ≤ CN by (U1). The cross term E C collects the martingale interaction of v with the innovations and its interaction with theqtruncation tail;qon the event it is small in operator norm: by (U2) and (U4), ∥C(ϑ̂)∥F ≤ C

(1 + Q̂)Λ+1+

which with the crude bounds Q̂ ≤ CN and Qw ≤ Cρ2m ∗ N (U3) gives ∥C(ϑ̂)∥F ≤ C( m N ρ∗ ), whence e := ∥E C ∥ ≤

p
∥C(ϑ̂) + C(ϑ̂)⊤ ∥ 2∥C(ϑ̂)∥F ≤ ≤ C Λ/N + ρm ≤ 12 ∗ 0 N cΣ λmin (S ) 50

Q̂ Qw , √ NΛ+

Scalable estimation of VARMA models

p √ for c small, m ≥ m1 , T ≥ T0 (as before, Λ/N ≤ 2c). The mechanism is now visible: (52) allows the log-determinant to grow by at most 2PB /N , while the signal term forces growth of order tr E V ≳ Q̂/N that the small cross term cannot cancel. We make this quantitative in three elementary steps. Step 1 (absorb the cross term): E C ⪰ −eI, so by Loewner monotonicity of log det over positive-definite matrices, log det(I + E V + E C ) ≥ log det (1 − e)I + E V

=

d X

log(1 − e + λi ),

i=1

where λ1 , . . . , λd ≥ 0 are the eigenvalues of E V (the shifted matrix is positive definite since e ≤ 21 ). Step 2 (separate the shift): 1 − e + λi ≥ (1 − e)(1 + λi ) (expanding, this is just eλi ≥ 0) and log(1 − e) ≥ −2e on [0, 12 ], so d X

d X

log(1 − e + λi ) ≥

i=1

log(1 + λi ) − 2de.

i=1

Step 3 (curvature): by concavity, log(1 + x) ≥ x log(1 + b)/b ≥ x/(1 + b) for 0 ≤ x ≤ b; taking b = maxi λi ≤ tr E V , d X

P

log(1 + λi ) ≥

i=1

tr E V Q̂ tr E V i λi ≥ ≥ , ≥ 1 + maxi λi 1 + tr E V C CN

the penultimate inequality because tr E V ≤ C as shown above. Combining the three steps with (52) and re-inserting the sharp form of ∥C(ϑ̂)∥F into 2de ≤ 4d∥C(ϑ̂)∥F /(N cΣ ), q q  Q̂ 2PB + 2η C ≤ + 2de ≤ PB + η + (1 + Q̂)Λ + 1 + Q̂ ρ2m N . ∗ CN N N √ √ √ Multiplying through by CN and applying ab ≤ a4 + b (valid since ( a/2 − b)2 ≥ 0) to each square root (with a = 1 + Q̂, b = C 2 Λ for the first and a = Q̂, b = C 2 ρ2m ∗ N for the second, so the a-contributions total 21 Q̂ + 14 and are absorbed into the left-hand side) yields

2m Q̂ ≤ C Λ + ρ2m ∗ N + η + PB + 1 ≤ C Λ + ρ∗ N + η , which is (51). The conclusion then follows exactly as in Appendix C.3.

□

C.5 Proof of Corollary 7, and optimality For any J ≥ 1, by Cauchy–Schwarz on the first J lags and Lemma 11(i) on the tail (both ϑ̂, ϑ0 ∈ KB ), 1/2 2Ḡρ J+1 X √ X ∗ ∥∆j (ϑ̂)∥2F + . (53) ∥∆j (ϑ̂)∥ ≤ J 1 − ρ∗ j≥1

j≥1

(ii) If m ≥ log T / log(1/ρ∗ ) then ρ2m ≤ T −2 ≤ Λδ /T , and on the event ET of the state∗ ment also η/T ≤ C ′ npar log T /T ≤ C ′ Λδ /T . On the intersection the theorem’s Pof ET with 2 ≤ CΛ /T ; takevent, which has probability at least 1 − δ − δ ′ , (36) gives ∥∆ ∥ j δ j √ √ F ing J := ⌈log T /(2 log(1/ρ∗ ))⌉, so that ρJ∗ ≤ T −1/2 and J ≤ C log T , (53) is at most 51

Paulin and Elvira

p p
Λδ log T /T + T −1/2 ≤ C Λδ log T /T . (i) Fix δ ∈ (0, 1/2), ϵ > 0 and t > 0. Under C m → ∞ and m ≤ T /2 the standing conditions of the theorem hold for all large T (for fixed δ, Λδ = o(T )), and βT := C(Λδ /T + ρ2m ∗ + t) → C t. On the intersection of the theorem’s −1/2 event with {η/T ≤ t}, the bound (36) is at most βT , and (53) with J = JT := ⌈βT ⌉ is at 1/4 most CβT + CρJ∗ T , which is smaller than ϵ for all large T once t is chosen small enough (depending on ϵ), since ρ∗ < 1 is fixed. Hence X  lim sup P ∥∆j ∥ > ϵ ≤ δ + lim sup P(η/T > t) = δ, T

T

j

P

using η/T → 0; since δ was arbitrary, the sum tends to 0 in probability.

□

Remark 13 (The T −1/2 dependence is unavoidable) The T −1/2 dependence in Corollary 7(ii) cannot be improved, already in a one-parameter Gaussian submodel; we make no claim of minimax optimality in npar for the summed operator-norm loss, which would require a matching dimension-dependent packing bound. Take Θ = 0, Σ √ = I d , Φ1 = γ e1 e⊤ 1 (so Π1 = Φ1 ), and the two hypotheses γ ∈ {0, h} with h := 1/ T , each observed from its stationary initial law. The Gaussian chain rule bounds the Kullback–Leibler P divergence between the two path laws by 12 Tt=2 E0 [(hzt−1,1 )2 ] plus the initial-state term 1 2 2 )) ≤ h4 , hence by T h2 /2 + h4 ≤ 1. Pinsker’s inequality gives total 2 (−h − log(1 − hp Le Cam’s two-point method yields a testing error of at least variation p at most 1/2, so P (1 − 1/2)/2 ≥ 0.14. Since j≥1 ∥Π̂j − Πj ∥ ≥ ∥Π̂1 − Π1 ∥ and the two candidate values √ P of Π1 are h apart, every estimator incurs j≥1 ∥Π̂j − Πj ∥ ≥ 1/(2 T ) with probability at least 0.14 under one of the two models. C.6 Extension to the multiplicative seasonal VARMA Proof [Proof of Corollary 9] Only Lemma 11 and the well-specification argument interact with the parametrization; we verify these, after which every step of Appendices C.1–C.5 applies verbatim with npar = (p + q + P + Q)d2 , which enters the proof only through the net cardinality and hence Λδ . (i) Expanded stability and uniform geometry. For ϑ ∈ KB write Φ̄(z; ϑ) := Φ(s) (z s ; ϑ) Φ(z; ϑ),

Θ̄(z; ϑ) := Θ(s) (z s ; ϑ) Θ(z; ϑ)

for the expanded operators of (35), and A(z; ϑ) := Θ̄(z; ϑ)−1 Φ̄(z; ϑ). The four per-factor PAR maps are smooth, so, exactly as in the proof of Lemma 11,  ρ̄ := max max spr(CΦ ), spr(CΦ(s) ), spr(CΘ ), spr(CΘ(s) ) ϑ∈KB

is attained and < 1. Factorwise stability implies expanded stability because determinants multiply, det Θ̄(z; ϑ) = det Θ(s) (z s ; ϑ) det Θ(z; ϑ), and the seasonal factor’s zeros satisfy |z|s ≥ 1/ρ̄, while the regular factor’s zeros have modulus at least 1/ρ̄ ≥ (1/ρ̄)1/s (as s ≥ 1);
1/s hence every zero of det Θ̄ has modulus at least (1/ρ̄)1/s > 1, and ρ∗ := (1+ ρ̄)/2 = ρ̃ 1/s gives a common analytic disk of radius 1/ρ∗ > 1 for (z, ϑ) 7→ A(z; ϑ) on KB . The fourfactor PAR map is smooth, so the Cauchy-estimate proof of Lemma 11 applies unchanged with this ρ∗ , the constants depending additionally on s through the polynomial degrees. 52

Scalable estimation of VARMA models

(ii) Well-specification and true-model constants. Zero-padding each factor to the fitted orders leaves the four factor polynomials, hence the expanded operators and the transfer function, unchanged, and each factor’s PAR map is a bijection onto its stable set (Proposition 1, together with the padding observation at the end of Appendix A), so the padded true coordinates ϑ0 are well defined and lie in KB by assumption. The expanded operators are the ordered products Φ̄(z) = Φ(s) (z s )Φ(z) and Θ̄(z) = Θ(s) (z s )Θ(z) of (35); the factors need not commute for d > 1, so the reversed products are in general different matrix polynomials, but the determinants agree with either ordering, det Φ̄(z) = det Φ(s) (z s ) det Φ(z), and |z| ≤ 1 ⇒ |z s | ≤ 1 shows that the expanded true operators (an ordinary VARMA whose AR and MA polynomials have degrees at most p̄ + sP̄ and q̄ + sQ̄; the leading coefficients (s) (s) −ΦP̄ Φp̄ and −ΘQ̄ Θq̄ can vanish for d > 1) satisfy Assumption 1, so the constants µmin , c0 , C0 of Lemma 12 are defined for the true seasonal process exactly as before. (iii) Everything else. The seasonal RLS and JMAP losses are (27) and (34) evaluated through the expanded Π1:m (ϑ), with the group-ridge penalty on the four unconstrained stacks, again bounded by PB on KB . The decomposition, the events, and the assemblies of Appendices C.2–C.5 use the parametrization only through Lemma 11, the dimension npar , and the true-model constants, and therefore apply verbatim. Finally, log(1/ρ∗ ) = log(1/ρ̃)/s, so the condition m ≥ log T / log(1/ρ∗ ) of Corollary 7(ii) reads m ≥ s log T / log(1/ρ̃).

Appendix D. Computational cost All timings below were measured on the runs reported in Section 7, on an Apple M4 Pro (14 cores, 48 GB). Every timed unit ran with nothing else computing on the machine; a unit is one model on the real-data panels and one (cell, seed, method) triple on the synthetic benchmark. The units run strictly one after another, and the parallelism sits inside each unit: the ridge-grid cross-validation runs one process per grid point, and the rollingwindow refits one process per window. Each figure is wall-clock seconds, so the numbers are comparable across the estimators reported here but not across hardware or implementations. The complexity statements of Section 6 are the implementation-independent counterpart, and the comparisons in Section 7 rest on forecast accuracy rather than on these timings. Table 7 reports the synthetic sweep: each entry is the mean over the eight seeds and the three training lengths of a dimension, and the standard deviations, a few tenths of a second, show how stable the per-seed cost is. Cost grows with the dimension and is nearly flat in T , which is the behaviour the Parseval evaluation of Section 6 predicts: the sufficient statistics are formed once per window, after which each optimizer iteration is independent of T . bigtime is fast on these small, short-memory problems, and its cost shows the opposite T -profile, roughly doubling from T = 800 to 3200 at d = 40, as its least-squares path is refit over all T observations. On the real-data panels the cost is reported per model: its cross-validation plus all rolling-window refits (Table 8). Dominick’s is small enough that the entire suite of ten estimators takes under three minutes. The two hourly panels are heavier, mainly through 53

Paulin and Elvira

Table 7: Synthetic benchmark, wall-clock seconds per seed (mean ± standard deviation over the eight seeds and the three training lengths T at each dimension). A unit comprises the pilot, the eight-point ridge-grid cross-validation (one process per grid point, in parallel), and the refit at the selected σ, timed with nothing else running; bigtime’s unit is its own pipeline, the penalty-path fits on the cross-validation and full windows plus the selection. dense

RLS JMAP bigtime

sparse

d=10

d=20

d=40

d=10

d=20

d=40

2.6 ± 0.4 2.6 ± 0.0 0.4 ± 0.2

8.5 ± 0.5 8.5 ± 0.2 0.6 ± 0.1

36.6 ± 0.6 37.9 ± 0.5 1.7 ± 0.4

2.2 ± 0.3 2.6 ± 0.0 0.4 ± 0.0

8.2 ± 0.2 8.6 ± 0.2 0.6 ± 0.1

36.8 ± 0.5 37.8 ± 0.7 1.7 ± 0.4

Table 8: Real-data panels, wall-clock seconds per model: its cross-validation plus all rollingwindow refits (4 windows for Beijing, 6 for Singapore, a single frozen window for Dominick’s), each model timed with nothing else running. Within a model the ridge-grid CV runs one process per grid point and the refits one process per window. The bigtime entry covers its sparse-VAR and sparse-VARMA modes together, which share one set of fits and are therefore counted once in the total. Dominick’s

Beijing

Singapore

– – 43 36 29 28 – – 3 9 – 8 5 <1

169 193 – – 144 148 67 77 5 320 4 37 1514 3

201 225 – – 200 221 108 133 6 167 5 37 2383 3

seasonal VARMA RLS seasonal VARMA JMAP VARMAX RLS VARMAX JMAP dense VARMA RLS dense VARMA JMAP seasonal VAR RLS seasonal VAR JMAP ridge VAR(X) Bayesian VAR(X) componentwise AR componentwise ARMA(X) bigtime training mean, random walk

the truncation m = 1024: every objective evaluation forms Π at 1024 lags, so a seasonal VARMA there costs a few minutes where the synthetic fits at m = 20 cost seconds. Two features of these numbers matter for practice. RLS and JMAP cost essentially the same, since they differ only in the scalar form of the loss and share the sufficient statistics, and the eight ridge-grid fits are mutually independent, so the selection parallelizes to one process per grid point. Cost grows with the dimension but is nearly flat in T , as the sufficient-statistics evaluation predicts. 54

Scalable estimation of VARMA models

The bigtime baseline is by far the most expensive model on the hourly panels: its penalty paths at lags up to 96 are fit in R (each fit single-threaded, the candidate fits and windows running concurrently), and the sparse-VARMA fit at p = 48 dominates each window.

Appendix E. Connection with state-space models VARMA models are known to have a strong connection with state-space models. As explained in (Box et al., 2015, Section 14.6), this model can be represented in an equivalent state-space form. Let r = max(p, q + 1), and let Θq+1 , . . . , Θr := 0 and Ψp+1 , . . . , Ψr := 0. As explained (Box et al., 2015, Section 14.5.2), under Assumption 1, the VARMA process also has an infinite moving average representation that can be written as ∞ X Z t = Ψ(B)at := Ψj at−j j=0

where B is the backward shift operator. As explained in (Box et al., 2015, Section 14.4.1), the matrices Ψj can be determined from the relation Φ(B)Ψ(B) = Θ(B), and satisfy the recursion Ψj = Φ1 Ψj−1 + Φ2 Ψj−2 + . . . + Φp Ψj−p − Θj ,

j = 1, 2, . . . ,

(54)

where Ψ0 = I, Ψj = 0 for j < 0, and Θj = 0 for j > q. Let Y t be an r · d dimensional (unobservable) vector called the state vector. Let us denote Y t (1) := (Y t )1:d , . . . , Y t (r) := (Y t )(r−1)d+1:rd . Let H = [I, 0, . . . , 0] be projecting to the first d components, and let Z t := HY t = Y t (1). It is explained in Box et al. (2015) that if we define the stationary VAR process Y t = F Y t−1 + Ψat , where   0 I 0 ... 0 0 0 I ... 0     .. .. . . ..  F :=  ... . .  . .   0 0 0 ... I  Φr Φr−1 . . . . . . Φ1

(55) 

and



I  Ψ1      Ψ :=  ...  ,   Ψr−1  Ψr

(56)

then Z t = Y t (1) is a VARMA process satisfying (1). The state-space transition matrix F ∈ Rrd×rd has the AR coefficients Φ1:r in its bottom block-row; it is structurally similar to the pd × pd AR companion CΦ of Section 3.2, but acts on a higher-dimensional state Y t that absorbs the MA component into a recursive form. Let Γl = E(Z t−l Z ′t ) denote the lag-l covariance matrices of the VARMA process (by definition, Γ−l = Γ′l ). As explained in Section 14.4.2, these satisfy that Γl =

Γl =

p X j=1 p X

Γl−j Φ′j −

q X

Ψj−l ΣΘ′j

for

l = 0, . . . , q,

(57)

j=l

Γl−j Φ′j

for

l > q,

j=1

with the convention that Θ0 = −I. 55

(58)

Paulin and Elvira

Appendix F. Online versions of the pointwise estimators In deployment the training window slides: at step r = 0, 1, 2, . . . the model is refit on the window Wr := {r∆ + 1, . . . , r∆ + T } of length T with stride ∆ ≤ T , standardized per (r) coordinate as Z t := (Y t − µr ) ⊘ σr with (µr , σr ) computed from Wr and ⊘ denoting coordinate-wise division. Recomputing the sufficient-statistics pre-pass from scratch at every shift costs O(T m d2 ); the exact identities below reduce the per-shift cost to O(∆ m d2 ), independent of T . Because the standardization changes with r, so that every entry of Z (r) re-used from the previous window takes a new value, the updates are maintained for the raw series Y , and the standardization is applied afterwards in closed form. The discussion covers any objective whose data-dependence is exhausted by (S 0 , S 1 , S 2 ), hence both RLS and JMAP. raw raw Raw updates and standardization. Every entry of the raw triple (S raw 0 , S1 , S2 ) (Eq. (24) with Z replaced by Y ) is a sum of outer products over the N = T − m target rows of the window, so a shift by ∆ drops the first ∆ contributions and adds ∆ new ones: with lag blocks Lj (t) := Y m−j+t , the indices aligned along the union segment of the two windows, ∆−1 X X

N +∆−1 ′ raw Y Y − Y m+t Y ′m+t , S raw W = S W + m+t r r−1 m+t 0 0 t=0

t=N raw S raw 1 (j; Wr ) = S 1 (j; Wr−1 ) +

N +∆−1 X

∆−1 X

t=N

t=0

Lj (t) Y ′m+t −

raw S raw 2 (i, j; Wr ) = S 2 (i, j; Wr−1 ) +

Lj (t) Y ′m+t ,

N +∆−1 X

∆−1 X

t=N

t=0

Li (t) Lj (t)′ −

(59)

Li (t) Lj (t)′ ,

at O(∆ m d2 ) flops and O(m d2 ) memory traffic per shift; the FFT-cached autocovariance summary of Section 6 updatesPby the same incremental sums. P Maintaining in addition the tgt lag first-moment sums M r := t<N Y m+t and M r (j) := t<N Lj (t), together with the within-window sums of Y t and Y t ⊙ Y t defining (µr , σr ) — all updated by the same dropand-add pattern — the normalized statistics of Z (r) follow in closed form by expanding the standardization; representatively,   lag ′ lag ′ ′ −1 S 2 (Z (r) ; i, j) = diag(σr )−1 S raw 2 (i, j) − M r (i) µr − µr (M r (j)) + N µr µr diag(σr ) , (60) and analogously for S 0 and S 1 , with the target moments M tgt in place of the correspondr ing lag moments. No observation in the window is re-touched; the transformation costs O(m2 d2 ), or O(m d2 ) in the FFT representation, in which S 2 is never materialized. Warm start, exactness and memory. The loss for window r is then built exactly as in Sections 4.3–5 and minimized by the same unconstrained L-BFGS, initialized at the previous window’s optimum in the unconstrained PAR coordinates; feasibility is automatic through the map. Equations (59)–(60) are exact in real arithmetic, so the incrementally maintained statistics define the same mathematical objective as a from-scratch pre-pass; floating-point results can differ slightly, the two computations summing in different orders. 56

Scalable estimation of VARMA models

Warm-starting from the preceding fit is empirically effective, typically terminating in a few iterations, but we claim no uniform iteration complexity and no equality of the returned local optima; the only statistical approximation, the use of a finite window, is the usual slidingwindow modelling choice. The persistent state is O(m d2 ) floats in the FFT representation, independent of T . The experiments of Section 7 use few, long windows and refit each one afresh; the identities above are the incremental alternative when refits are frequent.

Appendix G. Cost of the seasonal and exogenous extensions Section 6 gives the leading O(m log m d2 +m d3 ) cost of the base VARMA loss; this appendix accounts for the full per-gradient cost, including the multiplicative-seasonal factors of Section 5.5 and the exogenous regressors of the VARMAX extension. Write ℓΦ = p + sP and ℓΘ = q + sQ for the orders of the expanded operators (Φ1:ℓΦ , Θ1:ℓΘ ), and let nfft = Θ(m) be the FFT length (the smallest power of two with nfft ≥ 2m). The data enter only through a one-time pre-pass that forms the sufficient statistics; this is the sole T -dependent step. Direct summation costs O(T m d2 ) flops. The implementation instead evaluates the lagged cross-products by blocked (overlap–save) FFT correlation with length-Θ(m) blocks (distinct from the length-nfft transforms below): the forward transforms of each block are shared across the d2 channel pairs, the per-block cross-spectra are accumulated, and a single length-Θ(m) inverse transform per pair then recovers all m lags, giving O(T log m d + T d2 + m log m d2 ) = O(T log m d2 ) for the pre-pass. Each subsequent valueand-gradient evaluation is independent of T and costs


O m log m d2 + m d3 + O (p2 +P 2 +q 2 +Q2 ) d3 + O (P (p+1)+Q(q+1)) d3 | {z } | {z } | {z } (i) FFT

(ii) PAR maps

2

+ O (ℓΦ +ℓΘ ) d + | {z } (iv) assembly

3

Od | {z }

.

(iii) factor convolution

(61)

(v) Cholesky

  Term (i) forms the truncated infinite-VAR coefficients Π1:m = I −Θ(z)−1 Φ(z) 1:m by FFT inversion and evaluates the S 2 quadratic form by a Parseval (frequency-domain) identity, both at length nfft . Term (ii) is the four Heaps/Ansley–Kohn PAR maps, each a Durbin– Levinson recursion costing O(k 2 d3 ) for a factor of order k. Term (iii) convolves the regular and seasonal factors into the expanded operators (K(k+1) block products for a k × K factor pair). Term (iv) assembles and zero-pads those expanded coefficients and evaluates the per-factor penalties. Term (v) is the Cholesky log-determinant of the Σ-marginalized loss (34), absent for the RLS loss of Section 4.3. The plain VARMA is the case P = Q = 0, for which (iii) vanishes and (ii) reduces to O((p2 +q 2 )d3 ). The orders (p, q, P, Q) and the seasonal period s enter only the subdominant terms (ii)–(iv): the leading term (i) is independent of all of them, because the degree-ℓΘ expanded MA polynomial is zero-padded to length nfft before transformation, so its FFT costs the same regardless of how many of its coefficients are nonzero. A seasonal VARMA therefore evaluates at the same leading O(m log m d2 + m d3 ) per-gradient cost as a plain VARMA of equal dimension; the seasonal structure adds only the O(s(P +Q) d2 ) growth of (iv) and the O((P 2 +Q2 )d3 ) of (ii)–(iii). These subdominant terms are dominated by (i) whenever m ≳ ℓΦ (already required for the truncation bound (21) to be tight), but (61) 57

Paulin and Elvira

holds for any (p, q, P, Q, m, d), including high-order regimes (e.g. a full-diurnal regular factor p = s) where term (ii) can rival (i). Reverse-mode automatic differentiation returns the gradient at the same asymptotic cost and stores no recursion intermediates, the FFT inversion having replaced the sequential Box–Jenkins scan. The optimisation itself runs over the (p+P +q+Q)d2 unconstrained PAR coordinates, an s-fold smaller vector than the (ℓΦ +ℓΘ )d2 of the equivalent unrestricted VARMA(ℓΦ , ℓΘ ). Exogenous regressors (VARMAX). The VARMAX and seasonal-VARMAX extensions of Section 5.5 add only lower-order work. The covariates enter the sufficient statistics through the cross-covariances of Z with x, formed in the same T -dependent pre-pass, and the exogenous coefficients Πx (z) = [Θ(s) (z s ) Θ(z)]−1 Ξ(z) are obtained by passing Ξ through the same FFT inversion already used to form Π1:m . The per-gradient cost is therefore that of (61) at the same leading order, the exogenous block adding only an O(m log m d dx + m d2 dx ) FFT-and-contraction term; in particular it remains independent of T and near-linear in m.

Appendix H. Extension to VARMAX models The framework of Sections 2.1–6 and the online updates of Appendix F extend to models with eXogenous inputs (VARMAX) by a single augmentation: the exogenous regressors enter the residual, the sufficient statistics, and the Fourier evaluation alongside the lagged outputs, while the PAR map T , the Σ-marginalisation, the four (M Φ , M Θ ) ridge penalties, the Box–Jenkins recursion (10)–(11) for Π1:m , and the unconstrained L-BFGS optimisation all carry over verbatim. We record only what changes. H.1 The truncated VARMAX(p, q, s, m) model Let Y t ∈ Rd be the target series and X t ∈ Rdx a vector of exogenous covariates observed in synchrony with Y t . The VARMAX(p, q, s) model is Yt =

p X

Φj Y t−j +

j=1

s X

B l X t−l + at −

l=0

q X

Θk at−k ,

(62)

k=1

iid

with at ∼ N (0, Σ), Φ1:p ∈ Rp×d×d , Θ1:q ∈ Rq×d×d , and exogenous regression coefficients B 0:s ∈ R(s+1)×d×dx . Stationarity is again controlled by Assumption 1 on (Φ, Θ) only; the regression block B 0:s plays no role in the spectral-radius conditions. Pre-multiplying (62) by Θ(B)−1 gives the VAR(∞) form X X e l X t−l + at , Yt= Πj Y t−j + B j≥1

l≥0

in which Π1:∞ are exactly the pure-VARMA coefficients of (10)–(11) and the exogenous weights satisfy the analogous Box–Jenkins recursion min(l,q)

e l = Bl + B

X

e l−k , Θk B

e l = 0 (l < 0). B l = 0 (l > s), B

k=1

58

(63)

Scalable estimation of VARMA models

Truncating both sums at lag m gives the working model Yt ≈

m X

Πj Y t−j +

j=1

m X

e l X t−l + at . B

(64)

l=0

The endogenous truncation bound (21) applies verbatim to Π1:m (Proposition 3). The e exogenous weights inherit the same geometric coefficient decay (B(z) = Θ(z)−1 B(z) has e l ∥ ≤ C ′ ρl for every ρ > spr(CΘ )), the MA inverse as its only non-polynomial factor, so ∥B ρ but we do not claim a stationary autocovariance bound, since {X t } may be deterministic, nonstationary, or merely observed; the dropped tail instead perturbs the one-step mean by P 3 e l>m B l X t−l , negligible at the m ∼ 10 we use whenever the covariates do not grow faster −l than ρ (in particular for bounded {X t }). The regression block B 0:s enters (62) linearly and is unconstrained. H.2 Augmented sufficient statistics and residual scatter In addition to the VARMA sufficient statistics S 0 , S 1 , S 2 of Section 4.1, the VARMAX residual scatter requires three exogenous cross-statistics computed over the same N = T −m target rows: S xy 1 (l) := S yx 2 (j, l) xx S 2 (l, l′ )

:= :=

′ t X t−l Y t , P Y t−j X ′t−l , Pt ′ t X t−l X t−l′ ,

P

l = 0, . . . , m, j = 1, . . . , m,

l = 0, . . . , m,

(65)

′

l, l = 0, . . . , m.

P P e With rt := Y t − j Πj Y t−j − l B l X t−l the truncated residual of (64), the residual P (X) ′ scatter R := t rt rt extends (30) additively: R(X) = R −

X

+

X

l

j,l

e l S xy (l) − B 1

X

′e′ S xy 1 (l) B l

l X
′ yx e l S xx (l, l′ ) B e ′′ , e l S yx (j, l)′ Π′ + e +B B Πj S 2 (j, l) B l 2 l j 2 l,l′

(66)

where R is the VARMA scatter (30) unchanged. The truncated least-squares criterion is (X) yx xx Lm = tr R(X) , which depends on the data only through {S 0 , S 1 , S 2 , S xy 1 , S 2 , S 2 }. H.3 Regularised least squares and Σ-marginalised MAP The regularisation of Sections 4.3–5 carries over with one addition: fixed-variance ridge penalties on the exogenous coefficients,
iid 2 (B l )ab ∼ N 0, σB, g(a,b) ,

l = 0, . . . , s,

(67)

2 2 with two group variances σB,direct , σB,indirect separating own-channel from cross-channel regressions (tied to a single scalar, or grouped by lag, if preferred). The RLS and Σmarginalised losses are then those of Eqs. (27) and (34) with R → R(X) of (66) and the

59

Paulin and Elvira

two B-ridges appended: LRLS,X = 12 L(X) m m + (4 M -ridges) +

X S

B,g

g

Lmarg,X = ν0 +N m 2

log det Ψ0 + R

(X)

2 2σB,g

, (68)

+ (the same six ridges),

(X) 2 a,b: g(a,b)=g (B l )ab over the two groups g ∈ {direct, indirect} and Lm = Θ (p+q)d2 +(s+1)d dx , unconstrained in The L-BFGS state is (M Φ 1:p , M 1:q , B 0:s ) ∈ R

with SB,g :=

Ps

l=0

P

tr R(X) . every coordinate; the optimisation of Sections 4.3–5 is otherwise unchanged. H.4 Fourier evaluation and cost

The block-Toeplitz Parseval evaluation of Section 6 applies to each new second-moment b xx block of (65): the symmetric exogenous quadratic S xx 2 exactly as for S 2 (caching C τ = P yx ′ t X t X t+τ , τ = 0, . . . , m − 1), and the non-symmetric cross block S 2 likewise but caching P yx xy ′ b = the cross-autocovariance C τ t Y t X t+τ for both signs of τ ; the linear cross term S 1 is already O(m d dx ) and needs no acceleration. The per-evaluation cost of LRLS,X or Lmarg,X m m is then
O m log m (d + dx )2 + m (d + dx )3 , (69) b yy , reducing to (42) when dx = 0, with a single O(T log m (d + dx )2 ) pre-pass (caching C b xx , C b yx ) amortised across all L-BFGS iterations and rolling-window refits. C H.5 Online estimation The rolling-window machinery of Appendix F extends with no structural change. The raw bundle gains the three exogenous blocks S xy,raw , S yx,raw , S xx,raw , each a sum of outer prod1 2 2 ucts to which the drop-and-add identity (59) applies verbatim, at per-shift cost O(∆ m (d + dx )2 ) independent of T ; the closed-form standardisation (60) applies separately to the Y and X blocks via the same conjugate identity (60). Warm-started L-BFGS (Appendix F) is unchanged, the warm start now carrying the unconstrained regression block B 0:s with no projection, at the Fourier cost (69) per refit.

60

Related documents

• Theory of Statistical Estimation : Contributions to Higher Order Asymptotics and Non-Regular Cases / by Masafumi Akahira#1165539
• Stochastic Epidemic Model Criticism with Neural Evidence Estimation#1037246
• Local Influenza Forecasts Outperform State-Level Forecasts in the United States#343645
• Robust Semiparametric Inference for Bayesian Additive Regression Trees#625481
• Robust Semiparametric Inference for Bayesian Additive Regression Trees#623208