Time Series

Real series carry trend, seasonality, and changing variance; the first job is always to strip those until the residual is stationary, because stationarity — constant mean, and a covariance depending only on lag — is the only regime where averaging over time substitutes for the repeated samples we never get. Once stationary, every model is one object, the backshift polynomial $\phi(B)X_t=\theta(B)Z_t$, and the roots of those polynomials decide causality, invertibility, the autocorrelation shape, and the spectrum. Running underneath is a second thread, the time ↔ frequency duality: the autocovariance and the spectral density hold the same information in two bases, and you transform to whichever makes the next step easy. Prerequisites cashed in: stationarity and Markov structure (track 5 §1, §6), MLE/CLT/χ² tests and the delta method (track 3), Toeplitz systems and projection (track 2), EM and the forward recursion (track 7 §6, §8).

1. The frame and stationarity

Picture to hold. A time series is one realization of a stochastic process indexed by time. You see a single path, not repeated samples — so you cannot estimate “the distribution at time $t$” by averaging many draws. The only thing that rescues estimation is an assumption that the process looks the same after a time shift: then averaging over time substitutes for averaging over realizations. That assumption is stationarity, and it is the hinge of the whole subject.

The working decomposition for a raw series:

$$X_t = m_t + s_t + Z_t$$

• $m_t$ = trend (slow secular movement),
• $s_t$ = seasonality (periodic, known period $d$, $s_{t+d}=s_t$),
• $Z_t$ = the stationary residual we actually model.

White noise is the baseline residual: mean zero, finite variance, uncorrelated across time.

• White noise (WN): $E Z_t = 0$, $\mathrm{Var}(Z_t)=\sigma^2<\infty$, $\mathrm{Cov}(Z_s,Z_t)=0$ for $s\neq t$.
• iid noise: WN and independent (stronger — kills nonlinear dependence too).
• Gaussian WN: iid noise with $Z_t\sim N(0,\sigma^2)$.

White noise has no exploitable structure — you cannot predict it from its past. So the entire modeling pipeline is: transform → fit → check the residuals are white noise. If they are, you have extracted all the linear signal; if they aren’t, there is structure left to model.

Strict vs weak stationarity

Strictly stationary: the joint distribution of $\{X_{t_1},\dots,X_{t_k}\}$ equals that of the shifted block $\{X_{t_1+h},\dots,X_{t_k+h}\}$ for every $h$ and every collection. The entire probabilistic law is shift-invariant.

That is far more than we can check or need. We weaken it to the first two moments:

Weakly (second-order) stationary: finite variance, and

• (i) $E X_t = \mu$ constant (no trend),
• (ii) $\mathrm{Cov}(X_t, X_{t+h}) = \gamma(h)$ depends on the lag $h$ only, not on $t$.

“Stationary” unqualified means weakly stationary — it is the realistic assumption. Two immediate consequences:

• The autocovariance collapses to a one-argument function $\gamma(h)$, and it is symmetric: $\gamma(-h)=\gamma(h)$ (since $\mathrm{Cov}(X_t,X_{t+h})=\mathrm{Cov}(X_{t+h},X_t)$).
• The autocorrelation is $\rho(h)=\gamma(h)/\gamma(0)$, with $\rho(0)=1$.

Strict ⟹ weak always (matching whole laws forces matching moments), but not conversely.

Counterexample (weak ⇏ strict). Build a series that is Gaussian WN on even $t$ and centered-and-rescaled $\chi^2$ WN on odd $t$, with the two arranged to share mean $0$ and variance $\sigma^2$. Every covariance is $0$ and the mean is constant, so it is weakly stationary — yet the marginal distribution differs between even and odd $t$, so it is not strictly stationary. Matching moments is not matching laws.

Special case (Gaussian). For a Gaussian process (every finite block is multivariate normal) the law is determined by mean and covariance, so weak ⟺ strict. This is why Gaussian assumptions are so convenient here.

Not stationary: a random walk (covariance grows with $t$), anything with a trend (mean depends on $t$), ARCH-type changing variance.

Diagnostic picture. Take two windows from the series and overlay their histograms. Stationary → the histograms overlap. Non-stationary (trend or changing variance) → they separate or spread.

Measuring dependence: ACF and the correlogram

Autocorrelation function $\rho(h)$ — the correlation of the series with itself $h$ steps back. We estimate it with the sample ACF:

$$r_k = \frac{\sum_{t=1}^{n-k}(X_t-\bar X)(X_{t+k}-\bar X)}{\sum_{t=1}^{n}(X_t-\bar X)^2}.$$

White-noise test. Under WN, for fixed lag $k$ and large $n$ (a CLT result — cross-link [[03_Statistical_Inference]]):

$$\sqrt{n}\,(r_1,\dots,r_k)^\top \xrightarrow{d} N(0, I), \qquad\text{i.e.}\qquad r_k \overset{\text{approx}}{\sim} N(0, 1/n).$$

So a single $r_k$ outside $\pm 1.96/\sqrt n$ is significant evidence against white noise at 5%. Plot the $r_k$ with bands at $\pm1.96/\sqrt n$ — that plot is the correlogram.

Numeric anchor. $n=100 \Rightarrow$ bands at $\pm 1.96/\sqrt{100}=\pm0.196\approx\pm0.2$. Lag 0 is always $1$.

Multiple-testing trap. With $k$ lags plotted, the expected number exceeding the bands under pure white noise is $E\Big[\sum_{i=1}^{k}\mathbf 1(|r_i|>1.96/\sqrt n)\Big]=\sum_{i=1}^k 0.05 = 0.05k.$ For $k=40$ lags you expect $2$ exceedances even if the data is genuinely white noise. Never read a single border-crossing as proof of structure; this is exactly why §5 needs a joint (portmanteau) test, not eyeballing.

Implications of §1.

• Stationarity is what makes “average over time = average over the ensemble” legitimate; without it, estimation has no ground to stand on.
• The modeling loop is transform until residuals pass the white-noise test. Everything downstream is machinery for that loop.
• $\gamma$ being a function of lag alone is what later lets prediction reduce to a Toeplitz linear system (§5) — cross-link [[02_Linear_Algebra]].

Core competency set §1.

• State weak stationarity (two conditions) and why $\gamma(-h)=\gamma(h)$.
• Reproduce the WN band $\pm1.96/\sqrt n$ and the expected-exceedances calculation.
• Give the weak-but-not-strict counterexample and the Gaussian weak⟺strict fact.

---

2. Getting to stationary: trend, season, variance

The residual $Z_t$ is what we model, so first remove $m_t$, $s_t$, and any mean-dependent variance.

Trend estimation

Three routes, increasing in assumption-freedom:

1. Parametric: posit a form, e.g. quadratic $\hat m_t=\hat\alpha+\hat\beta t+\hat\gamma t^2$ by least squares. Accurate if the form is right; gives clean forecasts; risks misspecification.
2. Smoothing/filtering: local average $\hat m_t=\frac{1}{2q+1}\sum_{j=-q}^{q}X_{t+j}$. Window $q$ is a bias–variance dial — small $q$ low bias/high variance, large $q$ the reverse. No forecasts at the endpoints, no clean extrapolation.
3. Isotonic (monotone trend): minimize $\sum_t (X_t-a_t)^2$ s.t. $a_1\le\cdots\le a_n$ — a convex program. Rarely used; can’t forecast; degenerate on already-monotone data (returns the data, residual $0$, learns nothing).

Caution (ML models extrapolate flat). A random forest fit to a trended series predicts a constant outside the training range — it averages the last values and cannot continue a trend. It sees no time structure. This is the recurring weakness of memory-based methods here.

Differencing — remove the trend without estimating it

Define the difference operator $\nabla X_t = X_t-X_{t-1}$.

Linear trend, worked. $X_t=at+Z_t$. Then $\nabla X_t = (at+Z_t)-(a(t-1)+Z_{t-1}) = a + Z_t - Z_{t-1}.$ Constant mean $a$, no trend. One difference kills a linear trend.

Quadratic trend, worked. $X_t=t^2$ (noiseless). $\nabla X_t = t^2-(t-1)^2 = 2t-1$ — still trending. Difference again: $\nabla^2 X_t = \nabla X_t - \nabla X_{t-1} = (2t-1)-(2(t-1)-1) = 2.$ A degree-$p$ polynomial trend needs $p$ differences.

Counterexample / caution (over-differencing). Differencing white noise introduces autocorrelation: if $X_t=Z_t$ then $\nabla X_t=Z_t-Z_{t-1}$ has $\mathrm{Cov}(\nabla X_t,\nabla X_{t-1})=-\sigma^2\neq 0$. So you can manufacture spurious structure by differencing too many times. Difference, check the ACF, stop at the lowest order that looks white.

Forecasting back through differencing. If $\nabla X_t$ has mean $c$, invert: $\hat X_t = c + X_{t-1}$. For $\nabla^2$: since $\nabla^2 X_t = X_t-2X_{t-1}+X_{t-2}$, forecast $\hat X_t = c + 2X_{t-1}-X_{t-2}$.

Stochastic trend. A random-walk-with-drift trend $m_t=m_{t-1}+\delta+W_t$ is not deterministic, but differencing still removes it — another reason to prefer differencing when you don’t want to commit to a trend form.

Seasonality

Model $X_t = s_t + Z_t$ with $s_t$ periodic of known period $d$ (monthly→$d{=}12$, weekly→$d{=}7$). Two approaches:

• Harmonic regression: $s_t = a_0 + \sum_{f=1}^{k}\big(a_f\cos(2\pi f t/d)+b_f\sin(2\pi f t/d)\big)$. Start at low frequency $f$ (the big oscillations), add higher $f$ for finer structure — a bias–variance dial again. Linear in $(a_f,b_f)$, so ordinary least squares fits it.
• Seasonal differencing: $\nabla_d X_t = X_t - X_{t-d}$. With $X_t=s_t+Z_t$ and $s_{t}=s_{t-d}$,

$$\nabla_d X_t = (s_t - s_{t-d}) + (Z_t - Z_{t-d}) = Z_t - Z_{t-d},$$

mean zero. You lose a full period $d$ of observations. For trend and season, seasonal-difference first, then ordinary-difference the result.

Variance-stabilizing transforms (delta method)

Heteroscedastic data — e.g. variance growing with the level — breaks stationarity through condition (i)/(ii) on the transformed scale. Suppose $\mathrm{Var}(X_t)=g(\mu_t)$ for known $g$. Find $f$ so that $Y_t=f(X_t)$ has roughly constant variance.

Derivation (delta method — cross-link [[03_Statistical_Inference]]). Taylor-expand $f$ about $\mu_t$:

$$f(X_t) \approx f(\mu_t) + f'(\mu_t)(X_t-\mu_t).$$

Here $f(\mu_t)$ is deterministic; only the last term is random. So

$$\mathrm{Var}(Y_t) \approx \big(f'(\mu_t)\big)^2\,\mathrm{Var}(X_t) = \big(f'(\mu_t)\big)^2 g(\mu_t).$$

Choose $f$ to make the right side constant.

Numeric anchors.

• $\mathrm{Var}(X_t)=C\mu_t$ (count-like). Try $f(x)=\sqrt x$: $\big(\tfrac{1}{2\sqrt{\mu}}\big)^2 C\mu = C/4$, constant. → square-root transform, standard for counts.
• $\mathrm{Var}(X_t)=C\mu_t^2$ (multiplicative / compounding). Try $f(x)=\log x$: $\big(\tfrac1\mu\big)^2 C\mu^2 = C$, constant. → log transform, standard for stock-like data.

Transforms are not unique ($f(x)=2\sqrt x-1$ works as well as $\sqrt x$); pick the simplest.

Implications of §2.

• Pipeline order: variance-stabilize → de-trend (difference or regress) → de-season → check residual ACF.
• Differencing trades estimation for robustness — you remove trend without committing to its shape — but over-differencing injects negative autocorrelation, so it has its own cost.
• The square-root/log dichotomy (variance ∝ mean vs ∝ mean²) is worth memorizing; it recurs in GLMs (Poisson vs lognormal).

Core competency set §2.

• Difference a linear and a quadratic trend by hand; invert a one-step and two-step difference for forecasting.
• Derive the variance-stabilizing condition $(f')^2 g = \text{const}$ and produce $\sqrt{\cdot}$ / $\log$ for the two canonical $g$‘s.
• Explain why over-differencing is harmful and how seasonal differencing works.

---

3. MA and AR models — the backshift polynomial

Now the residual is stationary. The two atoms are moving-average (build the present from past noise) and autoregressive (build the present from past values). The unifying device is the backshift operator $B$: $BX_t=X_{t-1}$, $B^j X_t = X_{t-j}$.

Moving average MA(q)

$$X_t = Z_t + \theta_1 Z_{t-1} + \cdots + \theta_q Z_{t-q} = \theta(B)Z_t,\qquad \theta(B)=1+\theta_1 B+\cdots+\theta_q B^q.$$

A finite linear combination of white-noise shocks (not a true average — just a linear model in the $Z$‘s).

MA(1) is weakly stationary — worked. $X_t=Z_t+\theta Z_{t-1}$, $\mathrm{Var}(Z_t)=\sigma^2$.

$$E X_t = 0,\qquad \mathrm{Var}(X_t) = \sigma^2 + \theta^2\sigma^2 = \sigma^2(1+\theta^2).$$ $$\mathrm{Cov}(X_t,X_{t-1}) = \mathrm{Cov}(Z_t+\theta Z_{t-1},\, Z_{t-1}+\theta Z_{t-2}) = \theta\sigma^2,$$ $$\mathrm{Cov}(X_t,X_{t-j}) = 0 \quad (j>1).$$

No $t$ appears → weakly stationary. Every MA(q) is weakly stationary (finitely many shocks, constant-in-$t$ covariances). The general autocovariance:

$$\gamma(h) = \sigma^2\sum_{j=0}^{q-h}\theta_j\theta_{j+h}\quad(0\le h\le q),\qquad \gamma(h)=0\ (h>q).$$

Signature in the ACF. MA(q) autocorrelations are nonzero up to lag $q$ then cut off sharply to zero. (Contrast AR, which decays exponentially without a hard cutoff — §3 AR.) Seeing a clean cutoff at lag $q$ is the fingerprint of MA(q).

Invertibility — pinning down the parameters

Two different MA models can generate the same process, so the parameters aren’t identified without a restriction.

The two-representations problem, worked (Gaussian). With $\sigma^2=1$, $X_t=Z_t+\theta Z_{t-1}$ has $\mathrm{Var}=1+\theta^2$, $\mathrm{Cov}_1=\theta$. Now set $\tilde Z_t=\theta Z_t$ (so $\mathrm{Var}(\tilde Z)=\theta^2$) and consider $\tilde X_t=\tilde Z_t+\tfrac1\theta\tilde Z_{t-1}$: $\mathrm{Var}(\tilde X_t)=\theta^2+1,\qquad \mathrm{Cov}(\tilde X_t,\tilde X_{t-1})=\tfrac1\theta\cdot\theta^2=\theta.$ Identical mean, variance, and covariances. For Gaussian noise that means identical law — two MA(1) models ($\theta$ and $1/\theta$) for one process.

The fix: require the shock at time $t$ to dominate, $|\theta|<1$, equivalently the polynomial has its root outside the unit circle:

$$|\theta|<1 \iff \theta(z)=1+\theta z\neq 0\ \ \forall\,|z|\le 1.$$

Definition. MA(q) is invertible if $\theta(z)\neq 0$ for all $|z|\le 1$ (all roots strictly outside the unit circle).

Theorem (why invertibility matters). A model is invertible iff the noise can be recovered as a convergent linear combination of the observations:

$$Z_t = \pi(B)X_t = \sum_{j\ge0}\pi_j X_{t-j},\qquad \sum_j|\pi_j|<\infty,\ \pi_0=1,$$

i.e. there is a stable $\mathrm{AR}(\infty)$ form $X_t=-\sum_{j\ge1}\pi_j X_{t-j}+Z_t$. This is exactly the form prediction needs.

Why we insist on it — worked. Non-invertible $X_t=Z_t+2Z_{t-1}$. Solving for the AR form by repeated substitution gives weights $1,2,4,8,16,\dots$ — they explode, so the infinite past has unbounded influence: useless for prediction. The invertible twin $\tilde X_t=\tilde Z_t+\tfrac12\tilde Z_{t-1}$ gives decaying weights $1,\tfrac12,\tfrac14,\dots$. Same process, but only the invertible representation is usable.

Numeric anchors (root test).

• $X_t=Z_t-Z_{t-1}+\tfrac14 Z_{t-2}$: $\theta(z)=1-z+\tfrac14 z^2=(1-\tfrac12 z)^2$, root $z=2>1$ → invertible.
• $X_t=Z_t-3Z_{t-1}+2Z_{t-2}$: $\theta(z)=1-3z+2z^2$, roots $\{1,\tfrac12\}$ → root on the boundary ⇒ not invertible.

Autoregressive AR(p)

$$X_t = \phi_1 X_{t-1}+\cdots+\phi_p X_{t-p} + Z_t,\qquad \phi(B)X_t = Z_t,\ \ \phi(B)=1-\phi_1 B-\cdots-\phi_p B^p.$$

Note the sign convention: the $\phi$‘s move to the same side as $X$ before forming the polynomial (opposite of MA, where $\theta(B)$ acts on the $Z$‘s).

The geometric-series inversion. Treat $\phi(B)$ as a number and divide:

$$\phi(B)X_t=Z_t \ \Rightarrow\ X_t = \frac{1}{\phi(B)}Z_t.$$

For AR(1), $\phi(B)=1-\phi B$, so if $|\phi|<1$ use $\frac{1}{1-x}=\sum_{j\ge0}x^j$:

$$X_t = (1+\phi B+\phi^2 B^2+\cdots)Z_t = \sum_{j\ge0}\phi^j Z_{t-j},$$

the $\mathrm{MA}(\infty)$ representation. (Rigorously: $B$ is a norm-1 shift operator, so the Neumann series converges iff $|\phi|<1$ — same unit-root condition as MA, in disguise.)

Lemma — the guess solves the recursion. With $X_t=\sum_{j\ge0}\phi^j Z_{t-j}$,

$$X_t = Z_t + \phi\sum_{j\ge0}\phi^j Z_{t-1-j} = Z_t + \phi X_{t-1}. \checkmark$$

Uniqueness — the remainder vanishes. Suppose $Y_t$ also solves $Y_t=\phi Y_{t-1}+Z_t$. Iterating:

$$Y_t = Z_t + \phi Z_{t-1} + \cdots + \phi^{k-1}Z_{t-k+1} + \phi^k Y_{t-k}.$$

The remainder’s size (stationary, so $E Y^2$ constant):

$$E[(\phi^k Y_{t-k})^2] = \phi^{2k}E[Y_0^2] \xrightarrow{k\to\infty} 0 \quad\text{iff } |\phi|<1.$$

So $Y_t=\sum_{j\ge0}\phi^j Z_{t-j}=X_t$: the unique stationary solution. (Non-stationary solutions exist — add $\alpha\phi^t$ for any $\alpha$ — but only one is stationary.)

AR(1) autocovariance via the MA(∞) form. Using $\gamma(h)=\sigma^2\sum_{j\ge0}\psi_j\psi_{j+h}$ with $\psi_j=\phi^j$:

$$\gamma(h)=\sigma^2\sum_{j\ge0}\phi^j\phi^{j+h}=\sigma^2\phi^h\sum_{j\ge0}\phi^{2j}=\frac{\sigma^2\phi^{h}}{1-\phi^2},\qquad \rho(h)=\phi^{h}.$$

Signature in the ACF. AR(1) autocorrelations decay exponentially $\rho(h)=\phi^h$ — no hard cutoff. $\phi>0$: smooth positive decay. $\phi<0$: decaying oscillation (alternating sign), a good visual cue that an AR term is needed.

Causality — depend only on the past

Mirror of invertibility, but for AR: we want $X_t$ expressible as a stable combination of past noise (an $\mathrm{MA}(\infty)$), not future noise.

Definition. AR(p) $\phi(B)X_t=Z_t$ is causal if $\phi(z)\neq0$ for all $|z|\le 1$ (roots outside the unit circle) — equivalently it has an $\mathrm{MA}(\infty)$ form $X_t=\sum_{j\ge0}\psi_j Z_{t-j}$.

Non-causal, worked. $X_t=2X_{t-1}+Z_t$, $\phi(z)=1-2z$, root $z=\tfrac12<1$ → not causal. Forward substitution gives weights $1,2,4,8,\dots$ (explode). Reversing time, $X_{t-1}=\tfrac12 X_t-\tfrac12 Z_t$, yields a convergent expansion — but in terms of future noise $Z_{t},Z_{t+1},\dots$. There is a stationary solution, but it depends on the future: useless.

Numeric anchors (root test).

• $X_t=0.9X_{t-1}+Z_t$: $\phi(z)=1-0.9z$, root $z=10/9>1$ → causal.
• Random walk $X_t=X_{t-1}+Z_t$: $\phi(z)=1-z$, root $z=1$ on the boundary → not causal (and not stationary). Fix by differencing: $\nabla X_t=Z_t$. When a root sits near 1, difference.
• $X_t=-X_{t-2}+Z_t$: $\phi(z)=1+z^2$, roots $\pm i$ on the unit circle → not causal.

The duality, stated cleanly:

process	the polynomial restriction	what it buys
MA(q)	invertible: $\theta(z)\neq 0$, $	z
AR(p)	causal: $\phi(z)\neq 0$, $	z

MA↔AR are two views of one model. $\mathrm{AR}(1)$ is an $\mathrm{MA}(\infty)$ (shown above); conversely an $\mathrm{MA}(\infty)$ with $\theta_j=\phi^j$ collapses to $\mathrm{AR}(1)$ by the same geometric-series trick run backwards. The notation you choose is whichever makes the next step easy: MA for the ACF, AR for prediction.

Implications of §3.

• Roots outside the unit circle is one condition wearing two names (causality for AR, invertibility for MA); both ensure the relevant infinite expansion converges.
• The ACF/PACF shapes (cutoff vs decay) are how you read the order off a plot — the practical payoff of all this algebra (§4 PACF).
• The geometric-series inversion is the workhorse you reuse for ARMA (§4) and for the spectral density (§7).

Core competency set §3.

• Show MA(1) is stationary and compute its $\gamma(h)$; state the MA(q) cutoff.
• Invert AR(1) to $\mathrm{MA}(\infty)$ via geometric series; prove existence (it solves the recursion) and uniqueness (remainder $\to0$).
• Derive $\rho(h)=\phi^h$ and $\gamma(h)=\sigma^2\phi^h/(1-\phi^2)$.
• Test causality/invertibility by the unit-root rule on three worked polynomials.

---

4. ARMA — combining, computing, and Yule–Walker

ARMA(p,q): $\phi(B)X_t=\theta(B)Z_t$. AR and MA together; the earlier models are special cases ARMA(0,0)=WN, ARMA(p,0)=AR, ARMA(0,q)=MA. Always assume $\phi$ and $\theta$ share no common roots (else you’re describing a simple process with needless parameters).

Common-factor counterexample. $X_t=X_{t-1}+Z_t-Z_{t-1}$ has $\phi(z)=1-z=\theta(z)$. It looks like ARMA(1,1), but substituting shows $X_t=Z_t$ — plain white noise, ARMA(0,0). Cancel common factors first. (Another: $\theta=(1-\tfrac13 z)(1-\tfrac12 z)$, $\phi=(1-\tfrac12 z)(1-\tfrac14 z)$ share $(1-\tfrac12 z)$ ⇒ reduce ARMA(2,2) to ARMA(1,1).)

Conditions. ARMA is causal iff $\phi(z)\neq0$ on $|z|\le1$ (⇒ $\mathrm{MA}(\infty)$), invertible iff $\theta(z)\neq0$ on $|z|\le1$ (⇒ $\mathrm{AR}(\infty)$).

Getting the MA(∞) form — polynomial division by geometric series. Want $\psi(z)=\theta(z)/\phi(z)$.

Worked, ARMA(1,1). $X_t=\tfrac13 X_{t-1}+Z_t-\tfrac12 Z_{t-1}$, so $\phi(z)=1-\tfrac13 z$, $\theta(z)=1-\tfrac12 z$. $\frac{\theta(z)}{\phi(z)}=\Big(1-\tfrac12 z\Big)\Big(1+\tfrac13 z+(\tfrac13 z)^2+\cdots\Big)$ $= 1 + \sum_{j\ge0} z^{j+1}\Big(\big(\tfrac13\big)^{j+1}-\tfrac12\big(\tfrac13\big)^{j}\Big),$ giving $X_t=Z_t+\sum_{j\ge0}\big((\tfrac13)^{j+1}-\tfrac12(\tfrac13)^j\big)Z_{t-j-1}$. (R: ARMAtoMA.)

Why bother with each representation? The MA(∞) gives the ACF; the AR(∞) gives predictions.

ACF of ARMA — the difference-equation route

Computing $\gamma$ from the $\mathrm{MA}(\infty)$ needs the whole infinite $\{\psi_j\}$ even for one lag — wasteful. Instead, take covariances of both sides of the ARMA equation against $X_{t-k}$.

Derivation. From $\phi(B)X_t=\theta(B)Z_t$, take $\mathrm{Cov}(\cdot, X_{t-k})$. LHS:

$$\mathrm{Cov}(X_t-\phi_1 X_{t-1}-\cdots-\phi_p X_{t-p},\,X_{t-k}) = \gamma(k)-\phi_1\gamma(k-1)-\cdots-\phi_p\gamma(k-p).$$

RHS (using that $Z_{t-j}$ correlates with $X_{t-k}$ only through the $\mathrm{MA}(\infty)$ weights $\psi$):

$$\mathrm{Cov}(\theta(B)Z_t,\,X_{t-k}) = \begin{cases}\sigma^2\sum_{j=k}^{q}\theta_j\psi_{j-k}, & 0\le k\le q,\\[2pt] 0, & k>q.\end{cases}$$

Equating gives a finite recursion for $\gamma$.

Yule–Walker (the AR(p) special case)

For pure AR(p) the RHS is $\sigma^2$ at $k=0$ and $0$ for $k>0$. Writing the recursion for $k=1,\dots,p$ and using $\gamma(-h)=\gamma(h)$:

$$\underbrace{\begin{pmatrix}\gamma(0)&\gamma(1)&\cdots&\gamma(p-1)\\ \gamma(1)&\gamma(0)&\cdots&\gamma(p-2)\\ \vdots& &\ddots&\vdots\\ \gamma(p-1)&\gamma(p-2)&\cdots&\gamma(0)\end{pmatrix}}_{\Gamma\ (\text{Toeplitz, symmetric})}\begin{pmatrix}\phi_1\\\phi_2\\\vdots\\\phi_p\end{pmatrix} = \begin{pmatrix}\gamma(1)\\\gamma(2)\\\vdots\\\gamma(p)\end{pmatrix},$$

i.e. $\Gamma\Phi=\gamma$. Then recover $\sigma^2$ from $k=0$: $\sigma^2=\gamma(0)-\phi_1\gamma(1)-\cdots-\phi_p\gamma(p)$.

$\Gamma$ is symmetric Toeplitz (constant along diagonals — cross-link [[02_Linear_Algebra]]) and full rank, hence invertible: the system always solves, and the Toeplitz structure makes it fast (Levinson–Durbin).

Worked AR(2) numeric. Sample autocorrelations $\hat\rho(1)=0.6$, $\hat\rho(2)=0.2$. Dividing the YW system by $\gamma(0)$ gives it in correlation form: $\begin{pmatrix}1&\rho(1)\\\rho(1)&1\end{pmatrix}\begin{pmatrix}\phi_1\\\phi_2\end{pmatrix}=\begin{pmatrix}\rho(1)\\\rho(2)\end{pmatrix}\ \Rightarrow\ \begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}\begin{pmatrix}\phi_1\\\phi_2\end{pmatrix}=\begin{pmatrix}0.6\\0.2\end{pmatrix}.$ Determinant $1-0.36=0.64$. By Cramer’s rule: $\phi_1=\frac{0.6-0.6\cdot0.2}{0.64}=\frac{0.48}{0.64}=0.75,\qquad \phi_2=\frac{0.2-0.6\cdot0.6}{0.64}=\frac{-0.16}{0.64}=-0.25.$ Then $\sigma_Z^2=\gamma(0)\big(1-\phi_1\rho(1)-\phi_2\rho(2)\big)=\gamma(0)(1-0.45+0.05)=0.60\,\gamma(0)$ — the AR(2) explains 40% of the variance. (Check causality: $\phi(z)=1-0.75z+0.25z^2$ has complex roots of modulus $\sqrt{1/0.25}=2>1$ ✓.)

When does Yule–Walker fail? When $\Gamma$ degenerates — all covariances equal (a noiseless, perfectly correlated series). Near-zero noise ⇒ near-singular $\Gamma$ ⇒ unstable estimates. Opposite of ordinary stats, where low noise helps.

PACF — reading AR order off a plot

Partial autocorrelation $\mathrm{pacf}(h)$ = the coefficient of $X_{t-h}$ in the best linear predictor of $X_t$ from $X_{t-1},\dots,X_{t-h}$ — the correlation of $X_t$ with $X_{t-h}$ after removing the intervening lags’ influence.

Why AR(p) PACF cuts off at $p$ — derivation. Regress $X_t$ on the previous $h$ lags; $\mathrm{pacf}(h)$ is the last coefficient $\phi_{hh}$.

• For $h=p$: the AR(p) equation $X_t=\phi_1 X_{t-1}+\cdots+\phi_p X_{t-p}+Z_t$ already is the best linear predictor from $p$ lags, with $Z_t\perp$ past. So the order-$p$ regression returns the true coefficients exactly, and $\phi_{pp}=\phi_p\neq0$.
• For $h>p$: write the best predictor from $h$ lags. Because $X_t=\sum_{j=1}^{p}\phi_j X_{t-j}+Z_t$ with $Z_t$ uncorrelated with all of $X_{t-1},\dots,X_{t-h}$ (causality ⇒ $Z_t$ independent of the past), the optimal weights put $\phi_1,\dots,\phi_p$ on the first $p$ lags and exactly zero on lags $p{+}1,\dots,h$ — adding older lags cannot reduce error once the $p$-lag predictor already achieves the irreducible $\mathrm{Var}(Z_t)$. Hence $\phi_{hh}=0$ for $h>p$. ∎

The mirror statement for MA(q) (ACF cuts off at $q$) is immediate from $\gamma(h)=0$ for $h>q$ (§3). The two cutoffs are duals:

	MA(q)	AR(p)
ACF	$0$ for $h>q$ (sharp cutoff)	decays to $0$ as $h\to\infty$
PACF	decays to $0$ as $h\to\infty$	$0$ for $h>p$ (sharp cutoff)

They are mirror images: ACF cuts off for MA; PACF cuts off for AR. A mixed ARMA shows neither clean cutoff — which is why you fit low-order AR or MA first and only escalate to ARMA if needed.

Numeric/diagnostic anchor. ACF cuts off after lag 2 → consider MA(2); but a single sharp PACF spike at lag 1 → an AR(1) may be the simpler truth. Start simple, check residuals.

Implications of §4.

• Common-factor cancellation enforces parsimony — the simplest equivalent model is the one you want.
• Yule–Walker turns AR estimation into a Toeplitz solve, fast and always solvable (away from the noiseless degeneracy).
• ACF/PACF cutoffs are the bridge from the algebra of §3 to actual model-order choice.

Core competency set §4.

• Divide $\theta/\phi$ by geometric series to get $\mathrm{MA}(\infty)$.
• Derive Yule–Walker by covariancing both sides; write $\Gamma\Phi=\gamma$ and note Toeplitz/invertible; solve a 2×2 AR(2) system by hand.
• Derive why AR(p) PACF cuts off at $p$ (causality ⇒ older lags get zero weight); state the full ACF/PACF cutoff table and use it to pick orders.
• Give the common-factor counterexample.

---

5. Fitting, prediction, and selection

Best linear prediction (the projection view)

With only weak stationarity we trust linear predictors (higher moments may drift; cross-link [[02_Linear_Algebra]] projection). Predict $Y$ from predictors $W_1,\dots,W_n$ by $\sum a_i W_i$; the optimal weights are

$$a^\star = \Delta^{-1}\zeta,\qquad \Delta_{ij}=\mathrm{Cov}(W_i,W_j),\ \ \zeta_i=\mathrm{Cov}(Y,W_i),$$

exactly $(X^\top X)^{-1}X^\top y$ with covariances in place of cross-products. For forecasting, the predictors are past values and $Y$ is the future value, so $X$ appears on both sides.

Under stationarity, $\Delta_{ij}=\gamma(i-j)$ (Toeplitz) and $\zeta_i=\gamma(i)$ — so prediction reduces to a Toeplitz solve.

One-step, worked. Predict $X_{t+1}$ from $X_t$: $\Delta=\gamma(0)$, $\zeta=\gamma(1)$, so $\hat X_{t+1}=X_t\,\gamma(1)/\gamma(0)=\rho(1)X_t$ — the lag-1 autocorrelation times the last value.

For a causal AR(p) with $n>p$, the best linear predictor is the model equation: $\hat X_{n+1}=\phi_1 X_n+\cdots+\phi_p X_{n+1-p}$. Gaussian noise ⇒ prediction interval $\hat X_{n+1}\pm 1.96\,\sigma$.

ARMA prediction: switch to the $\mathrm{AR}(\infty)$ form $X_t=Z_t-\sum_{j\ge1}\pi_j X_{t-j}$, truncate, and forecast $\hat X_t=-\sum_{j=1}^{J}\pi_j X_{t-j}$. Multi-step: iterate (predict one step, feed it back) — optimal by iterated expectations.

Long-horizon punchline. As the horizon grows, the prediction converges to the overall mean and the interval widens to $1.96\times$ the stationary SD (you’re taking powers of a matrix with eigenvalues $<1$ in modulus). If log-returns follow an invertible ARMA, the long-run forecast of log-returns is just their mean — no free lunch.

Parameter estimation — three methods

1. Yule–Walker / method of moments (AR only): match sample ACVF at lags $0..p$, solve $\Gamma\Phi=\gamma$ (§4). Linear, fast, but AR-only and not efficient for ARMA.
2. Conditional least squares (general ARMA): rewrite the noise recursively, e.g. MA(1) $Z_t=X_t-\mu-\theta Z_{t-1}$; set $Z_0=0$, recurse $\tilde Z_t$, minimize $\sum \tilde Z_t^2$. “Conditional” because it conditions on the made-up initial $Z_0=0$ — harmless asymptotically since the past’s influence decays exponentially. (ARMA(1,1) needs $X$-history too, so start at $t{=}2$.) Nonlinear for ARMA.
3. Maximum likelihood (Gaussian): maximize the exact likelihood from the autocovariance matrix. Same idea as conditional LS but handles the start-up exactly. Standard practice: warm-start MLE from conditional-LS estimates (R arima method="CSS-ML"). The three usually agree closely.

Model selection

• AIC $=-2\log\hat L + 2k$ (for ARIMA $k=p+q+2$, the $+2$ for the estimated mean and $\sigma$). Estimates out-of-sample error, asymptotically like CV → use for prediction. Tends to pick a slightly-too-large model.
• BIC $=-2\log\hat L + k\log n$. Heavier penalty ⇒ smaller models; model-selection consistent (recovers the true order as $n\to\infty$) → use to identify the true process. Small risk of choosing too-small a model, which can hurt prediction.
• Time-series CV. Fit on $X_1,\dots,X_t$, forecast $X_{t+1}$, record $e^\star_{t+1}=X_{t+1}-\hat X_{t+1}$; sweep $t=m,\dots,n-1$; average the squared errors. Train on the past, validate on the future only — ordinary CV would illegally train on the future to predict the past. Seasonal data: leave out whole periods; $k$-step goal: forecast $\hat X_{t+k}$.

Testing white noise — the portmanteau test

Eyeballing the correlogram fails the multiple-testing problem (§1). Use a joint statistic.

For white noise, $r_1,\dots,r_k\overset{iid}{\sim}N(0,1/n)$, so

$$Q = n\sum_{i=1}^{k} r_i^2 \ \xrightarrow{d}\ \chi^2_k.$$

After fitting an ARMA with estimated parameters, the degrees of freedom drop and we use the Ljung–Box–Pierce statistic

$$\tilde Q = n(n+2)\sum_{i=1}^{k}\frac{\hat r_i^2}{n-i}\ \xrightarrow{d}\ \chi^2_{k-p-q}.$$

Reject “residuals are white / model adequate” if $\tilde Q$ exceeds the $(1-\alpha)$ quantile (one-sided). Subtract $p+q$ because parameters were estimated from the data (cross-link [[03_Statistical_Inference]] — same df-deduction as in goodness-of-fit). Choose max lag $k$ carefully: too small misses far seasonal spikes; too large breaks the $\chi^2$ approximation. Pair the test with the residual ACF and standardized-residual plots — the test can miss what the eye catches and vice versa.

Implications of §5.

• Prediction is a Toeplitz projection; AR makes it trivial because the model is the predictor.
• AIC-for-prediction vs BIC-for-truth is the one-line rule to memorize.
• “Validate only on the future” is the non-negotiable difference from ordinary CV.

Core competency set §5.

• Write the best-linear-predictor formula and specialize to one-step ($\rho(1)X_t$) and causal AR.
• Contrast Yule–Walker / conditional-LS / MLE and say when each applies.
• State AIC vs BIC roles and the time-series CV protocol.
• Write the $Q$ and Ljung–Box statistics with the right df and explain the $-p-q$.

---

6. Extensions: ARIMA and seasonal models

Bundle differencing and seasonality into the ARMA machinery.

• ARIMA(p,d,q): difference $d$ times, then fit ARMA: $\phi(B)(1-B)^d X_t=\theta(B)Z_t$ (plus mean). Removes polynomial trends / random walks inside the model. (R arima, middle index $=d$.)
• Seasonal ARMA (period $s$): $\Phi(B^s)X_t=\Theta(B^s)Z_t$ — AR/MA terms only at the seasonal lags $s,2s,\dots$. Why not plain ARMA(2s,2s)? That would carry $2s$ AR + $2s$ MA parameters (e.g. $24+24$ for $s=12$) — unstable. S-ARMA fits 4.
• Multiplicative seasonal ARMA: $\Phi(B^s)\phi(B)X_t=\Theta(B^s)\theta(B)Z_t$ — captures short-range and seasonal dependence together with few parameters. Use when the ACF/PACF show spikes at the seasonal lags and in the small lags around them; the width of the surrounding spikes signals the order of the short-range terms.
• SARIMA: the full object — multiplicative seasonal ARMA with both ordinary ($\nabla^d$) and seasonal ($\nabla_s^D$) differencing:

$$\Phi(B^s)\,\phi(B)\,\nabla_s^D\,\nabla^d\,Y_t = \delta + \Theta(B^s)\,\theta(B)\,Z_t.$$

Fit step by step (plot → de-trend → ACF/PACF → low-order AR/MA → escalate); R sarima (package astsa) gives the diagnostics.

Core competency set §6. Read ARIMA$(p,d,q)$; explain why multiplicative S-ARMA beats a giant plain ARMA on parameter count/stability; write the SARIMA operator.

---

7. Frequency domain

Picture to hold. A stationary process is a superposition of random periodic components at different frequencies. The time domain (ACF) and frequency domain (spectrum) hold the same information in two bases; you Fourier-transform to whichever is easier — exactly the MA↔AR spirit, now for the whole series.

Why transform: optimizing $\cos(2\pi ft+\phi)$ over a phase $\phi$ is non-convex, but the identity $\cos(x+y)=\cos x\cos y-\sin x\sin y$ rewrites it as a linear combination of $\cos$ and $\sin$ (amplitudes are linear coefficients). Euler $e^{ix}=\cos x+i\sin x$ packages this in complex form.

Discrete Fourier transform and its inverse

$$b_j = \sum_{t=0}^{n-1} x_t\,e^{-2\pi i jt/n}\quad(\text{contribution of frequency } j/n),\qquad x_t = \frac1n\sum_{j=0}^{n-1} b_j\,e^{2\pi i jt/n}.$$

Inversion proof (geometric sum / Euler). Substitute $b_j$ into the inverse:

$$\frac1n\sum_{j=0}^{n-1}\Big(\sum_{s=0}^{n-1}x_s e^{-2\pi i js/n}\Big)e^{2\pi i jt/n} = \frac1n\sum_{s=0}^{n-1}x_s\sum_{j=0}^{n-1}\big(e^{2\pi i (t-s)/n}\big)^{j}.$$

The inner sum: if $t=s$ every term is $1$, total $n$; if $t\neq s$, geometric series

$$\sum_{j=0}^{n-1}\big(e^{2\pi i (t-s)/n}\big)^{j} = \frac{1-e^{2\pi i (t-s)}}{1-e^{2\pi i (t-s)/n}} = 0$$

(numerator $0$ since $e^{2\pi i\cdot\text{integer}}=1$). So the inner sum is $n\,\mathbf 1(t=s)$, and the whole thing is $\frac1n\cdot n\, x_t = x_t$. $\checkmark$

The Fourier frequencies $j/n$ give exact reconstruction; the basis vectors $u^{(j)}$ are orthogonal (same geometric-sum argument on $\langle u^{(l)},u^{(k)}\rangle$). The $b_j$ are typically sparse — a few frequencies carry the signal — which is why this is great for compression.

Periodogram and spectral density

Periodogram (sample power at frequency $j/n$):

$$I(j/n) = \frac{|b_j|^2}{n},\qquad j=1,\dots,\lfloor n/2\rfloor$$

(only to $n/2$ by symmetry $b_{n-j}=\bar b_j$). Divide by $n$ because white-noise coefficients scale like $\sqrt n$, so $|b_j|^2/n$ keeps noise at a constant floor while true signals grow with $n$.

White-noise periodogram is not flat-looking. It shows random spikes, and a fresh sample gives entirely different spikes — white noise is an equal-weight (high-variance) superposition of all frequencies. Don’t over-read individual periodogram peaks (the §1 lesson, frequency-side).

Time↔frequency duality (sample side). The periodogram is the DFT of the sample autocovariance:

$$I(j/n) = \sum_{h=-(n-1)}^{n-1}\hat\gamma(h)\,e^{-2\pi i jh/n}.$$

(Proof: expand $|b_j|^2=b_j\bar b_j$ after centering — legal since the constant-signal DFT is $0$ — collect the double sum by lag $h=t-s$, recognize $\hat\gamma(h)$.)

Spectral density (population limit): for ACVF $\gamma(h)$,

$$f(\lambda) = \sum_{h=-\infty}^{\infty}\gamma(h)\,e^{-2\pi i\lambda h},\qquad -\tfrac12\le\lambda\le\tfrac12,$$

the periodogram with true $\gamma$ and continuous $\lambda$. It is invertible:

$$\gamma(h)=\int_{-1/2}^{1/2} e^{2\pi i\lambda h}f(\lambda)\,d\lambda = \int_{-1/2}^{1/2}\cos(2\pi\lambda h)f(\lambda)\,d\lambda.$$

Inversion sketch (continuous orthogonality). Substitute the definition of $f$ and swap sum and integral:

$$\int_{-1/2}^{1/2} e^{2\pi i\lambda h}f(\lambda)\,d\lambda = \sum_{k}\gamma(k)\int_{-1/2}^{1/2} e^{2\pi i\lambda(h-k)}\,d\lambda.$$

The integral is the continuous analogue of the DFT orthogonality used above: $\int_{-1/2}^{1/2}e^{2\pi i\lambda m}d\lambda = \mathbf 1(m=0)$ (for integer $m$ it’s $\frac{\sin(\pi m)}{\pi m}=0$ unless $m=0$). So only the $k=h$ term survives, leaving $\gamma(h)$. The $\cos$ form follows because $\gamma$ is symmetric, killing the imaginary part. ∎ (Same orthogonality move as the DFT inversion, run in the continuum.)

So the spectral density and the ACVF carry identical information — switch freely. White noise: $\gamma(h)=\sigma^2\mathbf 1(h=0)$ ⇒ $f(\lambda)=\sigma^2$, flat (equal power at all frequencies — and clearly not a probability density).

Linear filters

An LTI filter $Y_t=\sum_j a_j X_{t-j}$ (differencing, smoothing, lagging are all special cases). Its effect on the spectrum is multiplicative:

$$f_Y(\lambda) = f_X(\lambda)\,|A(\lambda)|^2,\qquad A(\lambda)=\sum_j a_j e^{-2\pi i j\lambda}\ (\text{transfer function}),\ |A|^2=\text{power transfer}.$$

Derivation. From $\gamma_Y(h)=\sum_{j,k}a_j a_k\gamma_X(h-k+j)$, write each $\gamma_X$ as a spectral integral and factor:

$$\gamma_Y(h)=\int_{-1/2}^{1/2} e^{2\pi i\lambda h} f_X(\lambda)\Big(\sum_j a_j e^{2\pi i\lambda j}\Big)\Big(\sum_k a_k e^{-2\pi i\lambda k}\Big)d\lambda = \int e^{2\pi i\lambda h} f_X(\lambda)|A(\lambda)|^2 d\lambda,$$

and the integrand multiplying $e^{2\pi i\lambda h}$ must be $f_Y$. $\checkmark$

Worked filter. Lag-$s$ differencing $Y_t=X_t-X_{t-s}$: $|A(\lambda)|^2=4\sin^2(\pi s\lambda)$. For $s=1$ this is $0$ at $\lambda=0$ (kills constants) and rises toward high $\lambda$ — a high-pass filter. For $s=6$ it has zeros at $\lambda=0,\tfrac16,\tfrac26,\tfrac36$ — it removes exactly the period-6 seasonal frequencies.

Spectral density of ARMA. For causal ARMA,

$$f_X(\lambda)=\sigma^2\frac{|\theta(e^{-2\pi i\lambda})|^2}{|\phi(e^{-2\pi i\lambda})|^2}.$$

As $\lambda$ traces the unit circle, near a root of $\phi$ the denominator → $0$ and $f$ peaks (AR resonances); near a root of $\theta$ the numerator → $0$ and $f$ dips (MA nulls). This is the frequency-domain echo of “roots govern everything.”

Aliasing. Sampling too slowly folds high frequencies onto low ones (the helicopter-blade-stops illusion at 24 fps). Sample fast enough or high frequencies masquerade as low — critical for audio.

Process representation — building a series with a target spectrum

The spectral density isn’t just a diagnostic; it tells you how to construct a process with prescribed frequency content. Represent

$$X_t = \sum_{j=1}^{m}\big(A_j\cos(2\pi\lambda_j t)+B_j\sin(2\pi\lambda_j t)\big),$$

with $A_j,B_j$ uncorrelated, mean zero, $\mathrm{Var}(A_j)=\mathrm{Var}(B_j)=\sigma_j^2$. This is stationary with $\gamma(h)=\sum_{j=1}^m\sigma_j^2\cos(2\pi\lambda_j h)$.

Why it approximates a target $f$ — derivation. Choose $\lambda_j=\tfrac{j}{2m}$ and $\sigma_j^2=\tfrac{f(\lambda_j)}{m}$. Then

$$\gamma(h)=\sum_{j=1}^m \tfrac{f(\lambda_j)}{m}\cos(2\pi\lambda_j h) = \tfrac{2}{2m}\sum_{j=1}^m f\big(\tfrac{j}{2m}\big)\cos\big(2\pi h\tfrac{j}{2m}\big).$$

As $m\to\infty$ this Riemann sum (step $\tfrac{1}{2m}$ over $[0,\tfrac12]$) becomes

$$\gamma(h)\to 2\int_0^{1/2} f(\lambda)\cos(2\pi h\lambda)\,d\lambda = \int_{-1/2}^{1/2} f(\lambda)\cos(2\pi h\lambda)\,d\lambda$$

(last step by symmetry of $f$ and $\cos$) — exactly the inversion formula. So $X_t$ has spectral density $f$ in the limit. ∎ This makes concrete the opening picture: a stationary process is a superposition of random sinusoids, weighted by the spectrum.

Implications of §7.

• ACVF ⟷ spectrum is a lossless change of basis; pick the view where the operation is simple (filtering is multiplication in frequency).
• The same geometric-sum identity underlies DFT inversion, basis orthogonality, and the duality theorems — one move, many corollaries.
• ARMA roots show up as spectral peaks/dips, tying §3’s algebra to what you see in a periodogram.

Core competency set §7.

• Prove DFT inversion via the geometric sum; sketch the continuous-orthogonality version for the spectrum.
• State the ACVF↔spectrum pair and the white-noise flat spectrum.
• Derive $f_Y=f_X|A|^2$ and read a differencing filter’s power transfer.
• Write the ARMA spectral density and explain peaks/dips from roots.
• Build a process with target spectrum via random sinusoids ($\sigma_j^2=f(\lambda_j)/m$) and show the Riemann-sum limit.

---

8. State space and hidden Markov models

Picture to hold. A latent process $X_t$ evolves; you see only noisy measurements $Y_t$. Want past/present/future of $X$ (smoothing/filtering/prediction). GPS: a true trajectory, scattered measurements — recover the path. This generalizes everything above and links to [[07_Bayesian_Statistics]] (EM) and [[05_Stochastic_Processes]] (Markov).

Hidden Markov model assumptions. $X$ is a Markov chain (future ⊥ past given present $X_t$); given $X$, the $Y_t$ are conditionally independent (measurement errors independent, not systematic). The joint factorizes:

$$f(X_0^T,Y_1^T) = f_0(X_0)\prod_{t=1}^{T} f_t(X_t\mid X_{t-1})\,g_t(Y_t\mid X_t).$$

Draw as a DAG: $X_{t-1}\to X_t\to Y_t$.

Linear state space model: $X_t=G_t X_{t-1}+V_t$, $Y_t=H_t X_t+W_t$ (position/velocity + noisy location). A causal/invertible ARMA is a state space model — stack lags into the state.

AR(2) as state space, worked. $Y_t=\phi_1 Y_{t-1}+\phi_2 Y_{t-2}+W_t$. State $X_t=(Y_t,Y_{t-1})^\top$: $\begin{pmatrix}Y_t\\Y_{t-1}\end{pmatrix}=\begin{pmatrix}\phi_1&\phi_2\\1&0\end{pmatrix}\begin{pmatrix}Y_{t-1}\\Y_{t-2}\end{pmatrix}+\begin{pmatrix}W_t\\0\end{pmatrix},\qquad Y_t=(1\ 0)X_t.$ This is the state augmentation move — bundle history into the state so the dynamics become first-order Markov. (Cross-link [[07_Bayesian_Statistics]] augment-then-marginalize.)

Filtering and prediction (the forward recursion)

Let $\pi^{t+k\mid t}_j=P(X_{t+k}=j\mid y_1^t)$. Prediction ($k>0$), filtering ($k=0$), smoothing ($k<0$). With transition matrix $A$ ($A_{ij}=P(X_{t+1}=i\mid X_t=j)$) and emission matrix $B$ ($B_{oj}=P(Y_t=o\mid X_t=j)$):

Prediction step: $\pi^{t+1\mid t}=A\,\pi^{t\mid t}$.

Update step — derive via Bayes. We have the model as $P(Y\mid X)$, $P(X\mid X)$; filtering needs the reverse condition. Using $P(A\mid B,C)=\frac{P(B\mid A,C)P(A\mid C)}{P(B\mid C)}$ with $A=\{X_{t+1}{=}j\}$, $B=\{Y_{t+1}{=}y_{t+1}\}$, $C=\{Y_1^t\}$, and dropping the observed past inside the emission (HMM conditional independence):

$$\pi^{t+1\mid t+1}_j = \frac{B_{y_{t+1},\,j}\,\pi^{t+1\mid t}_j}{\sum_i B_{y_{t+1},\,i}\,\pi^{t+1\mid t}_i}.$$

So the recursion alternates predict ($\pi^{t\mid t}\to\pi^{t+1\mid t}$, multiply by $A$) and update ($\pi^{t+1\mid t}\to\pi^{t+1\mid t+1}$, reweight by the emission and renormalize) — the HMM forward algorithm. Start from a prior $\pi^0$ (often uniform). Predicting $Y$: $p^{t+k\mid t}=B\,A^k\,\pi^{t\mid t}$.

Smoothing — the backward pass

Filtering uses data up to $t$; smoothing uses all data $y_1^T$ to estimate a past state $X_t$ ($k<0$) — better, because the future is informative about where you were. Define forward (filtered, unnormalized) messages $\alpha_t(j)=P(X_t=j,\,y_1^t)$ from the forward pass, and backward messages

$$\beta_t(j) = P(y_{t+1}^T\mid X_t=j),$$

the likelihood of all future observations given the current state. The backward recursion runs $T\to1$, conditioning on the next state and marginalizing it out:

$$\beta_t(j) = \sum_{i} A_{i,j}\,B_{y_{t+1},\,i}\,\beta_{t+1}(i),\qquad \beta_T(j)=1.$$

(Read right to left: from state $j$, transition to each next state $i$ with prob $A_{i,j}$, emit $y_{t+1}$ with prob $B_{y_{t+1},i}$, then carry the rest of the future $\beta_{t+1}(i)$.) The smoothed posterior multiplies the two passes:

$$P(X_t=j\mid y_1^T) \;\propto\; \alpha_t(j)\,\beta_t(j)$$

— past evidence ($\alpha$) times future evidence ($\beta$), renormalized. This forward–backward algorithm is exactly the E-step quantity EM needs below; it is the discrete cousin of the Kalman smoother.

Most likely trajectory — Viterbi

Find $\arg\max_{x_0^T} P(X_0^T\mid Y_1^T)$. Take logs; the normalizer $\log P(Y_1^T)$ is constant in $x$, so

$$\log P(X_0^T\mid Y_1^T) = \log\pi^0(x_0) + \sum_{t=1}^{T}\big[\log A_{x_t,x_{t-1}} + \log B_{y_t,x_t}\big] + C.$$

Read each term as a cost of traversing states over time: a start cost, a transition cost per move, an emission cost per stay. Minimizing total cost = shortest path through the trellis. Viterbi does it in one forward pass (cheapest path to each state at each time) + a backward trace — $O(T\cdot K^2)$ instead of $K^T$ brute force.

Parameter estimation — EM / Baum–Welch

When $A,B,\pi^0$ are unknown, the likelihood

$$P_\theta(Y_1^T)=\sum_{x_0^T}\pi^0(x_0)\prod_{t=1}^{T}A_{x_t,x_{t-1}}B_{y_t,x_t}$$

sums over exponentially many latent paths — non-convex, generally non-unique. EM (Baum–Welch) alternates: (E) estimate the latent $X$ distribution given $Y$ and current $\theta$ (the forward–backward recursions); (M) update $\theta$. Cross-link [[07_Bayesian_Statistics]] — same EM as for mixtures, here on a chain. Caveats: can diverge or stick in local optima (restart from random seeds). Hard EM replaces the E-step distribution with the single most likely path (Viterbi) — a Viterbi step + an MLE step.

Implications of §8.

• State augmentation unifies AR/ARMA with Kalman-style filtering — “stack the history” turns any finite-memory model Markov.
• Filtering is recursive Bayes: predict-then-update, the discrete Kalman filter; it’s why HMMs scale to long series.
• Viterbi (max path) vs forward (sum over paths) is the recurring max-vs-sum duality; both are dynamic programming over the trellis.

Core competency set §8.

• State the HMM factorization and write an AR(2) in state space form.
• Derive the filter update from Bayes; state the predict/update alternation.
• Write the backward recursion and the smoothing product $P(X_t\mid y_1^T)\propto\alpha_t\beta_t$ (forward–backward).
• Write the Viterbi log-decomposition and explain the DP cost-path reading.
• Say what EM/Baum–Welch and hard EM do and their failure modes.

---

9. Memorize cold

Must be instant:

• Weak stationarity: constant mean; $\mathrm{Cov}(X_t,X_{t+h})=\gamma(h)$ (lag only); $\gamma(-h)=\gamma(h)$; $\rho(h)=\gamma(h)/\gamma(0)$.
• WN band $\pm1.96/\sqrt n$; expected exceedances $0.05k$.
• $\nabla$ kills linear trend; $\nabla^2$ kills quadratic; $\nabla_d$ kills period-$d$ season.
• Variance-stabilize: $\mathrm{Var}\propto\mu\Rightarrow\sqrt{\cdot}$; $\mathrm{Var}\propto\mu^2\Rightarrow\log$.
• MA(q): $X_t=\theta(B)Z_t$; ACF cuts off after $q$; invertible $\iff\theta(z)\neq0$ on $|z|\le1$.
• AR(p): $\phi(B)X_t=Z_t$; PACF cuts off after $p$; causal $\iff\phi(z)\neq0$ on $|z|\le1$.
• AR(1): $\rho(h)=\phi^h$, $\gamma(h)=\sigma^2\phi^h/(1-\phi^2)$, $\mathrm{MA}(\infty)=\sum\phi^j Z_{t-j}$.
• Geometric-series inversion: $\frac{1}{1-\phi B}=\sum_{j\ge0}\phi^j B^j$.
• Yule–Walker: $\Gamma\Phi=\gamma$, $\Gamma$ symmetric Toeplitz.
• ACF cuts off → MA order; PACF cuts off → AR order.
• AIC ($+2k$) for prediction; BIC ($+k\log n$) for the true model.
• Time-series CV: train past, validate future only.
• Ljung–Box $\tilde Q=n(n+2)\sum \hat r_i^2/(n-i)\to\chi^2_{k-p-q}$.
• Spectral density $f(\lambda)=\sum_h\gamma(h)e^{-2\pi i\lambda h}$; WN ⇒ flat $f=\sigma^2$; filter ⇒ $f_Y=f_X|A|^2$; ARMA $f=\sigma^2|\theta|^2/|\phi|^2$.
• HMM filter: predict $\pi^{t+1\mid t}=A\pi^{t\mid t}$, update reweight-by-emission-and-renormalize.
• HMM smoothing: $P(X_t\mid y_1^T)\propto\alpha_t(j)\beta_t(j)$ (forward × backward).

Rederivable on demand: MA(1) stationarity, AR(1) uniqueness (remainder $\to0$), Yule–Walker from covariancing both sides + the AR(2) solve, the PACF-cutoff argument, DFT inversion, $f_Y=f_X|A|^2$, the process-representation Riemann-sum limit, the filter recursion from Bayes, the backward recursion, Viterbi log-decomposition.

Named moves (cross-track glossary)

• Invert-via-geometric-series — treat $\phi(B)$ as a scalar, expand $1/(1-\phi B)$; converts AR↔MA and divides ARMA polynomials. (§3, §4)
• Roots-outside-the-unit-circle — one condition, two names: causal (AR) / invertible (MA); also locates spectral peaks/dips. (§3, §7)
• Difference-to-stationary — remove trend/season without estimating it; one ∇ per polynomial degree, $\nabla_d$ per season; beware over-differencing. (§2)
• Time↔frequency duality — ACVF and spectrum are the same info; transform to make the operation easy (filtering = multiplication). (§7)
• State augmentation (“stack the history”) — bundle lags into a Markov state; turns ARMA into a state-space/Kalman problem. Mirrors augment-then-marginalize in [[07_Bayesian_Statistics]]. (§8)
• Predict-then-update — recursive Bayes filtering; the discrete Kalman filter / HMM forward pass. (§8)
• Delta method — variance of a transform via first-order Taylor; here for variance stabilization, same tool as in [[03_Statistical_Inference]]. (§2)
• Recognize-the-shape — match ACF/PACF cutoff-vs-decay to model order; the time-series version of the [[07_Bayesian_Statistics]] move. (§4)