Yule Walker Estimation: AR Model in Matlab and Python

Yule Walker estimation fits an auto-regressive model to a time series. Suppose the series is modeled as an AR(N) process. The model parameters $a_1,\dots,a_N$ can be estimated by least squares, but the equations become cumbersome as the order $N$ grows. The Yule-Walker equations are a cleaner route. They tie the AR coefficients to the autocorrelation of the process, so the coefficients can be solved for once the autocorrelation is estimated.

This post derives the equations, then simulates an AR(3) process in Matlab and in Python.

Learning path: AR models from random signals

Work through these in order, then come back here:

1. An introduction to estimation theory: what makes an estimator good.
2. Simulation and analysis of white noise in Matlab: the noise that drives an AR model.
3. AutoCorrelation (Correlogram) and persistence: the autocorrelation that Yule-Walker uses.
4. Shaping randomness: AR, MA, and ARMA models: what the model is.
5. Why can’t I just use a matrix to solve ARMA?: why the least-squares route is hard.
6. Yule-Walker estimation (this post): solving for the parameters.
7. Generate color noise using an auto-regressive model: using the fitted model to generate data.

Yule-Walker equations: the idea

The Yule-Walker equations relate the AR parameters to the autocorrelation $r_{xx}[k]$ of $x[n]$. To use them:

• Estimate the autocorrelation $\hat{r}_{xx}[k]$ from the data $x[n]$.
• Solve the Yule-Walker equations for $a_k$ $(k=1,\dots,N)$ and for $\sigma^2$.

The model order $N$ must be chosen first. The equations do not reveal it. These methods help choose it:

• ACF and PACF. The autocorrelation and partial autocorrelation functions show whether the data look like an AR process, and suggest an order.
• Akaike Information Criterion (AIC) [1]. Balances goodness of fit against the number of parameters. A penalty term stops the order from growing without limit, since a higher order always fits the data better.
• Bayesian Information Criterion (BIC) [2]. Similar to AIC, with a heavier penalty.
• Cross-validation [3]. Split the series into subsets, estimate on one subset, and check the fit on the others.

Deriving the Yule-Walker equations

An AR(N) model driven by white noise $w[n]$ of variance $\sigma^2$ is

$$x[n] + a_1 x[n-1] + a_2 x[n-2] + \cdots + a_N x[n-N] = w[n] \quad (1)$$

With $a_0 = 1$, this has the compact form

$$\sum_{k=0}^{N} a_k\, x[n-k] = w[n], \quad a_0 = 1 \quad (2)$$

Multiply (2) by $x[n-l]$ and take expectations:

$$\sum_{k=0}^{N} a_k\, E{x[n-k]\,x[n-l]} = E{w[n]\,x[n-l]} \quad (3)$$

For a wide-sense stationary process, $E{x[n-k]\,x[n-l]} = r_{xx}[l-k]$, where $r_{xx}[m] = E{x[n+m]\,x[n]}$. Write the cross term as $r_{wx}[l] = E{w[n]\,x[n-l]}$. Then (3) becomes

$$\sum_{k=0}^{N} a_k\, r_{xx}[l-k] = r_{wx}[l] \quad (4)$$

The next step is to evaluate $r_{wx}[l]$. Shift (1) by $l$ and solve for $x[n-l]$:

$$x[n-l] = -\sum_{k=1}^{N} a_k\, x[n-k-l] + w[n-l] \quad (5)$$

Substitute (5) into the cross term:

$$r_{wx}[l] = E{w[n]\,x[n-l]} = -\sum_{k=1}^{N} a_k\, E{w[n]\,x[n-k-l]} + E{w[n]\,w[n-l]}$$

For $l \ge 0$ the first term vanishes. Each $x[n-k-l]$ with $k \ge 1$ depends only on noise samples before time $n$, and white noise is uncorrelated with its past. The second term is the noise autocorrelation, which is $\sigma^2$ at lag 0 and zero at all other lags. Therefore

$$r_{wx}[l] = \begin{cases} \sigma^2, & l = 0 \ 0, & l \ge 1 \end{cases} \quad (6)$$

Substitute (6) into (4) and keep only $l \ge 1$:

$$\sum_{k=0}^{N} a_k\, r_{xx}[l-k] = 0, \quad l \ge 1$$

Since $a_0 = 1$, this rearranges to the Yule-Walker equations:

$$\sum_{k=1}^{N} a_k\, r_{xx}[l-k] = -r_{xx}[l], \quad l = 1,\dots,N \quad (7)$$

The sum now starts at $k=1$. For a real process the autocorrelation is symmetric, $r_{xx}[-m] = r_{xx}[m]$, so negative lags can be replaced by positive ones.

Writing (7) for $l = 1,\dots,N$ gives $N$ linear equations in $N$ unknowns:

$$\bar{R}\,\bar{a} = -\bar{r}, \qquad \bar{R}[l,k] = r_{xx}[\lvert l-k \rvert], \quad \bar{a}[k] = a_k, \quad \bar{r}[l] = r_{xx}[l], \quad l,k = 1,\dots,N \quad (8)$$

In compact form:

$$\bar{R}\,\bar{a} = -\bar{r} \quad (8)$$

$$\bar{a} = -\bar{R}^{-1}\bar{r} \quad (9)$$

$\bar{R}$ is a symmetric Toeplitz matrix. This structure allows fast solvers such as the Levinson-Durbin recursion [4][5]. Matlab’s aryule uses that approach.

Setting $l = 0$ in (4), with $r_{wx}[0] = \sigma^2$ and $r_{xx}[-k] = r_{xx}[k]$, gives the noise variance once the $a_k$ are known:

$$\sigma^2 = \sum_{k=0}^{N} a_k\, r_{xx}[k] \quad (10)$$

In practice $r_{xx}[k]$ is replaced by a sample estimate computed from the data.

Simulation in Matlab

Generate an AR(3) process and pretend the model parameters are unknown:

$$x[n] = -0.9\,x[n-1] – 0.6\,x[n-2] – 0.5\,x[n-3] + w[n]$$

Here $w[n]$ is unit-variance Gaussian white noise. The data is generated with filter, which realizes the all-pole model with coefficients $a = [1, 0.9, 0.6, 0.5]$.

rng(1);                                 % fixed seed for repeatable results
a = [1 0.9 0.6 0.5];                    % true model parameters, a0 = 1
x = filter(1, a, randn(1000,1));        % AR(3) data from 1000 Gaussian samples
plot(x);
title('Sample data from AR(3) process');
xlabel('Sample index');
ylabel('Value');

Plot the periodogram for reference. The PSDs of the estimated models will be compared against it.

figure(2);
periodogram(x); hold on;                % reference PSD of the data

Next, fit three model orders and compare each model’s PSD with the reference:

N = [2, 3, 4];
for i = 1:numel(N)
    [a_est, p] = aryule(x, N(i));       % Yule-Walker estimate; p = sigma^2
    [H, w] = freqz(sqrt(p), a_est);     % model response: sqrt(p) / A(z)
    plot(w/pi, 20*log10(2*abs(H)/(2*pi)));  % one-sided PSD in dB
end
xlabel('Normalized frequency (times pi rad/sample)');
ylabel('One-sided PSD (dB/rad/sample)');
legend('PSD estimate of x', 'PSD of model output N=2', ...
    'PSD of model output N=3', 'PSD of model output N=4');

The table gives the Matlab estimates with rng(1). Other seeds give slightly different numbers.

Model order N	$\hat{a}_1$	$\hat{a}_2$	$\hat{a}_3$	$\hat{a}_4$	$\hat{\sigma}^2$
2	0.8340	0.1602	–	–	1.2856
3	0.9105	0.5584	0.4775	–	0.9925
4	0.9157	0.5645	0.4874	0.0110	0.9924

Table 1: Matlab estimates of the model parameters and noise variance. True values: $a_1 = 0.9$, $a_2 = 0.6$, $a_3 = 0.5$, $\sigma^2 = 1$.

Figure 3 plots the noise variance estimate for orders 1 to 6.

The $N=2$ fit leaves too much unexplained: its $\hat{\sigma}^2$ is 1.29 against a true 1. Going from $N=3$ to $N=4$ barely changes $\hat{\sigma}^2$ (0.9925 to 0.9924), and $\hat{a}_4 = 0.011$ is close to zero. So $N=3$ is the sensible order here.

Solving the Yule-Walker equations directly

aryule hides the steps. The snippet below solves (8) and (10) by hand for the same data. It gives the same coefficients as aryule.

N = 3;
r = xcorr(x, N, 'biased');
r = r(N+1:end);                         % r(k+1) holds lag k, k = 0..N
R = toeplitz(r(1:N));                   % R(l,k) = r(|l-k|), l,k = 1..N
rhs = r(2:N+1);                         % r(1)..r(N)
ak = -R \ rhs;                          % a1..aN, from (9)
aFull = [1; ak];
sigma2 = sum(aFull .* r(1:N+1));        % sigma^2, from (10)

Simulation in Python

The same method in NumPy and SciPy. solve_toeplitz solves the Toeplitz system (8) directly. The function computes the biased autocorrelation, solves for the coefficients, and gets the noise variance from (10).

import numpy as np
from scipy.linalg import solve_toeplitz
from scipy.signal import lfilter


def yule_walker(x, order):
    """Estimate AR(order) coefficients [1, a1..aN] and noise variance sigma^2."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    # Biased sample autocorrelation r[k] for lags k = 0..N
    r = np.array([np.dot(x[: n - k], x[k:]) / n for k in range(order + 1)])
    # Solve R a = -r, where R is the Toeplitz matrix built from r[0..N-1]
    a_tail = solve_toeplitz(r[:order], -r[1 : order + 1])
    a = np.concatenate(([1.0], a_tail))
    # Noise variance from the lag-0 Yule-Walker equation: sigma^2 = sum_k a_k r[k]
    sigma2 = np.dot(a, r)
    return a, sigma2


rng = np.random.default_rng(1)
a_true = np.array([1.0, 0.9, 0.6, 0.5])
x = lfilter([1.0], a_true, rng.standard_normal(1000))

for N in range(1, 7):
    a_hat, sigma2 = yule_walker(x, N)
    print(f"N={N}: a = {np.round(a_hat[1:], 4)}, sigma2 = {sigma2:.4f}")

Running this gives:

Model order N	$\hat{a}_1$	$\hat{a}_2$	$\hat{a}_3$	$\hat{a}_4$	$\hat{\sigma}^2$
1	0.6224	–	–	–	1.3497
2	0.7820	0.2564	–	–	1.2610
3	0.9044	0.6297	0.4773	–	0.9737
4	0.9107	0.6380	0.4892	0.0132	0.9736

Table 2: Python estimates with default_rng(1). NumPy and Matlab generate different random sequences for the same seed, so these differ from Table 1. The pattern is the same: $\hat{\sigma}^2$ levels off at N=3.

Compare the AR(3) model with the periodogram of the data:

Summary

Yule-Walker estimation turns an AR model into a small linear system built from the autocorrelation. Choose the model order with a criterion such as AIC or BIC, and check the choice against the data. In both simulations above, $N=3$ recovers the true model.

References

[1] Akaike H. (1998) Information Theory and an Extension of the Maximum Likelihood Principle. In: Parzen E., Tanabe K., Kitagawa G. (eds) Selected Papers of Hirotugu Akaike. Springer Series in Statistics. Springer, New York, NY. https://doi.org/10.1007/978-1-4612-1694-0_15

[2] McQuarrie, A. D. R., and Tsai, C.-L., 1998. Regression and Time Series Model Selection. World Scientific. https://doi.org/10.1142/3573

[3] Arlot, Sylvain, and Alain Celisse. “A survey of cross-validation procedures for model selection.” Statistics Surveys 4 (2010): 40-79. https://doi.org/10.1214/09-SS054

[4] J. Durbin. The fitting of time series models. Revue de l’Institut International de Statistique, 28:233-243, 1960. https://doi.org/10.2307/1401322

[5] G. H. Golub and C. F. Van Loan, “Matrix Computations”, Johns Hopkins University Press, third edition, 1996. https://openlibrary.org/isbn/9780801854149

Further reading: Autoregressive model (Wikipedia) and Levinson recursion (Wikipedia).

---

Build the simulations yourself. The GaussianWaves eBooks include complete, runnable simulation code, beyond what fits in a single blog post.

• Digital Modulations using Python: modem and BER simulations built from scratch with NumPy, SciPy, and Matplotlib. Paperback on Amazon, or PDF with instant download.
• Wireless Communication Systems in Matlab (Second Edition): end-to-end wireless simulations, from modulation and coding to OFDM and fading channels, with Matlab scripts.