Iterated Generalized Yule-Walker Estimator

• The iterated generalized Yule–Walker estimator is a moment-based method that replaces classical linear equations with iterative fixed-point or alternating updates to solve nonlinear/bilinear moment systems.
• It effectively handles endogeneity and nonlinearity in settings like periodic covariance extension and high-dimensional matrix-valued spatio-temporal autoregressions.
• The approach uses techniques such as spectral factorization, alternating least squares, and quasi-Newton iterations to ensure consistent parameter identification.

An iterated generalized Yule–Walker estimator is a moment-based estimator that replaces the classical linear Yule–Walker system by generalized covariance equations and then solves the resulting nonlinear or bilinear system by repeated updates of the unknown parameters. In the literature considered here, this label is explicit in two rather different settings: periodic or circulant rational covariance extension, where a nonlinear Yule–Walker equation for a denominator spectral factor is solved by a fixed-point or quasi-Newton recursion, and high-dimensional matrix-valued spatio-temporal autoregressions, where Yule–Walker or GMM equations are solved by alternating updates of row-side and column-side coefficient matrices (Picci et al., 2016, Dou et al., 14 Aug 2025). Related work uses generalized Yule–Walker equations without defining an iterated estimator, including rowwise least-squares estimation for heterogeneous spatio-temporal models, sparse penalized structural estimation, and regularized operator-valued estimation in Hilbert space (Dou et al., 2015, Reuvers et al., 2021, Niu et al., 6 Jun 2025).

1. Conceptual scope

Classical Yule–Walker estimation is based on linear autocovariance identities for scalar AR or standard VAR models. The generalized variants discussed in this literature depart from that setting in several ways. In structural spatio-temporal models, contemporaneous feedback creates endogeneity because the dependent variable appears on both sides of the equation, so ordinary least squares is inconsistent; generalized Yule–Walker equations use lagged covariance identities to identify structural coefficients instead (Dou et al., 2015, Reuvers et al., 2021). In periodic covariance extension, the unknown enters nonlinearly through the ratio $b/a$, so the Yule–Walker relation itself becomes nonlinear (Picci et al., 2016). In matrix-valued autoregressions, the coefficient is constrained by Kronecker or bilinear structure, so the moment equations are no longer an unconstrained linear system in a single coefficient matrix (Dou et al., 14 Aug 2025, Kołodziejski, 21 May 2025). In functional autoregression, the covariance object is a compact operator and direct inversion is ill-posed, so Yule–Walker estimation becomes a regularized inverse problem (Niu et al., 6 Jun 2025).

Within this broader class, the adjective “iterated” has a narrower meaning. It does not refer merely to numerical optimization. In the periodic covariance-extension paper, the estimator is explicitly a fixed-point recursion for the AR denominator polynomial, with a scaling update and a spectral-factor projection step (Picci et al., 2016). In the matrix-valued spatio-temporal paper, the estimator is explicit alternating re-estimation: given row-side parameters, update column-side parameters; given column-side parameters, update row-side parameters; then repeat until the identifiable Kronecker products stabilize (Dou et al., 14 Aug 2025). By contrast, several nearby papers are generalized Yule–Walker in substance but not iterative in this econometric sense (Dou et al., 2015, Reuvers et al., 2021, Niu et al., 6 Jun 2025).

2. Generalized Yule–Walker moment systems

The generalized Yule–Walker idea is to derive moment restrictions from lagged second moments that remain valid even when direct regression is invalid or structurally misspecified. For the structural spatio-temporal autoregression

$y_t = A y_t + B y_{t-1} + \varepsilon_t,$

post-multiplication by $y_{t-1}'$ and expectation yield

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$

or equivalently

$$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$

This is “generalized” because the model is structural, contains contemporaneous spatial interactions $Ay_t$, and the moment equations identify $A$ and $B$ rather than only reduced-form autoregressive coefficients (Reuvers et al., 2021).

For the heterogeneous diagonal-coefficient spatio-temporal model

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$

the matrix Yule–Walker equation becomes

$$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$

and the $y_t = A y_t + B y_{t-1} + \varepsilon_t,$0-th row yields an overidentified system in $y_t = A y_t + B y_{t-1} + \varepsilon_t,$1, estimated by least squares on sample covariance substitutions (Dou et al., 2015).

In the periodic covariance-extension problem, the generalized Yule–Walker equation is nonlinear: $y_t = A y_t + B y_{t-1} + \varepsilon_t,$2 where $y_t = A y_t + B y_{t-1} + \varepsilon_t,$3 is built from the inverse DFT of $y_t = A y_t + B y_{t-1} + \varepsilon_t,$4. The nonlinearity is the defining feature: unlike classical Yule–Walker, the right-hand side depends on the unknown $y_t = A y_t + B y_{t-1} + \varepsilon_t,$5 through $y_t = A y_t + B y_{t-1} + \varepsilon_t,$6 (Picci et al., 2016).

In the matrix-valued spatio-temporal setting,

$y_t = A y_t + B y_{t-1} + \varepsilon_t,$7

the Yule–Walker identities are blockwise: $y_t = A y_t + B y_{t-1} + \varepsilon_t,$8 These are moment equations of the form $y_t = A y_t + B y_{t-1} + \varepsilon_t,$9, and the paper explicitly characterizes the resulting estimator as GMM-type (Dou et al., 14 Aug 2025).

For ARH($y_{t-1}'$0) in Hilbert space,

$y_{t-1}'$1

the Yule–Walker equations are operator-valued: $y_{t-1}'$2 with $y_{t-1}'$3 and $y_{t-1}'$4. The compactness of $y_{t-1}'$5 makes $y_{t-1}'$6 unbounded, so generalized Yule–Walker estimation requires regularization rather than direct inversion (Niu et al., 6 Jun 2025).

3. Fixed-point and quasi-Newton iteration in periodic covariance extension

The clearest explicit form of an iterated generalized Yule–Walker estimator appears in the periodic or circulant rational covariance-extension problem. The data are finite covariance lags $y_{t-1}'$7, and the target is a rational discrete spectral density

$y_{t-1}'$8

with $y_{t-1}'$9 fixed and $\Sigma_1 = A\Sigma_1 + B\Sigma_0,$0 to be estimated through its outer spectral factor $\Sigma_1 = A\Sigma_1 + B\Sigma_0,$1 (Picci et al., 2016).

The central nonlinear Yule–Walker system is

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$2

where $\Sigma_1 = A\Sigma_1 + B\Sigma_0,$3 is the periodic impulse response determined by

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$4

Solving formally for $\Sigma_1 = A\Sigma_1 + B\Sigma_0,$5 yields the iteration

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$6

with

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$7

This is the paper’s explicit iterative generalized Yule–Walker scheme (Picci et al., 2016).

The actual algorithm is a projected iteration. After the update, the method performs spectral factorization of

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$8

and replaces the result by the outer spectral factor, because the raw recursion need not preserve the Schur constraint. Initialization is suggested by the maximum entropy or Levinson solution, and the stopping rule is

$\Sigma_1 = A\Sigma_1 + B\Sigma_0,$9

Computationally, $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$0 is obtained by evaluating $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$1 on $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$2, applying inverse FFT, and then building the Toeplitz matrix (Picci et al., 2016).

The paper’s conceptual contribution is the link between this iteration and the Lindquist–Picci generalized entropy formulation. Writing the dual functional in spectral-factor coordinates,

$$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$3

the generalized Yule–Walker residual is half the gradient: $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$4 Under $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$5, the update becomes

$$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$6

so the iteration is a quasi-Newton method with constant approximate Hessian $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$7 (Picci et al., 2016).

The convergence result is local rather than global. Under feasibility, the minimizer is unique, $$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$8 is locally strictly convex near the optimum, and the quasi-Newton descent converges locally to the unique stationary point. At convergence, the limiting ARMA model matches the prescribed covariance data. The same framework is then used for finite-interval smoothing, where the periodic ARMA approximation produces a block-circulant banded linear system and yields a simpler procedure than classical Riccati-based calculations (Picci et al., 2016).

4. Alternating generalized Yule–Walker iteration for matrix-valued spatio-temporal autoregressions

A second explicit form of iterated generalized Yule–Walker estimation is developed for high-dimensional matrix-valued time series

$$\Sigma_1' = \begin{bmatrix}\Sigma_1' & \Sigma_0\end{bmatrix}\begin{bmatrix}A & B\end{bmatrix}'.$$9

with the central model

$Ay_t$0

Here both row and column interaction matrices are unknown and banded, and endogeneity arises because $Ay_t$1 appears on the right-hand side together with the innovation $Ay_t$2 (Dou et al., 14 Aug 2025).

The generalized Yule–Walker backbone comes from multiplying column equations by lagged columns and taking expectations. For $Ay_t$3,

$Ay_t$4

with sample analogues

$Ay_t$5

Because the moment equations are bilinear in row-side and column-side parameters, the paper defines the estimator through alternating conditional least-squares or GMM updates rather than a single closed-form solution (Dou et al., 14 Aug 2025).

Conditional on $Ay_t$6, each row $Ay_t$7 of $Ay_t$8 is updated by a band-restricted least-squares regression,

$Ay_t$9

where only support coordinates allowed by the bandedness assumption are estimated. Conditional on $A$0, each row $A$1 of $A$2 is updated analogously,

$A$3

The algorithm alternates these two blocks, normalizes $A$4 and $A$5 to satisfy $A$6 and $A$7, and stops when

$A$8

is close to zero. The identifiable objects are the Kronecker products, and the stopping rule is defined accordingly (Dou et al., 14 Aug 2025).

Initialization is based on a nearest Kronecker product procedure. The model is first vectorized,

$A$9

and candidate $B$0 are obtained from an existing generalized Yule–Walker procedure for vector spatio-temporal models. The factors $B$1 and $B$2 are then extracted from nearest-Kronecker-product problems solved by singular value decomposition. The paper treats NKP chiefly as an initialization device and reports that iteration improves finite-sample accuracy relative to NKP alone (Dou et al., 14 Aug 2025).

A distinctive feature is dual-bandwidth estimation. Because there are separate unknown bandwidths $B$3 and $B$4, the paper introduces a two-step, ratio-based sequential procedure built from directional drops in residual sums of squares. First $B$5 is estimated by maximizing a ratio involving

$B$6

aggregated over coordinates; then $B$7 is estimated conditionally on $B$8. Under mixing, boundedness, identifiability, minimum edge-signal conditions, and projection-matrix regularity, the paper proves bandwidth selection consistency: $B$9 It also gives consistency and asymptotic normality for the generalized Yule–Walker estimator under growing dimensions (Dou et al., 14 Aug 2025).

Empirically, the iterative estimator outperforms the NKP estimator in estimating $$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$0 and $$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$1, both when bandwidths are known and when they are estimated. In the trading-volume application, the banded Yule–Walker model delivers the best forecasting and execution performance among the compared methods (Dou et al., 14 Aug 2025).

5. Relation to non-iterated generalized Yule–Walker estimators

The term “iterated generalized Yule–Walker estimator” should not be extended indiscriminately to all work built on generalized Yule–Walker equations. Several influential papers are directly relevant to generalized Yule–Walker estimation but do not define an iterated estimator.

The rowwise spatio-temporal model with unknown diagonal coefficients estimates, for each location $$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$2,

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$3

This is a one-shot least-squares solution to an overidentified Yule–Walker system. The paper also proposes a reduced-equation refinement based on informative-moment screening, but it does not introduce initialization, fixed-point updating, or repeated re-estimation until convergence (Dou et al., 2015).

The sparse structural spatio-temporal paper develops SPLASH, a penalized generalized Yule–Walker estimator based on

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$4

and the convex criterion

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$5

The paper explicitly states that it does not define an “iterated generalized Yule–Walker estimator” as an econometric object. The iteration appears only at the level of numerical optimization for the sparse-group-lasso-type problem (Reuvers et al., 2021).

For matrix autoregressive models,

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$6

a generalized Yule–Walker estimator is defined by the constrained Frobenius minimization

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$7

Its implementation uses L-BFGS-B. The same paper contains an explicitly iterative Burg method, but that procedure is not Yule–Walker (Kołodziejski, 21 May 2025).

In functional time series, the operator-valued estimator

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$8

is a regularized generalized Yule–Walker estimator. The paper mentions “iterative methods” only for solving kernel integral equations numerically; it does not introduce a recursive statistical estimator based on repeated Yule–Walker updating (Niu et al., 6 Jun 2025).

A further variant arises in the presence of mean-shift changepoints. There the proposed estimator is

$$\mathbf y_t = D(\boldsymbol\lambda_0)W\mathbf y_t + D(\boldsymbol\lambda_1)\mathbf y_{t-1} + D(\boldsymbol\lambda_2)W\mathbf y_{t-1} + \boldsymbol\varepsilon_t,$$9

where $$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$0 is built from autocorrelations of first differences. This is a modified or generalized Yule–Walker estimator tailored to changepoint contamination, but not an iterative one (Gallagher et al., 2021).

6. Identification, asymptotics, and numerical architecture

Across these literatures, generalized Yule–Walker estimation is driven by identification through covariance structure rather than direct regression. In structural spatio-temporal models, lagged variables eliminate contemporaneous innovation correlation; this is the basic device for overcoming endogeneity (Dou et al., 2015, Reuvers et al., 2021, Dou et al., 14 Aug 2025). In periodic covariance extension, identification is tied to feasibility, Schur factorization, and the bijection between spectral factors and positive denominator polynomials (Picci et al., 2016). In matrix-valued models, normalization such as $$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$1 is required because $$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$2, so only the Kronecker products are directly identifiable (Dou et al., 14 Aug 2025, Kołodziejski, 21 May 2025). In Hilbert space, identifiability relies on a Picard condition and injectivity of the covariance operator $$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$3 (Niu et al., 6 Jun 2025).

The asymptotic targets also differ. The periodic covariance-extension paper proves local convergence of the quasi-Newton or fixed-point algorithm to the unique minimizer of the generalized entropy dual, not a global convergence theorem (Picci et al., 2016). The matrix-valued spatio-temporal paper proves consistency and asymptotic normality for the generalized Yule–Walker estimator under growing dimensions, together with consistency of the sequential bandwidth selector (Dou et al., 14 Aug 2025). The functional paper proves consistency in Hilbert–Schmidt norm for the regularized estimator and asymptotic normality for predictors, but also shows that the autoregressive-operator estimator does not converge weakly in distribution under the operator norm topology because of the compactness of the autocovariance operator (Niu et al., 6 Jun 2025). The sparse structural paper derives finite-sample error bounds and estimation consistency for a penalized generalized Yule–Walker criterion built on banded covariance estimation and structured sparsity (Reuvers et al., 2021).

Numerically, the class is heterogeneous. The periodic estimator relies on FFT-based evaluation over the discrete unit circle, Toeplitz and circulant structure, and a projection step onto the outer spectral factor (Picci et al., 2016). The matrix-valued spatio-temporal estimator decomposes a large bilinear problem into many small rowwise regressions subject to hard banding, with NKP used only for initialization (Dou et al., 14 Aug 2025). SPLASH reduces high-dimensional structural estimation to sparse-group-lasso computation over a regularization path (Reuvers et al., 2021). Functional Yule–Walker estimation combines empirical covariance operators, truncation, and Tikhonov filtering

$$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$4

in an eigensystem expansion (Niu et al., 6 Jun 2025). The MAR Yule–Walker estimator uses nonlinear constrained optimization and is reported to become slow in high dimension (Kołodziejski, 21 May 2025).

7. Applications, interpretation, and limits of the term

The applications illustrate that “iterated generalized Yule–Walker” is not tied to a single disciplinary domain. In periodic covariance extension, the estimator supports finite-interval smoothing for periodically extended processes and is presented as providing a simpler procedure than classical Riccati-based calculations (Picci et al., 2016). In matrix-valued spatio-temporal analysis, the method is motivated by intraday trading volume curves, where rows are synchronized time buckets and columns are assets; the empirical results favor estimating both row and column dependence structures rather than fixing them ex ante (Dou et al., 14 Aug 2025).

The surrounding generalized Yule–Walker literature broadens the interpretation. Sparse structural estimation has been applied to satellite-measured $$\big(I-D(\boldsymbol\lambda_0)W\big)\Sigma_1 = \big(D(\boldsymbol\lambda_1)+D(\boldsymbol\lambda_2)W\big)\Sigma_0,$$5 concentrations in London, where forecast improvements and strong spatial interactions were reported, but that method is penalized rather than iterated (Reuvers et al., 2021). Functional Yule–Walker estimation has been applied to high-frequency wearable sensor data through an ARH(5) model tuned by 5-fold cross-validation, again without an iterated Yule–Walker estimator in the strict sense (Niu et al., 6 Jun 2025). Matrix autoregressive Yule–Walker estimation has been compared with VAR and Burg-type procedures, with similar MAE and RMSE but heavier computation in high-dimensional settings (Kołodziejski, 21 May 2025).

A persistent source of confusion is terminological. The phrase “iterated generalized Yule–Walker estimator” is explicit in the periodic covariance-extension algorithm and in the high-dimensional matrix-valued spatio-temporal paper (Picci et al., 2016, Dou et al., 14 Aug 2025). It is not the natural description of rowwise least-squares generalized Yule–Walker estimation, sparse penalized generalized Yule–Walker estimation, regularized operator-valued Yule–Walker estimation, or changepoint-robust difference-based Yule–Walker estimation (Dou et al., 2015, Reuvers et al., 2021, Niu et al., 6 Jun 2025, Gallagher et al., 2021). This suggests that the term denotes a specific algorithmic subclass of generalized Yule–Walker methods: one in which the moment equations are solved by repeated fixed-point, quasi-Newton, or alternating conditional updates, rather than by a single closed-form, penalized, or regularized fit.