Module 5: Model Diagnostics and Seasonality

Econ 6376 — Applied Time Series Econometrics

How to Use These Notes

This chapter has one main reading path and three optional layers.

Core material is the main text. It contains the concepts, notation, interpretations, and applied skills expected of everyone.
Deeper Dive sections explain why a result works or develop it more fully. They are useful, but they can be skipped on a first reading.
Technical Note sections state qualifications or formal details that matter for precise reasoning.
Looking Ahead sections introduce an idea that will be taught formally in a later module.

If you missed lecture, read the main text, run the Core code, and complete the Core Practice problems. Then return to the optional sections that address your questions or interests.

Core Learning Objectives

By the end of this module, you should be able to:

State the core principle of residual diagnostics: if the model is correct, residuals \(e_t\) behave like innovations \(\epsilon_t \sim \text{WN}(0, \sigma^2)\), and anything the residuals show that white noise shouldn’t is evidence that the model is wrong.
List the four things a clean residual diagnostic checks for: zero mean, no serial correlation, no obvious heteroskedasticity, approximate normality — and rank them by how much the course cares (serial correlation first, then heteroskedasticity as a preview of Module 14, then mean and normality).
Read a residual ACF/PACF pair against the Module 3 identification table and name the structure the model missed (extra AR, extra MA, seasonal, heteroskedastic).
Write down the Ljung-Box Q statistic, state its distribution under the null, and explain the degrees-of-freedom correction \(h - \text{fitdf}\), with fitdf = p + q for a non-seasonal ARMA fit and fitdf = p + q + P + Q for a SARIMA residual check, along with the choice of \(h\).
Run a Monte Carlo–style experiment that fits a non-seasonal ARMA to data from a seasonal DGP and show that the residual ACF and the Ljung-Box test catch the misspecification cleanly.
Articulate the framing — serial correlation in residuals is not a data property; it is a model property — and explain why that reframing is the whole justification for doing diagnostics.
Compare general-to-specific (Mizon, 1995) and specific-to-general modeling strategies, state why the course prefers general-to-specific in the presence of suspected serial correlation, and walk through one simplification step out loud.
Distinguish deterministic seasonality (seasonal dummies, sinusoids — a fixed pattern) from stochastic seasonality (seasonal differencing, seasonal ARMA — a pattern that drifts), and give one example of each from macro or weather data.
Apply seasonal differencing \(\Delta_s y_t = y_t - y_{t-s}\) and the combined \(\Delta \Delta_s y_t\) operator, and read seasonal ACF patterns at lags \(s, 2s, 3s, \ldots\) to distinguish stationary-with-seasonality from non-stationary-at-the-seasonal-frequency.
Write down the SARIMA\((p,d,q)(P,D,Q)_s\) specification, identify each block of the mixing console it lives on, and state the airline model ARIMA\((0,1,1)(0,1,1)_{12}\) as Box & Jenkins’s canonical monthly example.
Fit seasonal ARIMA candidates to UNRATE in R using forecast::Arima() and forecast::auto.arima(), re-run the Module 4 residual diagnostics on the seasonal fits, recognize the seasonal over-differencing red flag (\(\hat{\Theta}_1 \approx -1\)) on an already seasonally adjusted series, and verify that the lag-12 residual spike vanishes under the winning specification.

5.1 Where We Are on the Mixing Board

In Module 1 we wrote down the master equation for the whole course and immediately turned off most of it. Module by module we have been turning knobs back on. Module 3 brought the full AR(p) and MA(q) families to life; Module 4 showed how to pick an ARMA model using information criteria and the Box-Jenkins workflow. Through all of that, one block of the master equation has stayed silent: the seasonal block. In this module we turn it on.

Term	L3	L4	L5
\(\alpha\), \(\delta t\)	On	On	On
\(\phi_j y_{t-j}\) (all \(p\))	On	On	On
\(\theta_l \epsilon_{t-l}\) (all \(q\))	On	On	On
Seasonal block	Off	Off	On
\(\epsilon_t\)	On	On	On

By the end of this module, every knob on the master equation mixing board is on. Modules 6 onwards are about what you do with the full model — forecasting, multivariable extensions, VARs, volatility. The modeling vocabulary is essentially complete once we finish the seasonal band in this module.

But before we can turn on the seasonal block, we need to learn how to check a model — how to tell whether a fit is good enough or whether it is leaving structure on the table. That is the diagnostic toolkit. The module has two halves: Part A builds the diagnostic toolkit (residual analysis, the Ljung-Box Q test, the assumption-break), and Part B uses that toolkit to detect and fix the seasonal structure we have been quietly ignoring since Module 3. The two halves meet at the end, when a seasonal fit to UNRATE makes the diagnostic problem disappear.

5.1.1 The Residual-Spike Cliffhanger from L4

Module 4 ended on a cliffhanger. We fit an ARMA(1,2) to the differenced UNRATE series — the candidate that won every IC column in the Module 4 workflow — pulled its residuals, and plotted the residual ACF. There, sitting at lag 12, was a negative spike that clearly poked through the confidence bands, with an echo at lag 24. We said “that’s a clue, not a verdict, and we’ll deal with it next time.” This is next time.

Let us reproduce the fit from scratch so the residuals are on our screen. To pull the series live, store your FRED key in the environment (never hardcode it in course files) and run:

library(fredr)
fredr_set_key(Sys.getenv("FRED_KEY"))

unrate_raw <- fredr(series_id = "UNRATE",
                    observation_start = as.Date("1960-01-01"),
                    observation_end   = as.Date("2019-12-31"))
unrate   <- ts(unrate_raw$value, start = c(1960, 1), frequency = 12)
d_unrate <- diff(unrate)

So that these notes render identically for everyone — with or without an API key or an internet connection — the code below loads the same series from the course’s cached copy in data/UNRATE.csv (the same cache Modules 2–4 used) and applies the same 1960–2019 pre-COVID sample window Module 4 justified (the April 2020 outlier dominates every sample autocovariance; see Module 4, §4.7.1).

library(forecast); library(ggplot2); library(patchwork)
source("../helpers/simulators.R")   # arma_simulator() from Module 3

unrate_csv <- read.csv("../data/UNRATE.csv")  # cached FRED pull; see data/README.md
unrate_csv <- subset(unrate_csv,
                     observation_date >= "1960-01-01" &
                     observation_date <= "2019-12-01")

unrate   <- ts(unrate_csv$UNRATE, start = c(1960, 1), frequency = 12)
d_unrate <- diff(unrate)

# The candidate L4 picked: ARMA(1,2) on the differenced series
fit_l4 <- Arima(d_unrate, order = c(1, 0, 2))

e_hat <- residuals(fit_l4)

(autoplot(e_hat) + ggtitle("Residuals: ARMA(1,2) on diff(UNRATE)") + theme_bw()) /
  (ggAcf(e_hat, lag.max = 36) + ggtitle("Residual ACF") + theme_bw())

There it is again: the lag-12 spike. The residual ACF is mostly within the \(\pm 2/\sqrt{T}\) bands, but at lag 12 there is a clear, unambiguous breach — negative, around \(-0.15\) against a band of about \(\pm 0.075\) — with further breaches at lags 24 and 36. The model’s residuals are telling us something, and we did not yet have the vocabulary to listen. Now we do.

5.1.2 Two Questions, One Picture

The same plot motivates both halves of this module:

The diagnostic question: How do I know this model is wrong? What would I need to see for the residuals to pass? Is eyeballing the ACF enough, or is there a formal test? That is Part A.
The seasonality question: The thing the model is missing is clearly a 12-month cycle. What is the model class that has a 12-month cycle? How do I fit it, and how do I know it fixed the problem? That is Part B.

These are the same lesson in a trench coat. The diagnostic toolkit tells you the model is wrong; the seasonal toolkit tells you how to fix it.

5.2 What Residuals Are Supposed to Look Like

5.2.1 The Core Principle

Put the Module 1 notation slide back up and read it one more time:

\(\epsilon_t\) is the innovation — the true stochastic component of the DGP. Unobservable. Theoretical.
\(e_t\) is the residual — what you compute from a fitted model: \(e_t = y_t - \hat{y}_t\). Observable. Concrete.
If the model is correctly specified, \(e_t\) should behave like \(\epsilon_t\). That is, the residuals should look like a draw from white noise.
If they don’t, the model is wrong. There is nothing to fix in the data. There is something to fix in the model.

Everything in Part A is downstream of that one distinction. The diagnostic question is literally: do my residuals behave like white noise? Not are they white noise — they are a finite sample and no finite sample passes every test — but could they plausibly be?

Analogy — the post-game film room. Your prediction \(\hat{y}_t\) is the play you called. The residual \(e_t\) is what actually happened on the field minus what you predicted would happen. A good defensive coordinator does not care that the residual exists — the defense is going to do something unexpected. The coordinator cares about patterns in the residual. If every Tuesday’s game film shows the same unblocked blitz, you have a pattern. Patterns mean the playbook is missing a page.

5.2.2 The Four Things to Check

What does “behaves like white noise” mean concretely? Four properties, in the course’s priority order:

Zero mean. \(\bar{e} \approx 0\). Easy to check, easy to fix (add an intercept), almost never the real problem. Included for completeness.
No serial correlation. \(\hat{\rho}_k \approx 0\) for all \(k \geq 1\). This is the big one. The ACF of the residuals should be indistinguishable from the ACF of a white noise series. Every formal diagnostic test in this module is some version of this check.
No obvious heteroskedasticity. The conditional variance \(\text{Var}(e_t \mid \mathcal{F}_{t-1})\) should be roughly constant over time. Volatility clustering — quiet stretches followed by loud stretches — is a failure of this property. For Module 5 we eyeball this one.
Approximate normality. The weakest of the four. ARMA point estimates do not need Gaussian errors, but the standard errors and \(t\)-statistics do, asymptotically. Gross non-normality (heavy tails, bimodality) is a flag; mild non-normality is usually tolerable. A QQ plot is the tool.

The ranking matters. If you only have time to check one thing, check the residual ACF. If you have time for two, add the Ljung-Box. Heteroskedasticity and normality are plotting-level checks in this course until we hit Module 14.

Looking Ahead — Module 14: When the Variance Gets Its Own Model

Check number 3 is a preview of Module 14 (ARCH/GARCH), not a topic for this module. There, volatility clustering stops being a diagnostic nuisance and becomes the object of interest: the conditional variance itself gets a master equation. For Module 5 we eyeball; for Module 14 we test — and then we model.

5.2.3 Residual ACF Against the L3 Table

A residual ACF/PACF pair is a correlogram like any other — and the Module 3 identification table reads it. The content of the reading changes: in Module 3, the table told you what ARMA to fit. Here it tells you what ARMA structure your fit missed. Same tool, different question.

Residual pattern	Diagnosis
Clean — within the \(\pm 2/\sqrt{T}\) bands at every lag	Model passes the serial-correlation check
Spike at lag 1 only	Missing one more MA term (or one more AR term)
Geometric decay from lag 1	Missing AR structure
Cuts off at lag \(q^*\)	Missing \(q^*\) MA terms
Spikes at lags \(s, 2s, 3s\) (say \(s = 12\))	Missing seasonality — the Module 5 transition
Spikes that grow rather than shrink	Probably non-stationary residuals — refit with another difference

5.2.4 The Clean Case and the Broken Case, Side by Side

Before we bring in the formal test, let us see the difference between residuals from a correct fit and residuals from an underfit. Two fits, two diagnostic panels:

set.seed(1985)
# Clean case: simulate ARMA(1,1), fit ARMA(1,1)
y_clean <- arma_simulator(n = 500, phi = 0.6, theta = 0.4)
fit_clean <- Arima(y_clean, order = c(1, 0, 1))
e_clean <- residuals(fit_clean)

# Broken case: same series, fit AR(1) — deliberately underfit
fit_broken <- Arima(y_clean, order = c(1, 0, 0))
e_broken <- residuals(fit_broken)

(ggAcf(e_clean, lag.max = 24) + ggtitle("Residual ACF: correct fit") + theme_bw()) |
  (ggAcf(e_broken, lag.max = 24) + ggtitle("Residual ACF: missing MA term") + theme_bw())

Which one passes? The left one — the residual ACF sits within the bands at every lag, consistent with white noise. The right one has a spike at lag 1. That spike is the MA(1) term the model left on the table. If you had not known the truth and you only saw the right-hand ACF, the correct recommendation would be: add an MA term, refit, re-check. That is diagnostics — the loop.

5.3 The Ljung-Box Q Statistic

5.3.1 Why an Omnibus Test?

Eyeballing the residual ACF works, but it has two problems. First, the \(\pm 2/\sqrt{T}\) bands are a per-lag confidence band — at 20 lags you expect about one spike to poke through just by chance, and which lag it happens at changes the story you tell yourself. Second, the pattern-match against the identification table is a judgment call, and judgment calls do not belong in a homework grader. We need a single number that answers “are these residuals jointly consistent with white noise over the first \(h\) lags?” That number is the Ljung-Box \(Q\).

Think of \(Q\) as a joint \(F\)-test for the null hypothesis that all the first \(h\) residual autocorrelations are zero. If any of them are materially non-zero, \(Q\) gets large and the test rejects.

5.3.2 The Formula and Its Distribution

Ljung-Box statistic:

\[Q(h) = T(T+2) \sum_{k=1}^{h} \frac{\hat{\rho}_k^2}{T - k}\]

where \(\hat{\rho}_k\) is the sample autocorrelation of the residuals at lag \(k\), \(T\) is the sample size, and \(h\) is a user-chosen maximum lag.

Distribution under the null (residuals are white noise from a correctly specified ARMA or SARIMA fit):

\[Q(h) \stackrel{a}{\sim} \chi^2_{h - \text{fitdf}}\]

Two things to notice in that degrees-of-freedom line:

The dof is \(h\) minus the number of ARMA parameters you estimated, not just \(h\). You “spent” degrees of freedom fitting the model, and Ljung-Box accounts for that. Practically: in R, pass fitdf = p + q for a non-seasonal ARMA fit, or fitdf = p + q + P + Q for a SARIMA residual check.
If you are testing whether a raw series (not residuals, no model) is white noise, there is no fit, so fitdf = 0 and the dof is just \(h\).

5.3.3 Choosing \(h\)

Rules of thumb from the literature. The two most common:

Small-\(h\) rule: \(h \approx \log T\). Tidier for small samples. Misses longer-range structure.
Standard-\(h\) rule: \(h = 10\) for non-seasonal series, \(h = 20\) (or \(2s\)) for seasonal monthly series. Better power for seasonal detection.

For this course, use \(h = 10\) on non-seasonal residuals and \(h = 24\) on monthly seasonal residuals. Report both if you are in doubt. Do not go hunting through every \(h\) until one of them rejects — that is \(p\)-hacking with a chi-square.

Technical Note — Software Defaults for \(h\)

Box.test() has no default worth trusting (lag = 1) — always set lag and fitdf yourself. Hyndman’s checkresiduals() chooses \(h\) automatically: \(\min(10, T/5)\) for non-seasonal data and \(\min(2s, T/5)\) for seasonal data, with fitdf filled in from the model object. Those defaults agree with the course rule for the series sizes you will meet here, but when you report a Ljung-Box result, always say which \(h\) and which fitdf you used — two students can run “the” Ljung-Box test on the same residuals and get different \(p\)-values if they silently used different horizons.

5.3.4 Interpretation: The Inversion

Reading the \(p\)-value:

Large \(p\)-value: fail to reject. Residuals are plausibly white noise. This is the outcome you want — it means your model has handled the serial structure.
Small \(p\)-value: reject. Residuals are not plausibly white noise. Your model is missing something. Go look at the residual ACF/PACF to diagnose what.

This is an important inversion of the usual hypothesis-testing vibe. In regression, you want to reject the null (no effect). In diagnostics, you want to fail to reject the null (no structure). A large \(p\)-value on Ljung-Box is good news. Students coming from regression courses find this disorienting at first — lean into the disorientation, because it is the same conceptual shift that separates model-building from hypothesis-testing.

5.3.5 The R Workflow

Three ways to run the test, from most manual to most automated:

# Way 1: Box.test() with explicit arguments
# non-seasonal ARMA case: fitdf = p + q = 1 + 2 = 3
Box.test(e_hat, lag = 10, type = "Ljung-Box", fitdf = 3)


    Box-Ljung test

data:  e_hat
X-squared = 3.7698, df = 7, p-value = 0.8059

# Way 2: same, at the seasonal horizon
Box.test(e_hat, lag = 24, type = "Ljung-Box", fitdf = 3)


    Box-Ljung test

data:  e_hat
X-squared = 46.106, df = 21, p-value = 0.001237

# Way 3: checkresiduals() — the all-in-one from forecast
checkresiduals(fit_l4)


    Ljung-Box test

data:  Residuals from ARIMA(1,0,2) with non-zero mean
Q* = 46.106, df = 21, p-value = 0.001237

Model df: 3.   Total lags used: 24

checkresiduals() is the one-line diagnostic. It prints the model, the Ljung-Box result, and a three-panel plot: residual time series, residual ACF, and residual histogram with a normal overlay. For quick work it is the single most useful function in the forecast package after auto.arima() itself.

5.3.6 Application to the L4 UNRATE Fit

Look carefully at the two Box.test() results above, because they disagree — and the disagreement is the lesson.

At \(h = 10\), the test passes comfortably (the \(p\)-value is around 0.8). The ARMA(1,2) has genuinely absorbed the short-run dynamics; over the first ten lags, these residuals are indistinguishable from white noise. If you had only run the non-seasonal check, you would have signed off on this model.

At \(h = 24\), the test fails (the \(p\)-value is around 0.001). Once the test window is wide enough to see lags 12 and 24, the joint statistic picks up exactly the structure the residual ACF has been showing us since Section 5.1. This is why the course rule says to check the seasonal horizon on monthly data: a test that never looks at lag 12 cannot reject because of lag 12.

The formal test confirms what the eyeball already suspected: the model is leaving structure on the table. But what kind of structure? The residual ACF does not show geometric decay from lag 1, nor a spike at lag 1 — it shows spikes at lags 12, 24, and 36. That pattern is not a missing AR or MA term at the base frequency. It is a missing seasonal component. No non-seasonal ARMA is going to handle a 12-month cycle; the problem is not that we should switch from \((1,2)\) to \((2,2)\), it is that we should extend ARMA to its seasonal generalization. That is Part B.

5.4 The Assumption-Break: Catching a Seasonal DGP with Diagnostics

This section is the direct analog of the Dickey-Fuller power simulation in Module 2 and the IC Monte Carlo in Module 4: before you trust a diagnostic tool on real data, watch it behave on data you generated. If it catches a known misspecification cleanly, you can trust it on cases where you do not know the truth.

Analogy — the stain on a freshly-washed shirt. Serial correlation in the residuals means the wash didn’t work — not that the stain is interesting in itself. The object of interest is the washing machine (the model), not the stain (the residual structure). If the stain is still there after the wash cycle, you do not study the stain; you fix the washing machine.

5.4.1 The Experiment

Simulate a seasonal AR(1) DGP — a process with no non-seasonal ARMA structure at all. The only structure is a purely seasonal one at lag 12:

\[y_t = \Phi_1 y_{t-12} + \epsilon_t, \qquad \Phi_1 = 0.7, \quad \epsilon_t \sim N(0, 1), \quad T = 500\]

In SARIMA shorthand this is SARIMA\((0,0,0)(1,0,0)_{12}\). Note the capital \(\Phi_1\) — the course convention uses uppercase Greek for seasonal coefficients and lowercase for non-seasonal, a distinction that becomes more important later in this module.

Now fit a non-seasonal ARMA(1,1) to it — exactly the model a student without the Module 5 toolkit would try, following the Module 4 recipe blindly. Pull residuals. Run the diagnostics.

set.seed(2005)

# Simulate a pure seasonal AR(1), lag 12, Phi_1 = 0.7
n <- 500
eps <- rnorm(n + 12)
y_sar <- numeric(n + 12)
for (t in 13:(n + 12)) {
  y_sar[t] <- 0.7 * y_sar[t - 12] + eps[t]
}
y_sar <- ts(y_sar[-(1:12)], frequency = 12)

# Fit a non-seasonal ARMA(1,1) — deliberately the wrong class
fit_wrong <- Arima(y_sar, order = c(1, 0, 1))

# Diagnose
e_wrong <- residuals(fit_wrong)

(autoplot(e_wrong) + ggtitle("Residuals: non-seasonal ARMA(1,1) on seasonal AR(1)") + theme_bw()) /
  (ggAcf(e_wrong, lag.max = 36) + ggtitle("Residual ACF") + theme_bw())

Box.test(e_wrong, lag = 24, type = "Ljung-Box", fitdf = 2)


    Box-Ljung test

data:  e_wrong
X-squared = 377.59, df = 22, p-value < 2.2e-16

What you see: the residual series looks mostly fine to the eye — no obvious trends, no exploding variance. But the residual ACF has a clean, obvious spike at lag 12 (and another at 24), and Ljung-Box at \(h = 24\) rejects with a \(p\)-value that is numerically zero. The diagnostic tools flag the misspecification even though the fit itself reports “perfectly fine” ARMA parameters with sensible standard errors.

And to close the loop: change the fit to the correct class and the residuals clean up completely.

# The correct class: SARIMA(0,0,0)(1,0,0)_12
fit_right <- Arima(y_sar, order = c(0, 0, 0),
                   seasonal = list(order = c(1, 0, 0), period = 12))
coef(fit_right)   # Phi_1-hat approx 0.7 — the DGP value

      sar1  intercept 
0.70303063 0.09951039

Box.test(residuals(fit_right), lag = 24, type = "Ljung-Box", fitdf = 1)


    Box-Ljung test

data:  residuals(fit_right)
X-squared = 14.714, df = 23, p-value = 0.9046

A note on the hand-written simulation loop. The line inside that for loop — y_sar[t] <- 0.7 * y_sar[t - 12] + eps[t] — is a literal transcription of the seasonal AR(1) DGP into code, exactly the same pedagogical move we made in Module 3 when we built arma_simulator(). When you write the DGP out this way, there is no mystery about what you simulated.

For the problem set (and for the seasonal ACF gallery in Section 5.7.3), the course helpers include sarima_simulator(), which generalizes this loop to arbitrary multiplicative seasonal ARMA orders. Here it is in full — the canonical copy lives in helpers/simulators.R:

sarima_simulator <- function(
    n = 500,
    alpha = 0,
    phi = numeric(0),
    theta = numeric(0),
    Phi = numeric(0),
    Theta = numeric(0),
    s = 12,
    sigma = 1,
    burn_in = 200 + 5 * s) {
  # Simulate a (stationary) multiplicative seasonal ARMA process:
  #   Phi_p(L) * Phi_P(L^s) * y_t = alpha + Theta_q(L) * Theta_Q(L^s) * eps_t
  # Integration (d, D) is NOT applied here — difference/cumsum outside
  # if you need a non-stationary seasonal DGP.

  phi   <- as.numeric(phi);  theta <- as.numeric(theta)
  Phi   <- as.numeric(Phi);  Theta <- as.numeric(Theta)

  # Multiply two lag polynomials given as coefficient vectors (constant first)
  polymul <- function(a, b) {
    out <- rep(0, length(a) + length(b) - 1)
    for (i in seq_along(a)) {
      idx <- i:(i + length(b) - 1)
      out[idx] <- out[idx] + a[i] * b
    }
    out
  }

  # Spread seasonal coefficients onto lags s, 2s, ... (zeros in between)
  seas <- function(coefs) {
    if (!length(coefs)) return(numeric(0))
    out <- rep(0, s * length(coefs))
    out[s * seq_along(coefs)] <- coefs
    out
  }

  # Expand the multiplicative polynomials into single AR / MA coefficient
  # vectors, respecting the course sign conventions:
  #   AR: (1 - phi_1 L - ...)(1 - Phi_1 L^s - ...)   [minus signs]
  #   MA: (1 + theta_1 L + ...)(1 + Theta_1 L^s + ...) [plus signs]
  ar_poly <- polymul(c(1, -phi),  c(1, -seas(Phi)))
  ma_poly <- polymul(c(1,  theta), c(1,  seas(Theta)))
  phi_full   <- -ar_poly[-1]   # implied AR coefficients at every lag
  theta_full <-  ma_poly[-1]   # implied MA coefficients at every lag

  p <- length(phi_full)
  q <- length(theta_full)
  total <- burn_in + n
  eps <- rnorm(total, 0, sigma)
  y <- numeric(total)

  for (t in seq_len(total)) {
    ar_part <- 0
    if (p > 0 && t > p) {
      ar_part <- sum(phi_full * y[(t - 1):(t - p)])
    }
    ma_part <- 0
    if (q > 0 && t > q) {
      ma_part <- sum(theta_full * eps[(t - 1):(t - q)])
    }
    y[t] <- alpha + ar_part + ma_part + eps[t]
  }

  ts(y[(burn_in + 1):total], frequency = s)
}

The engine is the same for-loop as arma_simulator(); the only new work is expanding the multiplicative polynomial products into one long AR vector and one long MA vector before the loop starts. Section 5.7.1 and the Deeper Dive that follows it explain why those products are the heart of the SARIMA specification.

5.4.2 The Three Lessons

The diagnostic tools work. When the residuals contain real structure, Ljung-Box catches it and the residual ACF shows you what kind of structure. No prior knowledge of the DGP is required — the evidence is sitting in the residuals.
The fit itself will not volunteer the problem. fit_wrong prints tidy parameter estimates, reasonable standard errors, an AIC, a log-likelihood. None of those numbers says “wrong class.” If you stop at IC and declare a winner without checking residuals, you end up with a confidently wrong model. The Module 4 checklist had residual diagnostics as Step 7 for this reason. This is the Module 1 spurious-regression Monte Carlo paying off in a new setting: there, ignored dependence in the series produced confident \(t\)-statistics on a meaningless regression roughly three times out of four; here, ignored dependence in the residuals produces confident output from a wrong-class model. Both times the standard printout lies, and only a dependence-aware check catches it.
Serial correlation in the residuals is not a property of the data — it is a property of your model relative to the data. The DGP \(y_t\) is perfectly well-behaved; it is fully determined by \(y_{t-12}\) and white noise. The serial correlation you see in the residuals is being generated by the fit, not by nature. Change the fit (as we just did, to SARIMA\((0,0,0)(1,0,0)_{12}\)) and the residuals are clean. The “problem” lived in the model, not the series.

This is the framing that matters most. Students arrive in this course thinking of serial correlation as a data problem — something the data has, like a disease. The reframe is that serial correlation in residuals is a model report card — something the model produces when it is wrong. The data is fine. The model is wrong. Fix the model.

“Serial correlation is not a property of the data — it is a property of your model relative to the data. It means your model is wrong.”

That is the whole justification for doing residual diagnostics. If you take away one sentence from Part A, take that one.

5.5 General-to-Specific Modeling (Mizon, 1995)

5.5.1 Two Strategies

Now that we have a diagnostic toolkit, we need a strategy for using it. When diagnostics fail, there are two ways to navigate the space of candidate models:

General-to-specific (Mizon, 1995): Start with a big model — more lags, more terms than you think you need. Estimate. Check residuals. If they are clean, simplify: drop insignificant terms one at a time, re-estimating and re-checking after each drop. Stop when further simplification would break the residuals or materially hurt the IC. The final model is the smallest specification that still passes diagnostics.

Specific-to-general: Start with the smallest plausible model (say, AR(1)). Estimate. Check residuals. If they are not clean, add a term — another lag, an MA term — based on what the residual ACF suggests. Re-estimate and re-check. Stop when residuals pass.

Both are legitimate. Both appear in the literature. Both converge to similar answers on well-behaved data.

5.5.2 Why the Course Prefers General-to-Specific

Mizon’s (1995) argument, paraphrased: when you start small, the residuals in the intermediate steps are contaminated by the terms you have not yet added. Those contaminated residuals can push you toward the wrong next term — you are essentially doing variable selection on misspecified residuals, and the misspecification biases the selection. Starting general and simplifying means your residuals are (usually) clean at every step, so any dropping decision is based on evidence that is not poisoned by omitted structure.

Practical rule of thumb for this course: when in doubt, start one size bigger than you think is necessary, check the residuals, then trim. For UNRATE, that looks like starting with a seasonal specification such as ARIMA\((1,1,2)(1,0,1)_{12}\), confirming diagnostics pass, and then trimming any term whose removal leaves the diagnostics intact. (Section 5.8 runs exactly this play and finds that nothing can be trimmed.)

5.5.3 Honest Caveats

General-to-specific is not a magic bullet. If “general” is already misspecified — for instance, if you start with a non-seasonal model when the data are seasonal — starting big does not help. You cannot trim your way to a correct class. Mizon is about navigating within a correctly-chosen class. The class question is what Part B is for.

5.5.4 Connection to L4’s IC Workflow

Module 4 taught an implicit specific-to-general approach: start with a grid, score with IC, pick the minimizer. That is a fine first pass. The Module 5 upgrade is: after picking the IC winner, check residuals, and if they fail, use the failure to guide the next candidate. The loop — candidate \(\rightarrow\) IC \(\rightarrow\) diagnostics \(\rightarrow\) respecify — is the full Box-Jenkins loop; Module 4 stopped at IC, Module 5 closes the loop.

5.6 Two Kinds of Seasonality

5.6.1 The Seasonal Band on the Mixing Console

Back to the master equation. The lectures so far have treated \(p\) and \(q\) as lags at the base frequency (lag 1, 2, 3, …). A seasonal process has a second frequency to keep track of — the seasonal frequency. For monthly data, that is lag 12. For quarterly data, lag 4. For daily data with a weekly cycle, lag 7.

The two-band mixing console analogy. Think of the SARIMA model as a mixing console with two bands. The non-seasonal band handles ARMA at the base frequency (lag 1, 2, 3, …) — the same knobs we have been turning since Module 3. The seasonal band handles the same kind of structure at the seasonal frequency (lag \(s\), \(2s\), \(3s\), …). Both bands have AR knobs and MA knobs; both bands have integration (differencing). SARIMA is what you get when both bands are in play. This is a natural extension of the mixing-board image from Modules 1, 3, and 4 — we are just adding a second row of faders.

Analogy — the tide under the weather. Monthly weather is what you notice day-to-day — the non-seasonal band. Underneath the weather there is a longer, more predictable rhythm that repeats each year — the tide. If you model the weather and ignore the tide, you keep being blindsided by something that is, in fact, perfectly predictable. You just were not looking for it at the right frequency. That lag-12 spike in the UNRATE residuals is the tide we have been ignoring.

5.6.2 Deterministic Seasonality

Definition: the seasonal pattern is a fixed, repeating function of the calendar. Every July looks the same; every Monday is a bit quieter than every Tuesday. The pattern does not evolve — it is the same shape year after year.

How to model it:

Seasonal dummies: add \(s - 1\) indicator variables (one per month except the reference month) to the regression. Captures any shape, needs \(s - 1\) parameters. For monthly data, that is 11 extra parameters — not parsimonious, but flexible.
Sinusoids (harmonic regression): \(\sin(2\pi t / s)\), \(\cos(2\pi t / s)\), and higher harmonics. Fewer parameters for smooth cycles; worse for jagged ones.
Fourier terms via forecast::fourier(): the tidy R interface to the sinusoid approach. You choose how many Fourier pairs to include; the function generates the appropriate sine and cosine regressors.

When is deterministic seasonality the right model? When the pattern really is fixed and exogenous to the series. Weather patterns. Retail holidays. Tax filing deadlines. Things driven by the calendar, not by the system’s own dynamics. If every December looks the same because Christmas is always in December, seasonal dummies work fine.

5.6.3 Stochastic Seasonality

Definition: the seasonal pattern drifts over time. There is still a rhythm at the seasonal frequency, but the shape of the rhythm — its amplitude, its phase, the relative heights of the peaks and troughs — evolves. Last year’s December is informative about this year’s December, but it is not identical.

How to model it:

Seasonal differencing: \(\Delta_s y_t = y_t - y_{t-s}\). If the seasonal rhythm is doing a random walk at the seasonal frequency, subtracting last year’s observation from this year’s removes it, the same way \((1 - L) y_t\) removes a unit root at the base frequency.
Seasonal ARMA: \(\Phi_P(L^s)\) and \(\Theta_Q(L^s)\) — AR and MA polynomials in the seasonal lag. These model the dependence at the seasonal frequency, with or without seasonal differencing first.
SARIMA: the full package — both kinds of seasonality handling (seasonal differencing and seasonal ARMA) combined with the non-seasonal block from Modules 3-4.

When is stochastic seasonality the right model? When the pattern evolves. Macro seasonality in labor and consumption series is the canonical example: the shape of the December bump in hiring is not the same in 2008 as in 1978. The cycle is there, but it drifts with the economy. Most macro monthly series that students in this course will encounter have stochastic rather than deterministic seasonality.

5.6.4 Testing Which Kind: The Dummy-Residuals Workflow

There is no clean single test for deterministic vs stochastic seasonality, which is an honesty point worth stating out loud. The practical workflow is:

Fit a seasonal-dummy regression. Look at the residuals.
If the residuals are clean (no seasonal structure left), the seasonality was deterministic — the dummies absorbed it.
If the residuals still have seasonal structure (ACF spikes at \(s, 2s, 3s\) that the dummies did not capture), the seasonality is at least partly stochastic — move to seasonal ARMA terms, and to seasonal differencing if the seasonal spikes refuse to decay.

Rule of thumb for this course: if the data are macro, monthly, and showing visible seasonal structure after differencing at the base frequency, go stochastic. That is the empirical regularity in most series students will encounter in problem sets.

5.7 SARIMA\((p, d, q)(P, D, Q)_s\)

5.7.1 The Full Lag-Polynomial Specification

The SARIMA\((p, d, q)(P, D, Q)_s\) model in lag-polynomial form:

\[\Phi_p(L) \, \Phi_P(L^s) \, (1 - L)^d (1 - L^s)^D \, y_t = \Theta_q(L) \, \Theta_Q(L^s) \, \epsilon_t\]

Read that equation left to right in plain English:

\(\Phi_p(L) = 1 - \phi_1 L - \ldots - \phi_p L^p\) — the non-seasonal AR polynomial. Lags 1 through \(p\). Same notation as Modules 3-4.
\(\Phi_P(L^s) = 1 - \Phi_1 L^s - \ldots - \Phi_P L^{Ps}\) — the seasonal AR polynomial, using capital \(\Phi\)’s. Lags \(s, 2s, \ldots, Ps\).
\((1 - L)^d\) — non-seasonal differencing, applied \(d\) times. The Module 2 operator.
\((1 - L^s)^D\) — seasonal differencing, applied \(D\) times. New as of this module.
\(\Theta_q(L) = 1 + \theta_1 L + \ldots + \theta_q L^q\) — the non-seasonal MA polynomial. Course convention: MA with plus signs.
\(\Theta_Q(L^s) = 1 + \Theta_1 L^s + \ldots + \Theta_Q L^{Qs}\) — the seasonal MA polynomial, capital \(\Theta\)’s.

Notice that the seasonal and non-seasonal polynomials multiply rather than add. That multiplicative structure is what makes SARIMA parsimonious: the cross-frequency dynamics come for free as products of low-order terms instead of being estimated as separate free parameters. The Deeper Dive below works one expansion by hand.

5.7.2 Notation Conventions

Lowercase \((p, d, q)\) for non-seasonal orders. Uppercase \((P, D, Q)\) for seasonal orders. Lowercase Greek (\(\phi, \theta\)) for non-seasonal coefficients. Uppercase Greek (\(\Phi, \Theta\)) for seasonal coefficients. The subscript \(s\) is the period: \(s = 12\) for monthly, \(s = 4\) for quarterly, \(s = 7\) for daily-with-weekly-cycle. When you see ARIMA\((1, 1, 1)(0, 1, 1)_{12}\), the first triplet describes the base-frequency band, the second triplet describes the seasonal band, and the subscript says “monthly.”

One collision to watch: when the seasonal MA order appears alone in prose, the notation dictionary writes it \(Q_s\) to keep it distinct from the Ljung-Box \(Q\) statistic from Part A. Inside the \((p,d,q)(P,D,Q)_s\) shorthand, the position makes it unambiguous.

Deeper Dive — The Multiplicative Structure: A Worked Example

The polynomials multiply, and it is worth seeing once, by hand, what that buys. Consider the simplest case with both bands active:

\[\text{ARIMA}(1, 0, 0)(1, 0, 0)_{12}\]

Its AR side is \(\Phi_p(L) \cdot \Phi_P(L^s) = (1 - \phi_1 L)(1 - \Phi_1 L^{12})\). Expand:

\[(1 - \phi_1 L)(1 - \Phi_1 L^{12}) = 1 - \phi_1 L - \Phi_1 L^{12} + \phi_1 \Phi_1 L^{13}\]

Three things to notice:

The expansion generates coefficients at lags 1, 12, and 13. The lag-13 coefficient is not a free parameter — it is the product \(\phi_1 \Phi_1\), forced by the multiplicative structure.
You get lag-12 and lag-13 dynamics for two parameters instead of 13. This is what makes SARIMA parsimonious relative to a long non-seasonal ARMA that tries to capture lag-12 dynamics through sheer lag depth.
The multiplicative factorization says: the base-frequency dynamics and the seasonal-frequency dynamics are tied together through a parsimonious product structure, so the lag-13 cross-term is constrained by the lower-order terms rather than added as a separate free parameter. That is why the model stays compact. If you need a more flexible seasonal shape, you can move to a less restrictive seasonal specification — but for most macro applications, the multiplicative structure works well.

This expansion is exactly what sarima_simulator() from Section 5.4.1 automates: polymul() multiplies the polynomial coefficient vectors, and the simulator runs the master-equation loop on the expanded coefficients.

5.7.3 Seasonal ACF Patterns: Reading the Second Fingerprint

The Module 3 identification table extends to the seasonal band. Two failure modes show up at the seasonal frequency, exactly paralleling the base frequency. The four-panel gallery below — all simulated, so we know the truth in every panel — shows the progression from seasonal structure to stationarity.

set.seed(6376)

# Panel 1: stationary seasonal AR — SARIMA(0,0,0)(1,0,0)_12, Phi_1 = 0.8
y_p1 <- sarima_simulator(n = 480, Phi = 0.8, s = 12)

# Panel 2: seasonal random walk — y_t = y_{t-12} + eps_t (Phi_1 = 1)
n2 <- 480; e2 <- rnorm(n2 + 12); y2 <- numeric(n2 + 12)
for (t in 13:(n2 + 12)) y2[t] <- y2[t - 12] + e2[t]
y_p2 <- ts(y2[-(1:12)], frequency = 12)

# Panel 3: the panel-2 series after one seasonal difference
y_p3 <- diff(y_p2, lag = 12)

# Panel 4: an airline-type DGP, (1 - L)(1 - L^12) y_t = eps_t,
#          after the combined difference Delta Delta_12
n4 <- 480; e4 <- rnorm(n4 + 13); y4 <- numeric(n4 + 13)
for (t in 14:(n4 + 13)) y4[t] <- y4[t - 1] + y4[t - 12] - y4[t - 13] + e4[t]
y_p4 <- diff(diff(ts(y4[-(1:13)], frequency = 12)), lag = 12)

((ggAcf(y_p1, lag.max = 40) + ggtitle("1: Stationary seasonal AR (Phi = 0.8)") + theme_bw()) |
 (ggAcf(y_p2, lag.max = 40) + ggtitle("2: Seasonal random walk (Phi = 1)") + theme_bw())) /
((ggAcf(y_p3, lag.max = 40) + ggtitle("3: Panel 2 after seasonal differencing") + theme_bw()) |
 (ggAcf(y_p4, lag.max = 40) + ggtitle("4: Airline-type DGP after combined differencing") + theme_bw()))

Panel 1: Stationary with seasonal structure (seasonal ARMA present, seasonal differencing not needed). The ACF has spikes at lags \(s, 2s, 3s\) that decay as you move to higher multiples (here roughly \(0.73, 0.51, 0.34\) — geometric decay in the seasonal lag). This is the seasonal analog of a stationary AR: the seasonal memory fades with distance. The PACF gives the usual seasonal AR/MA identification at those lags, but it is noisier and less reliable — for seasonal orders, lean on the ACF.

Panel 2: Non-stationary at the seasonal frequency (seasonal differencing needed). The ACF spikes at \(s, 2s, 3s\) do not decay — they stay near one even at large multiples. This is the seasonal analog of a random walk’s ACF failing to decay at the base frequency. The series has a unit root at the seasonal frequency, and the fix is the same as for a base-frequency unit root: difference it away.

Panel 3: After seasonal differencing \(\Delta_s\). Apply \(\Delta_s y_t = (1 - L^s) y_t = y_t - y_{t-s}\) and re-inspect the ACF. If the seasonal spikes are now decaying (or gone, as here — the seasonal random walk differences down to pure white noise), the seasonal differencing did its job. If they are still flat, apply another \(\Delta_s\) — but \(D > 1\) is rare in practice.

Panel 4: After combined differencing \(\Delta \Delta_s\). A series can have both a unit root at the base frequency and non-stationarity at the seasonal frequency. The fix is \(\Delta \Delta_s y_t = (1 - L)(1 - L^s) y_t\) — difference once at each frequency. This is what the airline model does. After both differences, the ACF is clean enough to read off any residual MA structure (here, none — the DGP’s innovations were white).

5.7.4 Seasonal Differencing and the Combined \(\Delta \Delta_s\)

The seasonal difference operator is a direct analog of the base-frequency difference operator from Module 2:

\[\Delta_s y_t = (1 - L^s) y_t = y_t - y_{t-s}\]

For monthly data with \(s = 12\), this subtracts January of last year from January of this year, February of last year from February of this year, and so on. If the seasonal pattern is doing a random walk — this year’s December is last year’s December plus a shock — the seasonal difference removes it.

The combined operator handles both kinds of non-stationarity at once:

\[\Delta \Delta_s y_t = (1 - L)(1 - L^s) y_t\]

Expand \((1 - L)(1 - L^s)\) for \(s = 12\):

\[(1 - L)(1 - L^{12}) = 1 - L - L^{12} + L^{13}\]

So \(\Delta \Delta_{12} y_t = y_t - y_{t-1} - y_{t-12} + y_{t-13}\). This is the operator that takes a series with both a unit root and a seasonal unit root and (hopefully) produces something stationary. It is the differencing structure that sits inside the airline model.

5.7.5 The Airline Model

Box & Jenkins (1970) built their canonical seasonal example around monthly international airline passenger data (the AirPassengers series in R, still living in datasets). They proposed:

\[\text{ARIMA}(0, 1, 1)(0, 1, 1)_{12}\]

Written out in full:

\[(1 - L)(1 - L^{12}) y_t = (1 + \theta_1 L)(1 + \Theta_1 L^{12}) \epsilon_t\]

Four ingredients. Difference once at the base frequency. Difference once at the seasonal frequency. Add one non-seasonal MA term. Add one seasonal MA term. That is it — two estimated parameters, two differences.

Why the airline model is the default:

Parsimony: two free parameters regardless of \(s\). Compare to a seasonal-dummy regression with 11 parameters or a long ARMA that tries to capture lag-12 dynamics through sheer lag depth.
Robustness: works well for a remarkable range of monthly series. Box and Jenkins found it on airline data; Hyndman’s forecast package uses it as a near-default starting point for seasonal monthly series; it keeps showing up in academic benchmarking.
Historical weight: when a colleague says “I fit an airline model,” they mean this exact specification. It is a shibboleth of applied time series.

If your monthly series looks seasonal and you need a first-pass SARIMA in a meeting in two minutes, fit the airline model. You will be right often enough to get through the meeting. One important qualifier, which Section 5.8 demonstrates the hard way: the airline model assumes the series needs a seasonal difference. On a series that has already been seasonally adjusted at the source, \(D = 1\) is one difference too many.

5.7.6 The Airline Model on `AirPassengers`

Before we return to UNRATE, let us close the equation-to-code loop on the canonical example:

autoplot(AirPassengers) +
  ggtitle("Monthly international airline passengers, 1949-1960") + theme_bw()

# The textbook uses a log transform first (variance stabilization)
log_ap <- log(AirPassengers)

# The airline model
fit_ap <- Arima(log_ap, order = c(0, 1, 1),
                seasonal = list(order = c(0, 1, 1), period = 12))
summary(fit_ap)

Series: log_ap 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.4018  -0.5569
s.e.   0.0896   0.0731

sigma^2 = 0.001371:  log likelihood = 244.7
AIC=-483.4   AICc=-483.21   BIC=-474.77

Training set error measures:
                       ME       RMSE        MAE        MPE      MAPE      MASE
Training set 0.0005730622 0.03504883 0.02626034 0.01098898 0.4752815 0.2169522
                   ACF1
Training set 0.01443892

checkresiduals(fit_ap)


    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
Q* = 26.446, df = 22, p-value = 0.233

Model df: 2.   Total lags used: 24

Exactly as advertised: two estimated parameters (\(\hat{\theta}_1 \approx -0.40\) and \(\hat{\Theta}_1 \approx -0.56\), both negative and both comfortably inside the invertibility region), clean residual diagnostics, and a Ljung-Box \(p\)-value above conventional thresholds. This is Box and Jenkins’s canonical example — the textbook answer comes out, and the diagnostics pass.

5.7.7 The R Workflow for SARIMA

The R tools are the same as Module 4, with a seasonal argument added:

# Manual specification
Arima(y, order = c(p, d, q),
      seasonal = list(order = c(P, D, Q), period = s))

# Automatic search — auto.arima searches over both non-seasonal
# AND seasonal orders when the frequency of the ts() object is > 1
auto.arima(y)           # searches seasonal automatically for monthly/quarterly ts
auto.arima(y, seasonal = FALSE)  # force non-seasonal search only

A brief aside on forecast::seasadj(): this function takes an stl() decomposition and returns the seasonally-adjusted component. It is a descriptive tool — useful when you want to look at an underlying trend without the seasonal cycle — not a model-fitting tool. You may encounter it on the problem set; it is not the primary focus of this module.

5.8 Closing the Loop on UNRATE

This is the four-step closer. Take the real series from Section 5.1, fit seasonal candidates, watch the lag-12 spike disappear, confirm diagnostics pass. Every thread in this module — diagnostics, seasonality, SARIMA — lands here. And, this being real data, it lands with a twist that no simulation in this module prepared you for. That is what Step 4 is for.

5.8.1 Propose Candidates

We have the residual ACF of the Module 4 non-seasonal fit showing spikes at lags 12, 24, and 36, and Ljung-Box failing at \(h = 24\). Time to propose seasonal candidates. Following Section 5.7.5’s advice, the airline model goes on the list first — always propose, always fit. The other hand-fit candidates keep the non-seasonal ARMA(1,2) block that won Module 4 and turn on the seasonal band in increasing sizes, per the general-to-specific instinct. And we let auto.arima() weigh in without human guidance.

#	Model	Non-seasonal	Seasonal	Parameters
1	Airline	ARIMA\((0,1,1)\)	\((0,1,1)_{12}\)	2
2	ARIMA\((1,1,2)(1,0,0)_{12}\)	ARIMA\((1,1,2)\)	\((1,0,0)_{12}\)	4
3	ARIMA\((1,1,2)(0,0,1)_{12}\)	ARIMA\((1,1,2)\)	\((0,0,1)_{12}\)	4
4	ARIMA\((1,1,2)(1,0,1)_{12}\)	ARIMA\((1,1,2)\)	\((1,0,1)_{12}\)	5
5	`auto.arima(unrate)`	(selected automatically)	(selected automatically)	varies

Notice the deliberate contrast in the seasonal band: candidate 1 seasonally differences (\(D = 1\)); candidates 2–4 leave \(D = 0\) and model the seasonal dependence with seasonal ARMA terms instead. Section 5.6.3 said both are tools for stochastic seasonality — this table is how you let the data pick between them.

5.8.2 Fit and Compare

fit_air  <- Arima(unrate, order = c(0, 1, 1),
                  seasonal = list(order = c(0, 1, 1), period = 12))
fit_sar  <- Arima(unrate, order = c(1, 1, 2),
                  seasonal = list(order = c(1, 0, 0), period = 12))
fit_sma  <- Arima(unrate, order = c(1, 1, 2),
                  seasonal = list(order = c(0, 0, 1), period = 12))
fit_full <- Arima(unrate, order = c(1, 1, 2),
                  seasonal = list(order = c(1, 0, 1), period = 12))
fit_auto <- auto.arima(unrate)

fit_auto     # what did the algorithm pick?

Series: unrate 
ARIMA(2,1,2)(1,0,2)[12] with drift 

Coefficients:
         ar1      ar2      ma1     ma2    sar1     sma1     sma2    drift
      1.1717  -0.2581  -1.1927  0.4346  0.3809  -0.6243  -0.0919  -0.0014
s.e.  0.1758   0.1674   0.1643  0.1346  0.1322   0.1316   0.0613   0.0077

sigma^2 = 0.02477:  log likelihood = 311.61
AIC=-605.23   AICc=-604.97   BIC=-564.03

Start with the algorithm’s verdict: auto.arima() picks a SARIMA with seasonal AR and seasonal MA terms — and no seasonal difference (\(D = 0\)). Hold that thought while we look at the airline model:

coef(fit_air)

        ma1        sma1 
 0.06768543 -0.99999280

Look at \(\hat{\Theta}_1\): it is \(-1.000\), pinned exactly on the non-invertibility boundary. You have seen this red flag before — it is the Module 3 over-differencing signature (an MA coefficient driven to \(-1\) trying to undo one difference too many), and it is exactly what Step 6 of the Module 4 checklist told you to look for (“any \(\hat{\theta}\) stuck near \(\pm 1\)?”). The seasonal band is telling us that seasonally differencing UNRATE was one difference too many at the seasonal frequency.

Why? Read the FRED page for UNRATE: the series is published seasonally adjusted. The BLS has already removed the stable seasonal cycle at the source. What is left in the residuals of our Module 4 fit is the faint, drifting seasonal correlation that the adjustment filter did not fully absorb — real enough to fail Ljung-Box, far too mild to justify \(D = 1\). Force a seasonal difference onto an already-adjusted series and you create a seasonal unit root in the MA polynomial, which the estimator dutifully reports by slamming \(\hat{\Theta}_1\) to \(-1\). The airline model is the right default for a raw seasonal series like AirPassengers; it is the wrong tool for a seasonally adjusted one. Diagnostics caught that for us.

So the live comparison is among the \(D = 0\) candidates. The airline model also stays out of the IC comparison for a second, more general reason you need in your toolkit: information criteria only compare models of the same data, and a \(D = 1\) model’s likelihood is evaluated on the seasonally differenced series — a different, shorter series than the one the \(D = 0\) models use — so its AIC is not on the same scale. (The Technical Note below develops this.) Here is the IC table, with the Module 4 non-seasonal winner included as the baseline (refit on the levels with \(d = 1\) so all rows share the same differencing and their likelihoods are comparable):

fit_base <- Arima(unrate, order = c(1, 1, 2))   # the L4 winner, no seasonal band

ic_table <- data.frame(
  model = c("(1,1,2) non-seasonal [L4]",
            "(1,1,2)(1,0,0)[12]",
            "(1,1,2)(0,0,1)[12]",
            "(1,1,2)(1,0,1)[12]",
            "auto.arima"),
  AIC   = c(fit_base$aic,  fit_sar$aic,  fit_sma$aic,  fit_full$aic,  fit_auto$aic),
  AICc  = c(fit_base$aicc, fit_sar$aicc, fit_sma$aicc, fit_full$aicc, fit_auto$aicc),
  BIC   = c(fit_base$bic,  fit_sar$bic,  fit_sma$bic,  fit_full$bic,  fit_auto$bic)
)
ic_table[order(ic_table$AIC), ]

                      model       AIC      AICc       BIC
4        (1,1,2)(1,0,1)[12] -607.4004 -607.2825 -579.9333
5                auto.arima -605.2284 -604.9745 -564.0277
3        (1,1,2)(0,0,1)[12] -575.1249 -575.0407 -552.2356
2        (1,1,2)(1,0,0)[12] -564.2364 -564.1522 -541.3471
1 (1,1,2) non-seasonal [L4] -549.3327 -549.2767 -531.0213

Three readings from this table. First, the IC gap between the seasonal fits and the non-seasonal baseline is large — tens of AIC points — which is the quantitative evidence that the seasonal band is doing real work. Second, adding only half the seasonal band (candidates 2 and 3) buys some of the improvement, but the full seasonal ARMA(1,1) band buys much more: the drifting leftover seasonality apparently needs both a seasonal AR and a seasonal MA term to be captured. Third, our hand-built candidate 4 edges out auto.arima()’s pick on AIC and AICc, and beats it comfortably on BIC — the stepwise search landed nearby but spent three more parameters getting there. A very good intern; read its work.

Technical Note — Comparing IC Across Differencing Orders

The airline model is missing from that table on purpose. Information criteria compare models of the same data. A model with \(D = 1\) has its likelihood evaluated on the seasonally differenced series — a different (and shorter) dataset than the \(D = 0\) models use. Its AIC is not on the same scale, and sorting it into the table would be comparing the heights of mountains on different planets. The same warning applies across different \(d\). When candidates disagree about differencing, compare them on residual diagnostics and out-of-sample forecast performance (Module 6), not on IC. Within the table above, every row has \(d = 1, D = 0\), so the comparison is legitimate.

5.8.3 Diagnostics on the Winner

The IC winner is ARIMA\((1,1,2)(1,0,1)_{12}\) — on all three criteria, so there is no AIC-vs-BIC judgment call to adjudicate this time. Now the part Module 4 could not do: check it.

checkresiduals(fit_full)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,2)(1,0,1)[12]
Q* = 21.176, df = 19, p-value = 0.3272

Model df: 5.   Total lags used: 24

# By hand, at both course horizons, with fitdf = p + q + P + Q = 1 + 2 + 1 + 1 = 5
Box.test(residuals(fit_full), lag = 10, type = "Ljung-Box", fitdf = 5)


    Box-Ljung test

data:  residuals(fit_full)
X-squared = 4.0009, df = 5, p-value = 0.5493

Box.test(residuals(fit_full), lag = 24, type = "Ljung-Box", fitdf = 5)


    Box-Ljung test

data:  residuals(fit_full)
X-squared = 21.176, df = 19, p-value = 0.3272

The outcome: the lag-12 spike that has been haunting this module from the first section is gone. The residual ACF is within bands at the seasonal lags. The Ljung-Box \(p\)-value is comfortably above 0.05 at both \(h = 10\) and \(h = 24\). The residual histogram is roughly symmetric with modestly heavy tails — noted, tolerable (check four of Section 5.2.2 is the weakest priority).

One more general-to-specific beat before we declare victory: can we trim? Candidates 2 and 3 are the trimmed versions of the winner — each drops one seasonal term. Run the residual check on each (each trimmed fit estimates \(p + q + P + Q = 4\) ARMA parameters, so fitdf = 4):

# Candidate 2: drop the seasonal MA term -> (1,1,2)(1,0,0)[12]
Box.test(residuals(fit_sar), lag = 24, type = "Ljung-Box", fitdf = 4)


    Box-Ljung test

data:  residuals(fit_sar)
X-squared = 43.169, df = 20, p-value = 0.001942

# Candidate 3: drop the seasonal AR term -> (1,1,2)(0,0,1)[12]
Box.test(residuals(fit_sma), lag = 24, type = "Ljung-Box", fitdf = 4)


    Box-Ljung test

data:  residuals(fit_sma)
X-squared = 39.71, df = 20, p-value = 0.005434

Both fail Ljung-Box at \(h = 24\). Dropping either seasonal term breaks the residuals, so the simplification stops here. That is Mizon’s procedure executed to completion: the winner is the smallest specification that still passes.

5.8.4 Before and After

Put the two pictures side by side — the broken residual ACF from Section 5.1 (lag-12 spike) next to the clean residual ACF from the seasonal fit (no spike). Same series. Two images.

(ggAcf(residuals(fit_l4), lag.max = 36) +
   ggtitle("Before: ARMA(1,2) on diff(UNRATE)") + theme_bw()) |
  (ggAcf(residuals(fit_full), lag.max = 36) +
     ggtitle("After: ARIMA(1,1,2)(1,0,1)[12]") + theme_bw())

Everything we did in this module is between these two pictures. The left one asked the question; the right one answers it. When someone asks you what residual diagnostics are for, hand them this pair of plots.

5.8.5 Honest Caveats

Two of them, both worth internalizing.

“Clean” does not mean “perfect.” On a long macro series with financial crises and a structural break or two, you can usually get Ljung-Box to fail at some \(h\) if you look hard enough. The course’s standard is: at \(h = 10\) and \(h = 24\) with the right fitdf, using the default checkresiduals() call, the model passes. If you want residuals that are white noise at every conceivable lag and every conceivable subsample, you will not get one on macro data, and you should lower your standards to “residuals look roughly like white noise, model handles the main structure, anomalies are explicable.” That is what “model passes diagnostics” means in practice.

Know what your data has already been through. The single most consequential fact in this whole section — UNRATE is seasonally adjusted at the source — is written on the FRED series page, not in any correlogram. The diagnostics caught the consequence (\(\hat{\Theta}_1 = -1\)), but reading the data documentation would have predicted it. Raw seasonal series (retail sales not seasonally adjusted, AirPassengers) are where the full airline-model machinery, \(D = 1\) included, earns its keep.

Common Pitfalls and Misconceptions

“My data has serial correlation, so I need to fix the data.” Backwards. Serial correlation in residuals is a property of your model relative to the data. The data is fine; the model is missing structure. Fix the model — add the term the residual ACF points to.
“Small Ljung-Box \(p\)-value, great, significant!” The inversion. In diagnostics you want to fail to reject: a large \(p\)-value means the residuals are plausibly white noise. A small \(p\)-value means your model failed the check.
Forgetting fitdf. Testing ARMA residuals with fitdf = 0 uses the wrong null distribution and overstates the \(p\)-value. Use fitdf = p + q (non-seasonal) or p + q + P + Q (seasonal). On a raw series with no model, fitdf = 0 is correct.
Hunting across \(h\) until something rejects. That is \(p\)-hacking with a chi-square. Pick the course horizons (\(h = 10\) and \(h = 24\) for monthly data) before you look.
“The fit looks fine — good parameters, good standard errors, good AIC — so the model is fine.” The assumption-break showed a wrong-class model printing perfectly tidy output. Nothing in the estimation printout flags misspecification; only the residuals do. IC ranks candidates, diagnostics vet them.
Checking only the short horizon on monthly data. The UNRATE fit passed Ljung-Box at \(h = 10\) and failed at \(h = 24\). A test that never looks at lag 12 cannot reject because of lag 12.
Seasonally differencing everything monthly. If the series is already seasonally adjusted (most headline FRED macro series are — check the series page), \(D = 1\) over-differences at the seasonal frequency, and the estimator tells you by pinning \(\hat{\Theta}_1\) at \(-1\). Treat a seasonal MA estimate stuck at \(\pm 1\) as a red flag, not a parameter.
“Start general enough and you can always trim to the right model.” General-to-specific navigates within a correctly chosen class. If your “general” model is non-seasonal and the data are seasonal, no amount of trimming fixes the class error.
“Passing diagnostics means the model is true.” It means the model is adequate for the structure these tests can see, at these horizons, in this sample. That is the honest claim — and for forecasting purposes (Module 6), it is the claim that matters.

Connection to Enders

Residual diagnostics and model adequacy checking: Enders Chapter 2, pp. 80-90
The Ljung-Box Q statistic: Enders Chapter 2, pp. 82-84
Seasonal processes and SARIMA: Enders Chapter 2, pp. 90-105
The airline model: Enders Chapter 2, pp. 97-99
General-to-specific vs specific-to-general: Mizon (1995), “Progressive Modelling of Macroeconomic Time Series: The LSE Methodology”

A convention warning when you cross-reference. As in Modules 3 and 4, Enders writes AR coefficients as \(a_i\) where we write \(\phi_j\), and works in the characteristic-root-inside-the-unit-circle framing where the course states stationarity and invertibility as roots of \(\Phi(L)\) and \(\Theta(L)\) lying outside the unit circle — same condition, reciprocal roots. Hyndman & Athanasopoulos cover seasonal ARIMA from the practitioner’s angle (their Chapter 9), including checkresiduals() and the seasonal-differencing decision; note that their \((P,D,Q)_m\) notation uses \(m\) where we use \(s\) for the seasonal period. Box & Jenkins (1970) is the original source for both the diagnostic-checking step and the airline model.

Practice Problems

Core Practice

Deliberate underfit. Simulate an ARMA(1,1) with \(\phi = 0.6\), \(\theta = 0.4\), and \(T = 500\). Fit an AR(1) to it (deliberately underfit). Run checkresiduals() and report the Ljung-Box \(p\)-value. What does the residual ACF show? What would you add to the model? Explain your reasoning in terms of the L3 identification table applied to the residuals.
Seasonal misspecification. Simulate a SARIMA\((0,0,0)(1,0,0)_{12}\) with \(\Phi_1 = 0.7\) and \(T = 500\) (use sarima_simulator() or the hand loop from Section 5.4.1). Fit an ARMA(1,1) — deliberately the wrong class. Run the Ljung-Box test at \(h = 24\) with the correct fitdf. Report the \(p\)-value and explain which part of the residual ACF carries the smoking gun. Fit the correct SARIMA\((0,0,0)(1,0,0)_{12}\) and show that the residual ACF cleans up.
Monte Carlo rejection rate. Repeat Problem 2 as a Monte Carlo over 200 replications. In each replication: simulate the SARIMA\((0,0,0)(1,0,0)_{12}\) DGP, fit the wrong ARMA(1,1), run Ljung-Box at \(h = 24\), and record whether the test rejects at the 5% level. Report the fraction of replications in which Ljung-Box rejects. Write one paragraph interpreting the result: is this a “powerful” test against seasonal misspecification? How does this compare to the rejection rate you would expect under the null (correctly specified model)? In your paragraph, connect the exercise to the Module 1 spurious-regression Monte Carlo — both use a rejection rate to quantify what ignoring dependence costs your inference.
AirPassengers deep dive. Load AirPassengers. Log-transform. Fit the airline model ARIMA\((0,1,1)(0,1,1)_{12}\). Report \(\hat{\theta}_1\) and \(\hat{\Theta}_1\) and the Ljung-Box \(p\)-value. Then fit the more general ARIMA\((1,1,1)(1,1,1)_{12}\) and compare on AIC, BIC, and residual diagnostics. Which model do you prefer and why? Does the general-to-specific principle apply here?
UNRATE seasonal selection. For the FRED UNRATE series (1960–2019, as in Section 5.8), fit the following candidates: the airline model, ARIMA\((1,1,2)(1,0,0)_{12}\), ARIMA\((1,1,2)(0,0,1)_{12}\), and ARIMA\((1,1,2)(1,0,1)_{12}\). Also run auto.arima(unrate). Report the airline model’s \(\hat{\Theta}_1\) and explain, in two to three sentences, what it tells you and why it happens on this particular series. Build an IC table for the \(D = 0\) candidates (state why the airline model is excluded from it). Run checkresiduals() on the IC winner. Produce a before/after figure: the residual ACF of the Module 4 non-seasonal fit side-by-side with the residual ACF of your seasonal choice.
Conceptual. In two to three sentences, explain the phrase “serial correlation is a property of the model, not of the data.” Give a concrete example where this reframing changes the right course of action (i.e., what a student might do wrong if they think of serial correlation as a data property, and what they should do instead).

Key Takeaways

Residuals should behave like the innovations you assumed in the DGP. The \(\epsilon_t\) vs \(e_t\) distinction from Module 1 is the entire conceptual foundation of diagnostics. If \(e_t\) does not look like white noise, the model is wrong — there is nothing to fix in the data.
Serial correlation in residuals is a property of the model, not of the data. When residuals show structure, look at what structure, then change the model.
Ljung-Box is the omnibus check. \(Q(h) \stackrel{a}{\sim} \chi^2_{h - \text{fitdf}}\). Large \(p\)-value is good news. Use \(h = 10\) non-seasonal, \(h = 24\) seasonal, with fitdf set correctly — and check both horizons on monthly data.
General-to-specific beats specific-to-general when you suspect the initial model is too small. Start one size bigger; trim after diagnostics pass; stop before diagnostics break.
Seasonality comes in two flavors. Deterministic (fixed calendar effects, seasonal dummies or Fourier) and stochastic (drifting cycles, seasonal ARMA and/or seasonal differencing). Macro monthly series are usually stochastic — but check whether the series was already seasonally adjusted at the source before reaching for \(D = 1\).
SARIMA\((p, d, q)(P, D, Q)_s\) is the mixing console with both bands on. Non-seasonal block at the base frequency, seasonal block at lag \(s\), multiplicative polynomials.
The airline model ARIMA\((0,1,1)(0,1,1)_{12}\) is the default first swing for a raw monthly seasonal series. Two parameters, two differences, robust. On a seasonally adjusted series, its \(D = 1\) is one difference too many — and \(\hat{\Theta}_1 \to -1\) is how the fit tells you.
Every knob on the master equation mixing board is now on. The fitting toolkit built since Module 1 is complete.

Looking Ahead — Module 6

You have a model. You have checked it. Now the question is what you use it for.

Module 5 hands Module 6 a diagnostic-clean univariate model — the ARIMA\((1,1,2)(1,0,1)_{12}\) on UNRATE that passed Ljung-Box at both horizons. Module 6 (Forecasting Fundamentals) picks up exactly there:

Forecasting: the recursive computation of \(\hat{y}_{t+h \mid t}\), point forecasts, interval forecasts, density forecasts.
Forecast uncertainty: why the interval grows with the horizon, how far out is “too far” — uncertainty bands built from the \(\sigma^2\) we have been tracking since Module 1.
Train/test splits for time series: why random splits are wrong, why rolling windows are right.
The practical workflow: fit \(\rightarrow\) split \(\rightarrow\) forecast \(\rightarrow\) evaluate. Everything we did in this module was fit + check. Modules 6 and 7 are forecast + evaluate.

The modeling vocabulary is complete. From here on, every lecture is about what you do with a fitted model, or about lifting the univariate framework to multiple variables, or about relaxing the constant-variance assumption. The fitting toolkit we have been building since Module 1 is finished as of this module.