cultivars.multivariate.reduced_form.moving_average#
The vector autoregressive moving average, estimated without a likelihood search.
A VARMA(p, q),
buys parsimony with a harder estimation problem: the moving-average innovations are unobserved, so least squares cannot be applied directly, and the Gaussian likelihood is a surface over \(k^2 (p + q)\) parameters that is neither cheap nor – for an unrestricted model – uniquely maximized. This module takes the Hannan-Rissanen route instead: recover the innovations from a long autoregression of order \(h \sim \sqrt{T}\), then estimate the VARMA by regressions in which those estimated innovations stand in for the true ones, and once more after re-deriving them through the model. Three passes, every one of them a least-squares solve, no numerical optimizer anywhere.
Two commitments shape the surface. First, the estimator’s conditioning is stated rather than hidden. The standard errors treat the lagged innovations in the design as observed regressors, so they are conditional in the Hannan-Rissanen sense and modestly understate uncertainty; the summary says so on every result. The third pass runs only when the second-stage moving-average block is invertible – re-deriving innovations through a non-invertible block would amplify the approximation error instead of reducing it – and the result records whether it ran. Second, the identification caveat is structural rather than numerical and is repeated on the result: distinct \((A, M)\) pairs can generate identical second moments, and the standard resolution – echelon-form restrictions indexed by Kronecker indices – is deliberately not imposed here. What survives the ambiguity is exactly what most users come for: for a stable, invertible representation the moving-average matrices \(\Psi_h = M_h + \sum_i A_i \Psi_{h-i}\), and with them the forecasts, impulse responses, and variance decompositions, are invariant across observationally equivalent parameterizations. Overriding the moving-average recursion is the single change that makes the inherited propagation surface correct for a VARMA, because every member of that surface reads \(\Psi_h\) rather than the autoregressive companion.
Layout. VARMA extends _VectorAutoRegressionModel in
_internals with the moving-average order, a long-order rule and a
three-pass fit() built from the base’s _least_squares and
_gaussian_moments: _stage_regression lays out
[deterministic | own lags | innovation lags] with lag_matrix and
deterministic_columns from _core, _unpack splits a pass’s
coefficients into the three blocks, and _recursive_innovations rebuilds
the innovations for the refinement; lag_order_selection is refused
because (p, q) candidates are not nested designs. VARMAResult
is a full _VectorResult with the summary, comparison, inference and
propagation mixins, overriding ma_representation and forecast and
adding the invertibility check; the identification and conditioning notes
it prints are _VARMA_IDENTIFICATION_NOTE and _HR_CONDITIONAL_NOTE in
_core. The pure autoregression is
vector_autoregression; the
scalar model with exact maximum likelihood is
box_jenkins.
References
Hannan, E. J., & Rissanen, J. (1982). Recursive estimation of mixed autoregressive-moving average order. Biometrika, 69(1), 81-94.
Hannan, E. J., & Kavalieris, L. (1984). Multivariate linear time series models. Advances in Applied Probability, 16(3), 492-561.
Dufour, J.-M., & Jouini, T. (2005). Asymptotic distribution of a simple linear estimator for VARMA models in echelon form. In Statistical Modeling and Analysis for Complex Data Problems (pp. 209-240). Springer.
Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer. Chapters 11-12.
Example
A bivariate VARMA(1, 1) whose moving-average block carries a spillover from the first innovation into the second equation. The three passes recover both stacks, a VAR(1) on the same data cannot represent the moving-average part and leaves autocorrelated residuals, and the one-step forecast carries the last innovation across the sample boundary:
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> A = np.array([[0.5, 0.1], [0.0, 0.4]])
>>> M = np.array([[0.4, 0.0], [0.2, 0.3]])
>>> u = rng.standard_normal((402, 2))
>>> y = np.zeros((402, 2))
>>> for t in range(1, 402):
... y[t] = A @ y[t - 1] + u[t] + M @ u[t - 1]
>>> res = VARMA(y[2:], order=(1, 1)).fit()
>>> bool(np.abs(res.coefficients[0] - A).max() < 0.15)
True
>>> bool(np.abs(res.ma_coefficients[0] - M).max() < 0.15)
True
>>> res.refined, res.is_stable, res.is_invertible
(True, True, True)
>>> from cultivars.multivariate.reduced_form.vector_autoregression import VAR
>>> plain = VAR(y[2:], order=1).fit()
>>> lag1 = lambda e: np.corrcoef(e[1:, 0], e[:-1, 0])[0, 1]
>>> bool(abs(lag1(plain.resid)) > 2 * abs(lag1(res.resid)))
True
>>> res.forecast(2).shape
(2, 2)
Classes#
Vector autoregressive moving average, estimated by Hannan-Rissanen. |
|
A fitted vector autoregressive moving average. |