Making VARMA Usable Again: Fixed-Size Statistics, Guaranteed Stability, No T in the Inner Loop

Paper: 2608.06340 Authors: Daniel Paulin, Victor Elvira Categories: stat.ML, cs.LG, math.ST

The Gap

VARMA is the textbook-correct model for multivariate time series and almost nobody uses it. The reason is not ignorance, it is three concrete failures that compound.

First, cost per iteration. Exact Gaussian likelihood for a VARMA goes through a Kalman filter or the innovations algorithm: one full pass over all $T$ observations, at $O(T d^3)$, for every single parameter proposal. Optimizers need hundreds to thousands of those.

Second, the feasible set is horrible. Stationarity needs the AR polynomial’s roots outside the unit circle; invertibility needs the same of the MA polynomial. That region is non-convex, has no cheap projection, and its boundary is exactly where the likelihood likes to wander. Conditional MLE routinely walks out and returns a non-invertible fit. A non-invertible MA is not a slightly worse model, it is a broken one: the recursion that reconstructs innovations from observations diverges, so the forecasts diverge with it.

Third, identifiability. VARMA parameters are identified only up to an equivalence class, and the classical fix — echelon forms, Kronecker indices, scalar component models — is bookkeeping that practitioners refuse to maintain.

So the field routed around the problem. Fit a VAR with many lags (Lütkepohl’s default; Hannan–Rissanen’s two-stage version makes it cheap) and accept that you now estimate $p d^2$ parameters to imitate what an MA term does with $d^2$. Or put a Minnesota prior on it and call it a Bayesian VAR. Or fit $d$ univariate ARMAs and give up cross-series dynamics. Or impose sparsity on VARMA coefficients and lose the likelihood’s information along the way. Every one of these is a concession that the MA term, which is where the parsimony lives, is not worth the fight.

PROBLEM: VARMA is the parsimonious model, unusable past d ~ 5
   |
   +-- prior art .................. what it gives up
   |     VAR(p), long lags .......... p*d^2 params, no MA parsimony
   |     conditional / exact MLE .... O(T d^3) per iter, escapes to
   |                                  non-invertible fits
   |     echelon-form MLE ........... identification bookkeeping
   |     component-wise ARMA ........ no cross-series dynamics
   |     sparse VARMA ............... convex-ish, discards likelihood
   v
ASSUMPTION: the estimand is the VAR(oo) predictive representation,
            not the VARMA parameters themselves
   |        (so non-identifiability stops mattering)
   v
METHOD: [PACF reparam: stable by construction]
      + [Gaussian prior, separate diag / off-diag scales]
      + [loss touches data only via G_hat(0..L), Parseval/FFT]
   |
   v
EVIDENCE: near-parametric consistency for the VAR(oo) fit (fixed d)
        + d = 10..40 simulations stay near oracle forecast error
        + retail / meteo / air-quality: >= VAR, BVAR, ARMA, sparse-VARMA
   |
   v
CONCLUSION: per-iteration cost free of T, fits always invertible,
            VARMA lands in the regime where people defaulted to VAR

The Increment

One sentence: Before, fitting a VARMA meant a pass over the data per optimizer step and a real chance of a divergent fit; after, the data is compressed once into a fixed-size block of autocovariances and every step is a stable-by-construction algebra problem whose cost never sees $T$ again.

Core Mechanism

The framework has two clearly separated halves, and the separation is the whole idea. The data half runs once. Sample autocovariance matrices $\hat\{G\}(0), \hat\{G\}(1), \ldots, \hat\{G\}(L)$ up to a truncation lag $L$ are computed via FFT at near-linear cost in $T$. That block is a fixed-size object: $(L+1)d^2$ numbers regardless of whether $T$ is a thousand or ten million. Nothing downstream ever looks at $y_t$ again.

The model half runs every iteration and never touches $T$. Unconstrained parameter matrices are pushed through a partial-autocorrelation map — the multivariate generalization of the Ansley–Kohn / Barndorff-Nielsen–Schou transform — which is a bijection onto the stationary region for the AR side and onto the invertible region for the MA side. Every point in the unconstrained space lands on a valid model. The optimizer literally cannot produce a non-invertible fit. From those polynomials the framework computes the truncated $VAR(\infty)$ coefficients $\Pi_1, \ldots, \Pi_L$ implied by $\Pi(z) = I - B(z)^\{-1\}A(z)$, and it does so in the frequency domain: evaluate the transfer function on $L$ Fourier frequencies, invert a $d \times d$ matrix at each, inverse-FFT back. Near-linear in $L$.

Then Parseval closes the loop. The least-squares criterion for a linear predictor is a quadratic form in the autocovariances and the coefficients, and Parseval’s identity says that lag-domain sum equals a frequency-domain integral, so it can be assembled at the same near-linear cost. That gives the implied residual covariance $S(\theta)$, from which two estimators fall out: a regularized least-squares fit that scores $S(\theta)$ against a weight matrix plus the Gaussian penalty, and a covariance-marginalized MAP that integrates the innovation covariance out under a conjugate prior, producing a $\log\det(S(\theta) + \Psi)$ objective that self-regularizes against degenerate covariance estimates. The priors are Minnesota-flavored: one scale for diagonal entries (a series’ own history), a tighter one for off-diagonal (everyone else’s).

  ONCE, cost ~ O(T log T)             EVERY ITERATION, cost free of T
  -------------------------           --------------------------------
  y_1 .. y_T                          theta (unconstrained matrices)
      |                                        |
      v  FFT                                   v  PACF map
  G_hat(0), G_hat(1), .. G_hat(L)      A(z), B(z), roots outside the
      |                                 unit circle BY CONSTRUCTION
      |  fixed size: (L+1) d^2                  |
      |                                         v  FFT over L freqs
      |                                 Pi(z) = I - B(z)^\{-1\} A(z)
      |                                         |  inverse FFT
      |                                         v
      |                                 Pi_1 .. Pi_L  (VAR(oo) trunc.)
      |                                         |
      +------> [ Parseval quadratic form ] <----+
                          |
                          v
                  S(theta) = implied residual covariance
                          |
            +-------------+--------------+
            v                            v
   RLS: tr(W S(theta))          MAP: log det(S(theta) + Psi)
            \                            /
             +------------+-------------+
                          |  + N(0, tau_diag) / N(0, tau_off) penalty
                          v
                  autodiff gradient step, repeat

The metaphor: a tailor who casts a mannequin. The old way is a tailor who calls the customer back in for every pin: try the jacket, mark it, customer goes home, tailor adjusts, customer returns. Hundreds of fittings, and each fitting costs a whole visit — that is $O(T)$ per iteration.

This paper’s tailor takes one careful set of measurements and casts a plaster mannequin. That cast is $\hat\{G\}(0..L)$: it is not the customer, but it contains everything about the customer that a jacket can respond to. Every subsequent fitting happens against the mannequin, at a cost that has nothing to do with how far away the customer lives.

The pattern-cutting rules are the PACF reparametrization: this tailor works from a block pattern whose adjustments are jointed like a ball-and-socket, so no combination of dial settings can produce a garment that a human body cannot enter. The old tailor could — and regularly did — pin a shape that was geometrically impossible to wear, which is what a non-invertible fit is. The house block itself is the prior: absent evidence, drift toward a conventional cut, and be more conservative about the unusual seams (off-diagonal, cross-series) than the obvious ones (diagonal, own history). And Parseval is the tailor’s trick of putting the mannequin on a turntable and taking one sweep of measurements as it rotates, instead of measuring every seam separately — a geometric identity guarantees the two totals agree, and the sweep is far faster.

Key Concepts

• Sufficient statistics for a linear predictor: To score how well a linear forecast does, you never need the data — only its second moments. Concretely: predict tomorrow’s temperature as a weighted sum of the last three days. The expected squared error is a quadratic function of those three weights whose coefficients are just the covariances among four numbers. Ten years of data or ten thousand, once you have the covariance table the scoring is identical work. A truncated $VAR(\infty)$ loss is exactly this, at scale, which is why $\hat\{G\}(0..L)$ is all the data the optimizer needs.

• Reparametrizing a nasty constraint set into a box: Scalar intuition first. An AR(2) is stationary only when $(\phi_1, \phi_2)$ lies inside a particular triangle — an awkward region to optimize inside or sample from. But there is a smooth invertible map from the plain square $(-1,1)^2$ of *partial autocorrelations onto that triangle. Optimize freely in the square; every point you visit is stationary. The multivariate version replaces “absolute value below one” with “singular values below one” on matrices, and the paper applies it to the MA polynomial too, which buys invertibility. The payoff is not just convenience: it removes an entire failure mode, because the optimizer has nowhere bad to go.

• Why non-invertibility destroys forecasts: An MA term says today’s observation carries yesterday’s shock. To forecast you must run that backwards — infer the unobserved shocks from the observed series. If the MA polynomial’s roots sit inside the unit circle, that backwards recursion multiplies errors by something bigger than one at each step, so the inferred shocks blow up and the forecast blows up with them. This is why the paper’s claim that classical conditional MLE “returns non-invertible fits whose forecasts diverge” at $d = 40$ is not a rounding-error complaint. The fit is unusable, not merely suboptimal.

• Marginalizing the nuisance covariance instead of plugging it in: Classical practice profiles out the innovation covariance $\Sigma$ by substituting its conditional MLE, which in high dimension is badly conditioned or singular. Integrating $\Sigma$ out under a conjugate prior instead turns the objective into a $\log\det$ of a regularized residual covariance. Same computational shape, but the estimator now pays a price for pretending some direction has zero noise.

Framework Shift

Before (mainstream approach):          After (this paper):

y_1 .. y_T                             y_1 .. y_T
    |   re-read EVERY iteration             |   read ONCE
    v                                       v
[Kalman filter / innovations]          [ G_hat(0..L) ]  fixed size
    |   O(T d^3) per step                   |
    v                                       |   ( no T below this line )
[ likelihood(theta) ]                       v
    |                                  [theta] --PACF--> stable A(z),B(z)
    v                                       |
[constrained optimizer]                     v  FFT
    |  boundary is where                [ Pi_1 .. Pi_L ]
    |  the optimum hides                    |
    v                                       v  Parseval
[fit may be non-invertible]            [ loss(theta) ]
    |                                       |  fit always invertible
    v                                       v
[forecasts diverge at d=40]            [near-oracle forecasts, d=10..40]

    data and model are FUSED               data and model are SEPARATED
    (each score needs the series)          (score needs only the cast)

From re-reading the series to interrogating a fixed-size cast of it, the core shift is moving the data out of the optimization loop and moving the constraints out of the optimizer and into the parametrization.

Expert Assessment

Problem choice: A real gap, and an old one. VARMA has been “theoretically superior, practically abandoned” since the 1980s, and the abandonment was driven by exactly the three obstacles named here. Choosing to attack computation *and the feasible set and identifiability simultaneously is the right diagnosis — fixing any one alone leaves the model unusable. One honest caveat about the framing: “per-iteration cost independent of $T$” is the headline, but for retail and macro data $T$ is a few hundred and was never the binding constraint. The claim earns its keep in two narrower places the abstract mentions almost in passing: rolling-window refits, where you re-run the fit hundreds of times, and long high-frequency series. The thing practitioners will actually care about is the guaranteed-invertible fit at $d = 40$, and that comes from the reparametrization, not from the $T$-independence.

Method maturity: Clever assembly rather than a new primitive. The PACF stationarity map is Ansley–Kohn, recently revived for Bayesian VARs by Heaps; the frequency-domain evaluation of a quadratic form via Parseval is close kin to Whittle estimation, which has been $T$-independent-per-iteration for decades once you precompute the periodogram; the $VAR(\infty)$ target and the sample-autocovariance route to it are Hannan–Rissanen. The contribution is that the pieces fit together and the composite has theory. That is legitimate, but it puts a burden on the paper: it must state precisely what it buys over (a) a Whittle-likelihood VARMA on the same reparametrization, and (b) plain two-stage Hannan–Rissanen, which is cheap, consistent, and the obvious “simpler approach being overlooked.” From the abstract I cannot tell whether that comparison is made head-on, and it is the first thing I would look for in the paper. Second reservation: the theory is stated for fixed dimension, while the selling point is $d = 10$ to $40$. The regime the experiments live in is not the regime the theorems cover. That is a common and forgivable gap, but it should be said out loud rather than left for the reader to notice.

Experimental integrity: The baseline set is the right one — VAR, Bayesian VAR, component-wise ARMA, sparse VARMA — and the phrasing “match or beat” on three real datasets is refreshingly unoversold; papers that claim uniform wins on real forecasting data are usually hiding something. Two things to watch. The simulation comparison against conditional MLE is close to a straw man: the fact that an unconstrained optimizer leaves the invertible region is precisely the problem this method solves by construction, so of course it wins, and a reader should ask how the method fares against a *repaired CMLE (root-flipping, or CMLE run on the same PACF parametrization) to separate “we fixed the constraint handling” from “our loss is better.” And “stays close to the oracle” in simulation is generous by construction when the data-generating process is itself a VARMA — the interesting question is misspecification, which the abstract does not mention. I would also want $T$ reported alongside $d$; the difficulty of $d = 40$ depends entirely on it.

Writing quality: The abstract is unusually good — it names the obstacles, the ingredients, the theory, and the empirical scope without padding, which is rarer than it should be. The corner I suspect was cut is identifiability. The framework does not *solve non-identifiability, it sidesteps it by declaring the estimand to be the $VAR(\infty)$ predictive representation, which is identified even when the VARMA parameters are not. That is a good move and arguably the paper’s sharpest conceptual choice, but the abstract lists identifiability as an obstacle and then never explains that the answer is “we changed the target.” Rewriting that discussion, and placing the method explicitly on the map next to Whittle estimation and Hannan–Rissanen, is the single edit that would raise the whole paper a level.

Verdict: weak accept — a well-engineered, theoretically backed synthesis that genuinely makes an abandoned model class usable, held back from a stronger rating by theory that does not cover the advertised regime and an unclear delta against Whittle-type and Hannan–Rissanen estimators.

Takeaways

Four things transfer well beyond VARMA:

Split the loss into a data compression pass and a model pass. Any objective that touches the data only through second moments — linear predictors, Gaussian likelihoods, stationary-kernel GPs, linear state-space models — can be rewritten so the optimizer sees a fixed-size statistic instead of the dataset. Ask of any iterative fit you own: what is the smallest object the loss actually needs, and can I compute it once? Often it is much smaller than the data and nobody checked.

Prefer reparametrization to constrained optimization. If your feasible set is defined by a spectral or root condition — stable linear RNNs and SSM layers, controller synthesis, covariance structure, copula parameters — look for a smooth bijection from an unconstrained box onto the feasible set before you reach for projections, penalties, or rejection sampling. You get free autodiff, no boundary pathologies, and, importantly, you eliminate a failure mode rather than mitigating it. The PACF construction specifically is a reusable tool: it handles “all roots outside the unit circle” for matrix polynomials.

When parameters are unidentified, change the estimand to something that is. Rather than imposing canonical forms to make parameters unique, target the identified functional you actually care about — here the predictive filter. This framing works anywhere over-parametrized models meet inference: you often do not need parameter identifiability, you need identifiability of the thing you will use.

Integrate out nuisance covariances instead of plugging in their MLEs. The resulting $\log\det$ objectives cost the same and stop the estimator from claiming a direction has zero noise. Pair it with structured shrinkage: separate prior scales for diagonal and off-diagonal blocks is a one-line change that encodes “a series’ own past is more informative than its neighbors’,” and that asymmetry is real in almost every multivariate system.

One caution to steal along with the ideas: precomputed sufficient statistics are only sufficient for the loss you committed to. The moment you want a robust loss, a non-Gaussian likelihood, or heteroskedastic weighting, the compression may no longer be lossless. That constraint is the price of the speedup, and it is worth being explicit about before adopting the pattern.