ARIMA and Box–Jenkins

Exponential smoothing models components — level, trend, season. ARIMA models dependence itself: it asks what linear machine, fed white noise, would produce a series with exactly the correlations yours has. That inversion of viewpoint — from “what is the series made of” to “what process generates these correlations” — is the deepest idea in classical time series, and it comes with the field’s best diagnostic skill: reading ACF/PACF fingerprints. This lesson trains that skill on a simulator where you control the truth.

Two atoms: AR and MA

Autoregression — the present is a linear function of the past plus a shock:

Memory is recursive: a shock at time $t$ echoes forever, decaying geometrically (for AR(1), $\rho_k = \phi_1^k$ exactly). Stationarity demands the echo fade — for AR(1), $|\phi_1| < 1$; in general the roots of the AR polynomial must lie outside the unit circle, and at $\phi_1 = 1$ the process is the random walk, the unit root of lesson 1.

Moving average — the present is a blend of recent shocks:

Memory is finite: after $q$ steps a shock is gone completely, so the ACF cuts off dead at lag $q$ — the sharpest fingerprint in the field.

The PACF (partial autocorrelation — correlation of $y_t$ with $y_{t-k}$ after projecting out lags $1..k{-}1$, the Frisch–Waugh idea from the Projections lesson working a time-series job) mirrors the pattern. The identification table earning its keep for seventy years:

The backshift algebra, and the “I”

With the backshift operator $B y_t = y_{t-1}$, the whole family compresses:

where $(1-B)^d$ is lesson 1’s differencing, promoted to a model component: “integrated” of order $d$. The algebra is not just notation — it is how the equivalences come out: SES is ARIMA(0,1,1); Holt’s linear method is ARIMA(0,2,2); a seasonal factor $(1 - B^m)$ gives SARIMA, whose $(1,1,1)(1,1,1)_{12}$ workhorse fits most monthly business series. The two classical pillars are one algebra wearing two interfaces.

Box–Jenkins, then and now

The classical loop: identify ($d$ by unit-root tests and ACF decay; $p, q$ from the fingerprints) → estimate (MLE — the Probability module’s engine, with the Gaussian likelihood built by the state-space form) → check (residual ACF white? Ljung–Box test aggregate the lags?) → iterate. What auto.arima automated is the identify step, by AIC search over $(p, d, q)$ — and honesty requires saying it usually wins: the human fingerprint reader’s residual value is (a) sanity-checking the search on weird series, (b) reading residual diagnostics, which no AIC search does for you, and (c) knowing that near-cutoff cases are genuinely ambiguous (an AR(1) with φ = 0.5 at n = 100 is statistically hard to tell from ARMA(1,1) — the fingerprints blur, and forecast accuracy barely cares).

Boundaries, honestly drawn: ARIMA is linear with constant variance — volatility clustering needs GARCH, regime switches need state-space (next lesson), and nonlinear dependence needs the neural lesson. And like ETS, an ARIMA fitted to one series learns nothing from the thousand sibling series next to it — the limitation the global-model revolution attacks.

import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.stats.diagnostic import acorr_ljungbox

rng = np.random.default_rng(1)
y = ArmaProcess(ar=[1, -0.6], ma=[1]).generate_sample(400, distrvs=rng.standard_normal)

fit = ARIMA(y, order=(1, 0, 0)).fit()
print(fit.params[1])                                   # ≈ 0.6: MLE recovers φ
print(acorr_ljungbox(fit.resid, lags=[10]).lb_pvalue)  # large p: residuals white ✓

wrong = ARIMA(y, order=(0, 0, 1)).fit()                # deliberately misspecified
print(acorr_ljungbox(wrong.resid, lags=[10]).lb_pvalue)  # small p: it shows

Exercises