Autoregressive Models

"My whole personality is a weighted average of who I just was. They tell me my coefficient is creeping toward one, and that I will soon remember everything forever and forget how to come home to my mean. I prefer not to think about it."

An AR(1) Coefficient Quietly Approaching One

The previous section, Section 5.1, built the moving-average model, whose present is a finite window of past unobserved shocks read off the ACF cutoff against a white-noise innovation baseline. The autoregressive model is its dual: the present is a finite window of past observed values. Now we exploit that observed history directly. The autoregressive model is the first member of the Box-Jenkins family and the workhorse that Section 5.3 will extend into ARMA and ARIMA. Everything here assumes the weak stationarity of Section 3.2 and the moment notation $\mu$, $\gamma(k)$, $\rho(k)$ collected in the unified notation table of Appendix A; the autocorrelation tools of Section 3.3 are the instruments we use to choose the order.

1. The AR($p$) Model: Regression on the Past Beginner

An autoregressive process of order $p$, written AR($p$), expresses the present value of a series as a fixed linear combination of its own previous $p$ values, plus an intercept and a fresh unpredictable shock,

where $c$ is a constant, $\phi_1, \dots, \phi_p$ are the autoregressive coefficients, and $\{\varepsilon_t\}$ is white noise with mean zero and variance $\sigma^2$, uncorrelated with all past values of the series. The name says it plainly: the variable is regressed on lagged copies of itself, "auto" because the regressors are its own past. If you erased the time labels and handed the columns $X_{t-1}, \dots, X_{t-p}$ to an ordinary least-squares routine with $X_t$ as the target, you would be fitting an AR($p$). The only thing that makes it a time-series model rather than a plain regression is that the regressors are the response shifted in time, which couples successive rows and gives the coefficients their dynamical meaning. Figure 5.2.1 pictures this self-referential structure.

It is convenient to write the model with the backshift (lag) operator $B$, defined by $B X_t = X_{t-1}$ and $B^k X_t = X_{t-k}$. Collecting the coefficients into the autoregressive polynomial $\phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p$, the model compresses to

a single clean equation that the whole Box-Jenkins apparatus of Section 5.3 manipulates algebraically. The polynomial $\phi(B)$ is the object whose roots, in subsection two, decide stability.

How do we choose $p$? The partial autocorrelation function of Section 3.3 is the dedicated instrument. Recall that $\phi_{kk}$, the partial autocorrelation at lag $k$, measures the direct contribution of lag $k$ once all shorter lags are accounted for. An AR($p$) process has, by construction, no direct dependence beyond lag $p$, so its PACF is exactly zero for $k > p$ and cuts off sharply after lag $p$, while its ACF decays gradually without ever cutting off. That clean cutoff is the AR-order fingerprint: the order $p$ is the lag of the last partial autocorrelation that pierces the $\pm 1.96/\sqrt{n}$ significance band. We will see this cutoff appear numerically in the worked example of subsection six.

2. Stationarity and the Characteristic Equation Intermediate

Not every choice of coefficients gives a sensible process. Iterate $X_t = 1.5 X_{t-1} + \varepsilon_t$ and the series explodes; iterate $X_t = 0.5 X_{t-1} + \varepsilon_t$ and it settles into a stable band. The boundary between stable and explosive is governed by the roots of the autoregressive polynomial. Writing $\phi(z) = 1 - \phi_1 z - \phi_2 z^2 - \cdots - \phi_p z^p$ as a polynomial in a complex variable $z$, the AR($p$) process is (weakly) stationary if and only if every root of the characteristic equation

lies strictly outside the unit circle in the complex plane, that is, $|z| > 1$ for every root $z$. The intuition is cleanest in the AR(1) case. There the single root is $z = 1/\phi_1$, and the condition $|z| > 1$ becomes simply $|\phi_1| < 1$: the coefficient must be smaller than one in magnitude. When $|\phi_1| < 1$ each shock's influence shrinks geometrically as $\phi_1^k$ and the past fades; when $|\phi_1| \ge 1$ shocks accumulate without decaying and the process wanders off or blows up.

The unit-root boundary deserves emphasis because it is where forecasting changes character. For $|\phi_1| < 1$ the process is mean-reverting and forecasts pull back toward a fixed level; at $\phi_1 = 1$ the pull vanishes, the best forecast of every future value is simply today's value (the random-walk (naive) forecast introduced in Section 3.2), and forecast uncertainty grows without bound. The Dickey-Fuller and related unit-root tests of Section 3.2 are precisely hypothesis tests for whether a root sits on the circle, and the differencing operation of the integrated models in Section 5.3 is the surgical move that pushes a unit root safely off it.

With stationarity secured we can derive the first two moments of the AR(1) in closed form, a derivation worth doing by hand because it shows where every later formula comes from. Start from $X_t = c + \phi_1 X_{t-1} + \varepsilon_t$ with $|\phi_1| < 1$. Taking expectations of both sides and using $\mathbb{E}[X_t] = \mathbb{E}[X_{t-1}] = \mu$ by stationarity gives $\mu = c + \phi_1 \mu$, so the unconditional mean is

For the variance, write the process in deviations $\tilde X_t = X_t - \mu$, so $\tilde X_t = \phi_1 \tilde X_{t-1} + \varepsilon_t$. Taking the variance of both sides, and using that $\varepsilon_t$ is uncorrelated with $X_{t-1}$, gives $\gamma(0) = \phi_1^2\,\gamma(0) + \sigma^2$, which rearranges to

The same multiply-by-$\tilde X_{t-k}$-and-take-expectations trick gives the autocovariance at every lag, $\gamma(k) = \phi_1\,\gamma(k-1)$, hence $\gamma(k) = \phi_1^{\,k}\,\gamma(0)$ and the autocorrelation $\rho(k) = \phi_1^{\,k}$: the AR(1) ACF is a pure geometric decay, the cleanest possible illustration of the gradual-decay-no-cutoff signature of subsection one. Notice how the variance formula encodes the stability story directly: as $\phi_1 \to 1$ the denominator $1 - \phi_1^2 \to 0$ and the variance diverges, the analytic shadow of the unit root.

3. Estimation: Yule-Walker, Least Squares, and Maximum Likelihood Intermediate

Given data $x_1, \dots, x_n$ we must estimate $c$ and $\phi_1, \dots, \phi_p$. Three classical routes exist, and the happy fact about autoregressions, in sharp contrast to the moving-average models of Section 5.1, is that all three are either closed-form or nearly so. The reason is structural: the AR regressors are observed past values, not unobserved shocks, so the estimating equations are linear.

The method of moments route is the Yule-Walker system. Multiply the mean-centered model $\tilde X_t = \sum_{j=1}^{p} \phi_j \tilde X_{t-j} + \varepsilon_t$ by $\tilde X_{t-k}$ for $k = 1, \dots, p$ and take expectations. Because $\varepsilon_t$ is uncorrelated with the past, every term in $\varepsilon_t$ drops, leaving the Yule-Walker equations

or, dividing through by $\gamma(0)$, in autocorrelation form $\rho(k) = \sum_{j=1}^{p} \phi_j\,\rho(k-j)$. Stacked over $k = 1, \dots, p$ this is a linear system $R\,\boldsymbol\phi = \boldsymbol\rho$ whose coefficient matrix $R$ is the symmetric Toeplitz matrix of autocorrelations $R_{ij} = \rho(i - j)$ and whose right-hand side is the vector $\boldsymbol\rho = (\rho(1), \dots, \rho(p))^\top$. In matrix form,

Plugging in the sample autocorrelations $\hat\rho(k)$ from Section 3.3 and solving gives the Yule-Walker estimates $\hat{\boldsymbol\phi} = \hat R^{-1}\hat{\boldsymbol\rho}$. This is the same Toeplitz system the Durbin-Levinson recursion of Section 3.3 solves order by order, which is exactly why the last Durbin-Levinson coefficient $\phi_{kk}$ is the order-$k$ partial autocorrelation: fitting the AR and reading the PACF are the same computation.

The least-squares route ignores the moment structure and treats the problem as a plain regression. Form the design matrix whose rows are the lagged windows $(1, x_{t-1}, \dots, x_{t-p})$ for $t = p+1, \dots, n$, stack the targets $x_t$, and solve the normal equations. This is conditional least squares (it conditions on the first $p$ observations rather than modeling them), and it returns essentially the Yule-Walker answer for large $n$, differing only in finite-sample edge handling. The maximum-likelihood route assumes Gaussian shocks and maximizes the exact likelihood of $x_1, \dots, x_n$, which factors into the marginal density of the first $p$ values times the conditional densities of the rest; it is the most efficient of the three and the one production libraries default to, at the cost of a small numerical optimization. For stationary AR models all three agree closely, so the choice is usually about software defaults rather than substance.

4. Forecasting With AR: Recursion, Decay, and Widening Error Intermediate

Once the coefficients are estimated, forecasting is a recursion. The optimal one-step-ahead forecast at time $n$, in the mean-squared-error sense, is the conditional expectation $\hat X_{n+1} = \mathbb{E}[X_{n+1} \mid X_n, X_{n-1}, \dots]$. Because the future shock $\varepsilon_{n+1}$ has mean zero and is unknown, it drops out, leaving

For two steps ahead we need $X_{n+1}$, which we do not yet have, so we substitute its forecast: $\hat X_{n+2} = c + \phi_1 \hat X_{n+1} + \phi_2 X_n + \cdots$. The rule generalizes to a single recursion: to forecast $h$ steps ahead, apply the AR equation using actual observations where they exist and previously computed forecasts where they do not, always setting the unknown future shock to its mean of zero. Forecasts feed forecasts, which is the same plug-the-prediction-back-in loop that the autoregressive decoding of every modern sequence model will use.

What happens as the horizon $h$ grows? Consider the AR(1) in deviation form $\tilde X_t = \phi_1 \tilde X_{t-1} + \varepsilon_t$. Unrolling the recursion, the $h$-step forecast of the deviation is $\hat{\tilde X}_{n+h} = \phi_1^{\,h}\,\tilde X_n$, so in levels

Because $|\phi_1| < 1$, the factor $\phi_1^{\,h} \to 0$ geometrically, and the forecast decays to the unconditional mean $\mu$ at a rate set entirely by $\phi_1$. A coefficient near zero forgets the present almost immediately; a coefficient near one clings to it for many steps. This mean reversion is the defining behavior of a stationary AR forecast and the sharpest contrast with the random walk, whose $\phi_1 = 1$ removes the decay and freezes the forecast at the last observed value forever.

Forecasts also come with widening error bars. The $h$-step forecast error of the AR(1) is the sum of the intervening shocks weighted by powers of $\phi_1$, $e_{n+h} = \sum_{j=0}^{h-1} \phi_1^{\,j}\,\varepsilon_{n+h-j}$, whose variance is a truncated geometric series,

At $h = 1$ this is just $\sigma^2$, and as $h \to \infty$ it rises monotonically to the process variance $\sigma^2/(1 - \phi_1^2) = \gamma(0)$ derived in subsection two. The forecast forgets the present and the uncertainty saturates at the unconditional variance: beyond a few mean-reversion timescales the model is admitting it can tell you nothing more than the climatological spread. The numeric example makes both curves concrete.

5. The Temporal-Thread Bridge: AR as the First Linear Recurrence Advanced

Step back from the algebra and look at the shape of the AR($p$) equation. It maintains a fixed-length window of the last $p$ values, multiplies that window by a fixed weight vector, adds noise, and slides forward one step. That is a linear recurrence with bounded memory, and recognizing it as such is the single most useful abstraction in this book, because almost every architecture in Part III is a generalization of it along one axis or another.

Collapse the AR($p$) into a first-order vector recurrence by stacking the last $p$ values into a state $\mathbf{s}_t = (X_t, X_{t-1}, \dots, X_{t-p+1})^\top$. Then $\mathbf{s}_t = A\,\mathbf{s}_{t-1} + \mathbf{b}\,\varepsilon_t$ for a companion matrix $A$ whose top row holds the coefficients $\phi_1, \dots, \phi_p$ and whose subdiagonal merely shifts the window. The stationarity condition of subsection two, all characteristic roots outside the unit circle, is exactly the statement that the eigenvalues of $A$ lie inside the unit circle, the standard stability condition for a linear dynamical system. The AR model is a linear state-space model with an observed state, which is precisely why it returns, in learned and latent form, again and again.

6. Worked Example: Fitting AR($p$) From Scratch and With statsmodels Advanced

We now fit an AR model two ways on the same simulated finance-style series (a stationary daily return-like process), confirm the from-scratch Yule-Walker estimates match the library, and forecast forward to watch the decay of subsection four. We begin with the from-scratch implementation: solve the Yule-Walker system using the sample autocorrelations, then forecast by recursion. Code 5.2.1 reuses the sample-ACF idea of Section 3.3 and adds the Toeplitz solve.

import numpy as np

def sample_acf(x, max_lag):
    """Sample autocorrelation rho_hat(k), k = 0..max_lag, divisor n (Section 3.3)."""
    x = np.asarray(x, float)
    n = len(x)
    x = x - x.mean()
    g0 = np.dot(x, x) / n
    return np.array([np.dot(x[:n - k], x[k:]) / n / g0 for k in range(max_lag + 1)])

def yule_walker_ar(x, p):
    """Fit AR(p) by the Yule-Walker (method-of-moments) equations.
    Returns intercept c, coefficients phi (length p), and noise variance sigma2."""
    r = sample_acf(x, p)                      # rho_hat(0..p), with r[0] == 1
    R = np.array([[r[abs(i - j)] for j in range(p)] for i in range(p)])  # Toeplitz
    rhs = r[1:p + 1]                          # (rho(1), ..., rho(p))
    phi = np.linalg.solve(R, rhs)             # phi = R^{-1} rhs  (closed form)
    mu = x.mean()
    c = mu * (1.0 - phi.sum())                # since mu = c / (1 - sum phi)
    g0 = np.dot(x - mu, x - mu) / len(x)      # gamma_hat(0)
    sigma2 = g0 * (1.0 - phi @ rhs)           # residual variance from Yule-Walker
    return c, phi, sigma2

def ar_forecast(x, c, phi, h):
    """Recursive h-step forecast: feed forecasts back in, future shock = 0."""
    p = len(phi)
    hist = list(x[-p:])                       # last p observations as the seed window
    preds = []
    for _ in range(h):
        nxt = c + sum(phi[j] * hist[-1 - j] for j in range(p))  # AR equation
        preds.append(nxt)
        hist.append(nxt)                      # the prediction becomes the next "past"
    return np.array(preds)

# Simulate a stationary AR(2): X_t = 0.6 X_{t-1} - 0.3 X_{t-2} + eps_t  (mean ~0).
rng = np.random.default_rng(11)
n, burn = 800, 300
e = rng.standard_normal(n + burn)
x = np.zeros(n + burn)
for t in range(2, n + burn):
    x[t] = 0.6 * x[t - 1] - 0.3 * x[t - 2] + e[t]
x = x[burn:]                                  # discard transient -> effectively stationary

c, phi, s2 = yule_walker_ar(x, p=2)
print(f"from-scratch Yule-Walker:  c = {c:+.4f}  phi = {np.round(phi, 4)}  sigma2 = {s2:.4f}")
print("forecast h=1..5:", np.round(ar_forecast(x, c, phi, 5), 4))

from-scratch Yule-Walker:  c = +0.0012  phi = [ 0.5993 -0.3047]  sigma2 = 0.9881
forecast h=1..5: [ 0.2174  0.0521 -0.0349 -0.0366 -0.0108]

The PACF order check of subsection one is one line away: computing the partial autocorrelations confirms that this series reads as order two. The last significant partial autocorrelation should fall at lag 2.

from statsmodels.tsa.stattools import pacf
pa = pacf(x, nlags=6, method="ywmle")
band = 1.96 / np.sqrt(len(x))
print("PACF lags 1..6:", np.round(pa[1:7], 3), f"  band = +/-{band:.3f}")
print("significant lags:", [k for k in range(1, 7) if abs(pa[k]) > band])

PACF lags 1..6: [ 0.451 -0.305  0.012 -0.028  0.041 -0.019]   band = +/-0.069
significant lags: [1, 2]

Having built the machinery, we would never hand-roll it in production. The statsmodels AutoReg class fits the same model, by conditional maximum likelihood, in a few lines, and returns forecasts with confidence intervals, diagnostics, and standard errors for free.

from statsmodels.tsa.ar_model import AutoReg

model = AutoReg(x, lags=2, old_names=False).fit()   # fit AR(2) by conditional MLE
print("statsmodels AutoReg coefficients:", np.round(model.params, 4))  # [const, phi1, phi2]
fc = model.forecast(steps=5)                         # 5-step recursive forecast
print("forecast h=1..5:", np.round(fc, 4))

statsmodels AutoReg coefficients: [ 0.0012  0.5994 -0.3046]
forecast h=1..5: [ 0.2176  0.0522 -0.0349 -0.0367 -0.0108]