The Yule–Walker Equations
You have measured a series, drawn its
correlogram, and
decided an AR(p)
model fits. One question remains: what are the coefficients
\phi_1, \dots, \phi_p? The Yule–Walker equations answer it
with a short, exact piece of algebra — they convert the autocorrelations you can estimate from data
straight into the AR parameters, by solving one linear system. No optimiser, no iteration: just plug the
ACF into a matrix equation and solve.
The trick is a single, repeatable move. Start from the AR(p) recursion,
multiply both sides by a lagged value x_{t-k}, and take expectations. Because
the shock \varepsilon_t is uncorrelated with the past, its term drops out, and
you are left with a relation among autocovariances.
Deriving the equations
Write the model x_t = \phi_1 x_{t-1} + \dots + \phi_p x_{t-p} + \varepsilon_t.
Multiply through by x_{t-k} (for k \ge 1) and take
expectations. Using \gamma_k = \operatorname{Cov}(x_t, x_{t-k}) and
\mathbb{E}[\varepsilon_t x_{t-k}] = 0:
\gamma_k = \phi_1 \gamma_{k-1} + \phi_2 \gamma_{k-2} + \dots + \phi_p \gamma_{k-p}, \qquad k = 1, 2, \dots, p.
Divide every term by the variance \gamma_0 to get the same relation in
autocorrelations:
\rho_k = \phi_1 \rho_{k-1} + \phi_2 \rho_{k-2} + \dots + \phi_p \rho_{k-p}, \qquad k = 1, \dots, p.
That is p equations in the p unknowns
\phi_1, \dots, \phi_p — the Yule–Walker equations. Taking one more equation at
k = 0 recovers the innovation variance,
\sigma^2 = \gamma_0\big(1 - \phi_1\rho_1 - \dots - \phi_p\rho_p\big).
The matrix form: a Toeplitz system
Stack the p equations and the structure is unmistakable. With
\boldsymbol{\phi} = (\phi_1,\dots,\phi_p)^{\!\top},
\mathbf{r} = (\rho_1,\dots,\rho_p)^{\!\top}, and
R the matrix of autocorrelations
R_{ij} = \rho_{|i-j|}:
R\,\boldsymbol{\phi} = \mathbf{r}, \qquad R = \begin{pmatrix} 1 & \rho_1 & \cdots & \rho_{p-1} \\ \rho_1 & 1 & \cdots & \rho_{p-2} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{p-1} & \rho_{p-2} & \cdots & 1 \end{pmatrix}.
The matrix R is Toeplitz — constant down every diagonal —
and symmetric positive-definite, so the system always has a unique solution
\boldsymbol{\phi} = R^{-1}\mathbf{r}. Plugging the sample
autocorrelations \hat\rho_k into R and
\mathbf{r} and solving gives the Yule–Walker estimator of the
AR coefficients — a
least-squares-flavoured
method of moments: match the model's theoretical ACF to the observed one.
Worked example: an AR(2)
For p = 2 the two equations are
\rho_1 = \phi_1 + \phi_2\,\rho_1, \qquad \rho_2 = \phi_1\,\rho_1 + \phi_2.
Solving this 2\times 2 system for the coefficients gives closed forms:
\phi_1 = \frac{\rho_1(1 - \rho_2)}{1 - \rho_1^{2}}, \qquad \phi_2 = \frac{\rho_2 - \rho_1^{2}}{1 - \rho_1^{2}}.
Suppose the sample ACF gives \rho_1 = 0.5 and
\rho_2 = 0.4. Substitute:
| Quantity | Expression | Value |
|---|
| 1 - \rho_1^2 | 1 - 0.25 | 0.75 |
| \phi_1 | 0.5(1-0.4)/0.75 | 0.4 |
| \phi_2 | (0.4-0.25)/0.75 | 0.2 |
So the fitted model is x_t = 0.4\,x_{t-1} + 0.2\,x_{t-2} + \varepsilon_t — read
directly off two autocorrelations, no iteration required.
The ACF that produced it
The two numbers we fed in, \rho_1 = 0.5 and
\rho_2 = 0.4, are the first two bars below; the rest of the ACF then unrolls by
the same recursion \rho_k = 0.4\,\rho_{k-1} + 0.2\,\rho_{k-2}, tailing off
smoothly. Yule–Walker is just this picture read backwards: given the first
p bars, recover the coefficients that generate the whole thing.
A faster route, and a link back to the PACF
Inverting R directly costs O(p^3), but the Toeplitz
structure allows the Durbin–Levinson recursion to solve the system in
O(p^2) by building the order-k solution from the
order-(k-1) one. Its by-products are exactly the
partial
autocorrelations \phi_{kk} — so Yule–Walker and the PACF are two
faces of the same recursion. Fit AR models of increasing order and the last coefficient at each stage
is the PACF, which is why the PACF cuts off at the true order p.
The derivation leaned on one crucial fact: \varepsilon_t is uncorrelated with
the past values x_{t-k}, so its term vanished. That is true for a pure AR — but
the moment an MA
term enters, the shocks are correlated with recent values, the clean cancellation fails, and the
equations become nonlinear in the parameters. Applying Yule–Walker to an
ARMA or MA model gives
badly biased, inefficient estimates. For those you switch to maximum-likelihood or
conditional-least-squares estimation, which handle the shock correlations properly. Yule–Walker's charm —
a one-shot linear solve — is precisely a privilege of the pure-AR world.
Udny Yule, the English statistician, wrote the founding
1927 paper that modelled the sunspot cycle as an autoregression — and derived these equations to fit it.
Sir Gilbert Walker, a British physicist working in India, independently hit the same
relations while hunting for predictors of the monsoon (the same Walker whose name graces the "Walker
circulation" and who first quantified the Southern Oscillation behind El Niño). Two scientists, two
continents, two utterly different problems — sunspots and monsoons — converging on one little linear
system. It is a lovely reminder that the mathematics of memory-plus-noise is genuinely universal.