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_s$ for a SARIMA residual check ($P$ and $Q_s$ are the seasonal AR and MA orders, Section 5.7), 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 use the course’s three terms precisely: autocorrelation (a property of any series, measured by its ACF), serial correlation (the diagnostic verdict on a fitted model’s residuals), and residual seasonality (seasonal autocorrelation left in a published seasonally adjusted series — the adjustment procedure’s serial correlation).
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), place the seasonal unit root as the persistence boundary inside stochastic seasonality, explain why there is no clean “isn’t”, and apply the triage: deterministic → dummies; stochastic → SARIMA; “I need it gone, not modeled” → seasonal adjustment.
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.
Work the three treatment paths on the raw UNRATENSA series in R: seasonal dummies on the first difference (and read what their residual ACF leaves behind); SARIMA candidates with forecast::Arima() and forecast::auto.arima(), information criteria compared only among rows that share a differencing order, and diagnostics on the winner; and one seasonal-adjustment call whose output reproduces the published UNRATE.
State the differencing commitment from both sides: on an already adjusted series, $D = 1$ drives $\hat{\Theta}_1$ toward $-1$ (one seasonal difference too many); on a raw seasonal series, $D = 0$ drives $\hat{\Phi}_1$ toward $+1$ (one seasonal difference too few).
Describe seasonal adjustment at recognition level: what it does (trend-cycle + seasonal + irregular; remove the seasonal, publish the rest), who does it (statistical agencies running X-13ARIMA-SEATS), that most headline FRED macro series are its outputs, and why an adjusted series is a model’s output rather than the data.
On the published UNRATE, fit the seasonal candidates built on Module 4’s block, recognize the seasonal over-differencing red flag ($\hat{\Theta}_1 \approx -1$) under $D = 1$, verify that the lag-12 residual spike vanishes under ARIMA$(1,1,2)(1,0,1)_{12}$, and read what that model absorbed as residual seasonality — the model Module 6 inherits.

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 treat the seasonal structure we have been quietly ignoring since Module 3. Part B asks what seasonality is, walks three treatment paths — seasonal dummies, SARIMA, and seasonal adjustment — on the raw unemployment rate, and ends where Module 4 left off: on the published, seasonally adjusted UNRATE, with a seasonal fit that makes the diagnostic problem disappear and a name for what it absorbed.

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), include.mean = FALSE)  # Module 4's pick, same zero-mean policy

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())

# The seasonal-lag values behind the picture, and the band ggAcf() draws
acf_l4 <- acf(e_hat, lag.max = 36, plot = FALSE)$acf[-1]
round(c(lag12 = acf_l4[12], lag24 = acf_l4[24], lag36 = acf_l4[36],
        band = 1.96 / sqrt(length(e_hat))), 3)

 lag12  lag24  lag36   band 
-0.145 -0.125 -0.130  0.073

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, $\hat{\rho}_{12} = -0.145$ against the $\pm 0.073$ band R draws (ggAcf() uses $1.96/\sqrt{T}$ with $T = 719$; the course’s rule-of-thumb $\pm 2/\sqrt{T} \approx \pm 0.075$ is the same band to the eye) — with further breaches at lags 24 and 36 ($-0.125$ and $-0.130$). 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:

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.
Zero mean. $\bar{e} \approx 0$. Easy to check, easy to fix (add an intercept), almost never the real problem. Included for completeness.
Approximate normality. The weakest of the four. Gaussian innovations justify the default Gaussian likelihood and its model-based standard errors and $t$-statistics. Normality is not a universal requirement for ARMA large-sample inference: under suitable moment and regularity conditions, Gaussian quasi-MLE can remain consistent and asymptotically normal with non-Gaussian innovations, although robust covariance estimates may then be needed. Gross non-normality (heavy tails, bimodality) is still a flag because the software’s default model-based inference can be misleading; mild non-normality is often tolerable. A QQ plot is the first check.

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 2 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	First candidate diagnosis — not a verdict
Mostly within the $\pm 2/\sqrt{T}$ bands, with no systematic pattern	No visible evidence of remaining serial correlation at the inspected lags
Spike at lag 1 only	May suggest omitted short-run AR or MA structure
Geometric decay from lag 1	Suggests omitted low-order AR dynamics
Apparent cutoff after a few lags	Suggests omitted low-order MA dynamics
Spikes at lags $s, 2s, 3s$ (say $s = 12$)	Suggests omitted dependence at the seasonal frequency — the Module 5 transition
Large correlations that persist rather than decay	May indicate insufficient differencing or another source of non-stationarity; reassess the transformation

Treat the table as a way to propose the next candidate, not as automatic identification. A residual pattern can have more than one explanation in a finite sample. Propose a small number of plausible changes, refit them, compare parsimony and IC within a common model space, and then re-check the residuals. A diagnosis earns its place only if the revised fit removes the pattern without creating a new one.

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: apart from one marginal poke at lag 13 (about what chance gives across 24 lags; Section 5.3.1), the residual ACF sits inside the bands. The right one has a clear spike at lag 1 and a smaller, opposite-signed one at lag 2. Because we generated the data, we know those spikes come from the omitted MA(1) term. If you did not know the truth, the lag-1 spike would make an MA term a sensible first candidate, not a verdict: propose it alongside any plausible small AR alternative, refit, and re-check. That is diagnostics — the loop.

5.2.5 Three Terms the Course Keeps Apart

The word correlation is about to do three different jobs in this module, and the textbooks do not always tell you which one. The course draws a deliberate line. Learn it here, because every verdict in Part B is written in this vocabulary.

Term	Meaning	Belongs to
Autocorrelation (also dependence, memory, persistence)	A property of a series, measured by its ACF. Applies to any series — including residuals looked at as a series.	The data (or any series)
Serial correlation	The diagnostic verdict that a fitted model’s residuals $e_t$ show autocorrelation. A property of the model relative to the data. The course never says “the data has serial correlation.”	The model
Residual seasonality	Seasonal autocorrelation left in a published seasonally adjusted series. Under the rule above, it is the adjustment procedure’s serial correlation — the procedure is itself a model.	The adjustment procedure

The rule, in the instructor’s words:

“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.”

Read the three rows against the pictures you have already seen. The ACF of a raw series — UNRATE, or the seasonal AR(1) we will simulate in Section 5.4 — shows autocorrelation: that is what the series is like, and no model is involved. The lag-12 spike in the residuals of the Module 4 ARMA(1,2) is serial correlation: it is a report on that fit, and a different fit can make it disappear without touching the data. And when a spike at lag 12 survives in a series that a statistical agency has already seasonally adjusted, the course’s name for it is residual seasonality — a report on the agency’s procedure, which is a model like any other. That third row is the destination of Part B.

A textbook flag, stated once. Enders, Hamilton, and Hyndman & Athanasopoulos use autocorrelation and serial correlation interchangeably, and so does most applied practice. You should recognize the interchangeable usage in the wild. The course’s line is deliberate: it keeps “what the series is like” separate from “what your model failed to capture,” and that separation is the entire logic of diagnostics.

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_s$ 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. Its documented rule also makes the selected lag at least fitdf + 3, so the reference chi-squared distribution retains at least three degrees of freedom after fitted-model parameters are counted. 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.7597, df = 7, p-value = 0.807

# 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.115, df = 21, p-value = 0.001233

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


    Ljung-Box test

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

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.

Because Part B runs the two-horizon check on many fits, here is a small convenience function that reads fitdf off an Arima() object (its $arma slot stores $p, q, P, Q$ in that order) and reports both horizons at once. The canonical copy lives in helpers/diagnostics.R:

lb_both <- function(fit, h = c(10, 24)) {
  # Ljung-Box at each horizon in h on an Arima() fit's residuals, with
  # fitdf = p + q + P + Q read from the model object (fit$arma stores
  # p, q, P, Q, s, d, D in that order).
  fitdf <- sum(fit$arma[1:4])
  rows <- lapply(h, function(hh) {
    bt <- Box.test(residuals(fit), lag = hh, type = "Ljung-Box", fitdf = fitdf)
    data.frame(h = hh, fitdf = fitdf, Q = unname(bt$statistic), p = bt$p.value)
  })
  do.call(rbind, rows)
}

lb_both(fit_l4)   # the same two verdicts as Way 1 and Way 2 above

   h fitdf         Q           p
1 10     3  3.759659 0.807006394
2 24     3 46.115372 0.001233451

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 points to dependence at the seasonal frequency. A sufficiently long unrestricted AR or MA could mechanically place terms at lags 12, 24, and beyond, but that is not the parsimonious low-order, consecutive-lag ARMA search we have been conducting. The sensible next candidate is therefore a seasonal ARMA block rather than simply moving from $(1,2)$ to $(2,2)$. 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 — it has autocorrelation, which is what a time series is. 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.5 demonstrates this pruning step and finds that neither seasonal term can be trimmed within the tested reductions.

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 ran the Box-Jenkins checklist through Step 6: an IC-scored grid, then a sanity check on the winner. That is neither strategy, because it never looked at residuals. 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 deliberately stopped before Step 7, Module 5 closes the loop.

5.6 What Seasonality Is: Two Kinds, Four Traditions, Three Paths

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 Nobody Fully Agrees What It Is — But You Can See It

The two definitions above are the course’s working vocabulary. Now step back, because there is an honesty point that shapes everything in Part B: there is no single agreed-upon definition of seasonality. A century of applied statistics has produced at least four traditions, each with its own core object and its own representative papers:

Tradition	Core object	Representative work
Calendar means (regularity and repetition)	Seasonal means across a calendar partition — “every July is higher than the annual average”	Falkner (1924); Kallek (1978)
Seasonal unit roots and stability (the ARIMA polynomial)	Roots of the AR polynomial at the seasonal frequencies — “the seasonal pattern is a random walk”	Hylleberg, Engle, Granger & Yoo (1990); Canova & Hansen (1995)
Unobserved components (signal extraction)	A latent seasonal component $S_t$ identified by a restriction — “the series is trend-cycle plus seasonal plus irregular”	Hillmer & Tiao (1982); Bell & Hillmer (1984); X-13ARIMA-SEATS
Spectral (frequency domain)	A functional of the spectral density at the seasonal frequencies — “variance piles up at the annual cycle and its harmonics”	Nerlove (1964); Granger (1978)

Each tradition is internally consistent. Each answers a well-posed question. And none nests the others: a series can have seasonal unit roots and pass a calendar-means test, or show a clear spectral peak with no identifiable unobserved component. In practice, analysts consult several of them at once and quietly assemble a superset that nobody has written down. A line from the 1978 Census–NBER conference volume that tried to settle the question captures the situation:

“Some of us looking for answers regarding seasonal analysis may feel as if we’re in a dark room looking for a black cat. That’s bad, but it’s not as bad as it could be. Pity the philosophers who are in a dark room looking for a black cat that isn’t there.” — Zellner et al. (1978, p. 451)

The instructor’s standing position for this course: we do not fully agree on what seasonality is, but you can see it. Section 5.8.1 opens with a picture that makes the point — a raw monthly series whose ACF spikes at lags 12, 24, and 36 with no model fitted at all. Whatever definition you end up preferring, that picture is what every definition is trying to capture.

5.6.5 The Taxonomy: Where the Seasonal Unit Root Sits

The course’s working taxonomy has one more boundary than the two definitions above, and it is worth drawing explicitly:

Non-seasonal series.
Seasonal series, which split into:
- Deterministic seasonality (Section 5.6.2): fixed calendar effects.
- Stochastic seasonality (Section 5.6.3): drifting effects — and inside the stochastic box, at its extreme-persistence boundary, sit the seasonal unit roots.

The seasonal unit root is not a third species. It is the stochastic case pushed to the limit where the seasonal pattern never mean-reverts — exactly as the random walk is the AR(1) pushed to $\phi_1 = 1$ (Module 2). That is why the seasonal ACF gallery in Section 5.7.3 reads the way it does: decaying spikes at $s, 2s, 3s$ mean a stationary seasonal ARMA (mild persistence, inside the stochastic box); flat spikes mean a seasonal unit root (the boundary), and the fix is to difference at the seasonal frequency. The differencing decision $D$ is the decision about which side of that boundary the series is on, and Section 5.8.3 shows that the data push back from both sides when you get it wrong.

5.6.6 There Is No Clean “Isn’t”

Students want a checklist: “seasonality is X; calendar effects, trading-day variation, holidays, and outliers are not seasonality.” The course declines to give one, and the reason is worth understanding.

A definition is a line you draw. Once you have chosen a definition — seasonal means, seasonal unit roots, a latent component, a spectral peak — whatever lands inside that line is seasonal by construction, whether or not you would have called it seasonal beforehand. Trading-day variation is the standard example: the number of Mondays in a month varies with a rhythm close to (but not exactly at) the annual frequency, so under a spectral definition it leaks into the seasonal frequencies, and under a calendar-means definition it shows up as part of the monthly means. It is not that trading-day effects “are” or “are not” seasonality; it is that the definition decides, and different definitions decide differently. The same goes for moving holidays like Easter. These effects are not statistically separable from seasonality without an additional identifying assumption — which is precisely what X-13ARIMA-SEATS supplies when it pre-regresses them out (Section 5.8.4).

There is one safe non-example, and only one: monthly sampling alone does not make a series seasonal. A random walk observed every month has no seasonal structure of any kind, under any of the four definitions. Being monthly gives a series a seasonal frequency at which structure could live; it does not put structure there.

This is also why the Deeper Dive at the end of these notes (Section 5.9) starts by drawing a line explicitly — a partition of the frequency axis into “seasonal” and “non-seasonal” bins — and then asks a testable question relative to that line. The partition is the definition. It is declared before looking at the data, and everything that lands in a seasonal bin counts.

5.6.7 The Triage: Three Paths

Given a monthly series with visible seasonal autocorrelation, there are three things you can do about it. Which one you reach for depends on what you want the model for:

Diagnosis or need	Path	What it commits you to
The pattern is deterministic — fixed, exogenous, calendar-driven	Path 1: seasonal dummies (or Fourier terms) in a regression	A fixed shape. If the pattern drifts, the residuals will say so (Section 5.8.2)
The pattern is stochastic — it drifts, and you want to model it and forecast with it	Path 2: SARIMA — seasonal ARMA terms, with seasonal differencing if there is a seasonal unit root	A joint model of both bands of the mixing console (Section 5.8.3)
“I need it gone, not modeled” — you want a series you can read, compare across months, or feed to a non-seasonal model	Path 3: seasonal adjustment — decompose, remove the seasonal component, publish the rest	Trusting a procedure that is itself a model (Section 5.8.4)

The three paths are not rivals. Statistical agencies run path 3 to publish the series you read in the news; forecasters run path 2 on either the raw or the published series; and path 1 is the diagnostic you run first because it is the cheapest way to ask whether the pattern is fixed. Section 5.8 walks all three on one series.

5.6.8 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), that result is consistent with seasonal dummies being adequate for the observed sample. It is not proof that the underlying seasonality is permanently deterministic; re-check stability and out-of-sample performance when that distinction matters.
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.

This workflow is path 1 of the triage. Section 5.8.2 runs it on the raw unemployment rate, and the residuals are the lead-in to path 2.

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. Read the seasonal lags of the ACF and PACF as a pair, as in the Module 3 table: a seasonal AR(P) has ACF decay at s, 2s, 3s and a PACF that cuts off after lag Ps; a seasonal MA is the mirror image.

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 seasonal correlations remain, do not automatically apply a second $\Delta_s$. Renew the seasonal-unit-root diagnosis and consider seasonal ARMA structure, breaks, or changing seasonal patterns first. A second seasonal difference requires fresh evidence of another seasonal unit root; otherwise $D = 2$ risks seasonal over-differencing and an MA root at or near the unit circle.

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 a canonical benchmark:

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.
Benchmark value: it works well for many raw monthly series and gives a transparent two-parameter baseline against which richer seasonal models can be judged. auto.arima() does not treat the airline specification as a near-default; it chooses differencing orders using tests and searches across seasonal and non-seasonal ARIMA candidates under its configured limits.
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 a raw monthly series has evidence of both an ordinary and a seasonal unit root, the airline model is a useful first benchmark. It is not an automatic choice: its $d = 1$ and $D = 1$ both require evidence. Section 5.8 demonstrates both failure modes: on the raw unemployment rate the airline model’s seasonal difference is warranted but its non-seasonal block is too thin (Section 5.8.3), and on the series already seasonally adjusted at the source another seasonal difference is one difference too many (Section 5.8.5), so require renewed seasonal-unit-root evidence before setting $D = 1$.

5.7.6 The Airline Model on `AirPassengers`

Before we return to the unemployment rate, 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. The production tool statistical agencies use for seasonal adjustment is X-13ARIMA-SEATS, which Section 5.8.4 calls exactly once.

5.8 Closing the Loop: Three Paths on UNRATENSA, and an Epilogue on UNRATE

This is the four-step closer, and it has a twist built into its first line. Every UNRATE number in this course since Module 2 has come from a series that FRED publishes seasonally adjusted. There is a second series on FRED, UNRATENSA — the same civilian unemployment rate, not seasonally adjusted — and it is what the Bureau of Labor Statistics actually measures each month before any procedure touches it. This section takes the raw series, walks all three treatment paths from the triage (Section 5.6.7) on it, and then returns to the published series we have been modeling all semester — now able to say what the lag-12 spike from Module 4 was.

5.8.1 The Raw Series: You Can See It

Load UNRATENSA from the course cache exactly the way Section 5.1 loaded UNRATE, over the same 1960–2019 window (the cached file runs from 1948 to the present, with the same single release-gap NA in 2025-10 as UNRATE.csv; it is outside the window and never touched).

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

unratensa   <- ts(unratensa_csv$UNRATENSA, start = c(1960, 1), frequency = 12)
d_unratensa <- diff(unratensa)

dates <- seq(as.Date("1960-01-01"), by = "month", length.out = length(unratensa))

Now look at it — the level, and the ACF of its first difference (the transformation Module 2 settled on for the unemployment rate) out to lag 36. No model has been fitted.

(autoplot(unratensa) +
   ggtitle("UNRATENSA: civilian unemployment rate, not seasonally adjusted, 1960-2019") +
   ylab("Percent") + theme_bw()) /
  (ggAcf(d_unratensa, lag.max = 36) +
     ggtitle("ACF of diff(UNRATENSA) — no model fitted") + theme_bw())

# The seasonal-lag autocorrelations, and the band ggAcf() draws (1.96 / sqrt(T))
acf_nsa <- acf(d_unratensa, lag.max = 36, plot = FALSE)$acf[-1]
round(c(lag12 = acf_nsa[12], lag24 = acf_nsa[24], lag36 = acf_nsa[36],
        band = 1.96 / sqrt(length(d_unratensa))), 3)

lag12 lag24 lag36  band 
0.804 0.777 0.742 0.073

# For contrast: the same lag on the published (adjusted) series
round(acf(d_unrate, lag.max = 12, plot = FALSE)$acf[13], 3)

[1] -0.089

The sawtooth in the level plot is the annual hiring and layoff cycle — January and June jump, spring and fall fall — and it is visible without any statistics. In the ACF of the first difference, the spikes at lags 12, 24, and 36 are about 0.80, 0.78, and 0.74 against a band of $\pm 0.073$: ten times the band, and barely decaying from one year to the next. On the published UNRATE, the same lag-12 autocorrelation is $-0.09$.

Use the vocabulary from Section 5.2.5. What this picture shows is seasonal autocorrelation — a property of the raw series, measured by its ACF, with no model in sight. The lag-12 spike in the Module 4 residuals was serial correlation — a verdict on a fitted model. The two pictures look alike and mean different things, and Section 5.8.5 will connect them.

One more reading before we treat anything. Section 5.7.3 said decaying seasonal spikes mean a stationary seasonal ARMA and flat spikes mean a seasonal unit root. These spikes are close to flat. Hold that thought; it is the whole story of path 2.

5.8.2 Path 1: Seasonal Dummies

The cheapest thing to try is the deterministic model: eleven month dummies plus an intercept, on the first difference (consistent with the course’s treatment of the unemployment rate as $I(1)$ since Module 2). The regression assumes a fixed monthly pattern — the same January jump every year. If that assumption is right, the residuals will be clean at lags 12, 24, 36.

Building the dummy matrix is a three-line job worth having as a function. Here it is inline; the canonical copy lives in helpers/seasonality.R:

seasonal_dummies <- function(y, ref = 1) {
  # The (s - 1) seasonal indicator columns for a ts object y, with the
  # reference season `ref` dropped so a regression intercept is that
  # season's mean. Column names follow the calendar (Jan..Dec for s = 12,
  # Q1..Q4 for s = 4, s1..s_s otherwise).
  s      <- frequency(y)
  season <- as.numeric(cycle(y))
  labels <- if (s == 12) month.abb else if (s == 4) paste0("Q", 1:4) else paste0("s", seq_len(s))
  keep   <- setdiff(seq_len(s), ref)
  X <- sapply(keep, function(j) as.numeric(season == j))
  colnames(X) <- labels[keep]
  X
}

X_month <- seasonal_dummies(d_unratensa)          # reference month: January
fit_dum <- lm(d_unratensa ~ X_month)

# The fitted monthly pattern: mean change in the unemployment rate, by month
month_means <- c(coef(fit_dum)[1], coef(fit_dum)[1] + coef(fit_dum)[-1])
names(month_means) <- month.abb
round(month_means, 3)

   Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct    Nov 
 0.903 -0.068 -0.263 -0.517 -0.108  0.635 -0.157 -0.260 -0.183 -0.140  0.110 
   Dec 
 0.018

round(summary(fit_dum)$r.squared, 3)

[1] 0.721

ggplot(data.frame(month = factor(month.abb, levels = month.abb), effect = month_means),
       aes(month, effect)) +
  geom_col(fill = "steelblue") + geom_hline(yintercept = 0) +
  ggtitle("Path 1: the fixed monthly pattern in diff(UNRATENSA), 1960-2019") +
  ylab("Mean monthly change (pp)") + theme_bw()

The dummies absorb a lot: $R^2 \approx 0.72$ on a first-differenced series, with January averaging $+0.90$ percentage points, June $+0.64$, and April $-0.52$. Every one of the eleven dummies is many standard errors from zero. Now the diagnostic — the residuals, as a series, and their ACF:

e_dum <- ts(residuals(fit_dum), start = start(d_unratensa), frequency = 12)

ggAcf(e_dum, lag.max = 36) +
  ggtitle("Residual ACF: month dummies on diff(UNRATENSA)") + theme_bw()

acf_dum <- acf(e_dum, lag.max = 36, plot = FALSE)$acf[-1]
band    <- 1.96 / sqrt(length(e_dum))
round(c(lag12 = acf_dum[12], lag24 = acf_dum[24], lag36 = acf_dum[36]), 3)

lag12 lag24 lag36 
0.350 0.298 0.223

which(abs(acf_dum) > band)   # seasonal lags 12, 24, 36, plus low lags 2, 3, 5, 7, 11 and a few others

 [1]  2  3  5  7 11 12 13 18 23 24 25 30 33 35 36

# Ljung-Box on a regression residual: no ARMA parameters were fitted, so fitdf = 0
Box.test(e_dum, lag = 24, type = "Ljung-Box", fitdf = 0)


    Box-Ljung test

data:  e_dum
X-squared = 261.53, df = 24, p-value < 2.2e-16

The verdict, honestly. The dummies took out the fixed part and left substantial seasonal structure: the residual ACF at lags 12, 24, and 36 is about 0.35, 0.30, and 0.22 — three to five times the band — and Ljung-Box at $h = 24$ rejects with a $p$-value that is numerically zero. Two things about that test need saying. First, it rejects for non-seasonal reasons too: several low lags — 2, 3, 5, 7, 11 — also sit outside the band, because the first difference of the unemployment rate has the short-run ARMA dynamics Module 4 modeled, and a dummy regression has no ARMA terms at all. Second, therefore, the seasonal-lag values are the evidence to look at, not the omnibus $p$-value.

The clean way to isolate the seasonal verdict is to give the regression the non-seasonal dynamics Module 4 found and see what is left. Arima() takes the dummies as xreg and fits the ARMA(1,2) block around them:

fit_dum_arma <- Arima(d_unratensa, order = c(1, 0, 2), xreg = X_month)
round(coef(fit_dum_arma)[1:3], 3)   # the ARMA block: phi_1, theta_1, theta_2

   ar1    ma1    ma2 
 0.683 -0.757  0.213

e_dum_arma <- residuals(fit_dum_arma)
ggAcf(e_dum_arma, lag.max = 36) +
  ggtitle("Residual ACF: month dummies + ARMA(1,2) errors") + theme_bw()

acf_dum_arma <- acf(e_dum_arma, lag.max = 36, plot = FALSE)$acf[-1]
round(c(lag12 = acf_dum_arma[12], lag24 = acf_dum_arma[24], lag36 = acf_dum_arma[36]), 3)

lag12 lag24 lag36 
0.364 0.318 0.231

which(abs(acf_dum_arma) > band)                     # lags outside the band

[1]  6 12 18 23 24 25 33 35 36

round(max(abs(acf_dum_arma[-c(12, 24, 36)])), 3)    # the largest non-seasonal breach

[1] 0.108

# fitdf = p + q = 3 for the ARMA block; the dummies are regression coefficients
Box.test(e_dum_arma, lag = 10, type = "Ljung-Box", fitdf = 3)


    Box-Ljung test

data:  e_dum_arma
X-squared = 8.4022, df = 7, p-value = 0.2985

Box.test(e_dum_arma, lag = 24, type = "Ljung-Box", fitdf = 3)


    Box-Ljung test

data:  e_dum_arma
X-squared = 211.05, df = 21, p-value < 2.2e-16

Now the picture is much clearer. With the ARMA(1,2) errors, the first five lags are inside the band and Ljung-Box passes at $h = 10$ ($p \approx 0.30$) — the short-run dynamics are handled. It fails at $h = 24$ with $p \approx 0$, and the seasonal lags 12, 24, and 36 tower over everything else (about 0.36, 0.32, 0.23). Six other lags poke just past the band, at 6, 18, 23, 25, 33 and 35, the largest of them about 0.11 in magnitude, but the seasonal lags do most of the rejecting. That is the dummy-residuals workflow of Section 5.6.8 returning its answer: the fixed part is absorbed; a drifting part remains.

Is the pattern really drifting, or is this just a noisy fixed pattern? One line of evidence settles it. Fit the same month effects separately on 1960–1989 and 1990–2019 and test whether they are equal:

half  <- factor(time(d_unratensa) >= 1990, labels = c("1960-89", "1990-2019"))
month <- factor(month.abb[cycle(d_unratensa)], levels = month.abb)

fit_common <- lm(d_unratensa ~ month + half)     # same month effects in both halves
fit_split  <- lm(d_unratensa ~ month * half)     # month effects free to differ
anova(fit_common, fit_split)                     # F test of equal month effects

Analysis of Variance Table

Model 1: d_unratensa ~ month + half
Model 2: d_unratensa ~ month * half
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    706 39.773                                  
2    695 33.219 11    6.5544 12.466 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

signif(anova(fit_common, fit_split)[2, "Pr(>F)"], 3)   # the p-value behind "< 2.2e-16"

[1] 1.15e-21

# What moved: the June and July effects by half
by_half <- tapply(d_unratensa, list(month, half), mean)
round(by_half[c("Jun", "Jul"), ], 2)

    1960-89 1990-2019
Jun    0.84      0.43
Jul   -0.33      0.02

The month effects are not the same in the two halves — $F(11, 695) \approx 12.5$, $p \approx 10^{-21}$. The June jump roughly halves (from about $+0.84$ to $+0.43$ percentage points) and the July effect changes sign (from about $-0.33$ to $+0.02$). The seasonal pattern of the U.S. unemployment rate in 2019 is not the seasonal pattern of 1965. A fixed-dummy model cannot say that; a stochastic model can. On to path 2.

One line on Fourier terms: replacing the eleven dummies with a few sine–cosine pairs from forecast::fourier() is more parsimonious and gives the same verdict here, because parsimony is not the problem — fixedness is.

5.8.3 Path 2: SARIMA

Path 2 models the seasonal band instead of assuming it fixed. We take the residual evidence from path 1 — seasonal spikes that barely decay — into the SARIMA candidate table. Four candidates, chosen to make one decision visible: is there a seasonal unit root, so that $D = 1$?

#	Model	Non-seasonal	Seasonal	Why it is on the table
1	Airline	ARIMA$(0,1,1)$	$(0,1,1)_{12}$	The two-difference, two-parameter benchmark (Section 5.7.5)
2	ARIMA$(1,1,2)(1,0,1)_{12}$	ARIMA$(1,1,2)$	$(1,0,1)_{12}$	Module 4’s block plus a seasonal ARMA band, no seasonal difference
3	`auto.arima(unratensa)`	(selected automatically)	(selected automatically)	The intern’s answer, unguided
4	ARIMA$(1,1,2)(0,1,1)_{12}$	ARIMA$(1,1,2)$	$(0,1,1)_{12}$	The airline’s seasonal band on Module 4’s non-seasonal block

fit_air_nsa  <- Arima(unratensa, order = c(0, 1, 1),
                      seasonal = list(order = c(0, 1, 1), period = 12))
fit_d0_nsa   <- Arima(unratensa, order = c(1, 1, 2),
                      seasonal = list(order = c(1, 0, 1), period = 12))
fit_auto_nsa <- auto.arima(unratensa)
fit_win_nsa  <- Arima(unratensa, order = c(1, 1, 2),
                      seasonal = list(order = c(0, 1, 1), period = 12))

Candidate 1, the airline model. On the raw series the seasonal difference is warranted, and the coefficients say so:

round(coef(fit_air_nsa), 3)

   ma1   sma1 
 0.066 -0.777

lb_both(fit_air_nsa)

   h fitdf        Q            p
1 10     2 77.52260 1.538769e-13
2 24     2 96.63739 2.481870e-11

# Where the rejection comes from: the low lags, not the seasonal lags
acf_air <- acf(residuals(fit_air_nsa), lag.max = 36, plot = FALSE)$acf[-1]
round(c(lag2 = acf_air[2], lag3 = acf_air[3], lag4 = acf_air[4], lag5 = acf_air[5],
        lag12 = acf_air[12], lag24 = acf_air[24], lag36 = acf_air[36]), 3)

  lag2   lag3   lag4   lag5  lag12  lag24  lag36 
 0.214  0.139  0.119  0.112  0.076  0.024 -0.032

$\hat{\Theta}_1 \approx -0.78$ is comfortably interior — nowhere near the $-1$ boundary — so there is no over-differencing signature at the seasonal frequency, and the residual ACF at the seasonal lags is unremarkable (at the band edge at lag 12, well inside it at 24 and 36). Yet the airline model fails Ljung-Box at both horizons, with $p$-values around $10^{-13}$ and $10^{-11}$. Look at where the rejection comes from: lags 2 through 5 of the residual ACF are about 0.21, 0.14, 0.12, 0.11. The $(0,1,1)$ non-seasonal block is too thin for this series — $\hat{\theta}_1 \approx 0.07$ is doing almost nothing — and the short-run ARMA dynamics Module 4 found are sitting in the residuals. The airline model got the seasonal band right and the non-seasonal band wrong.

Candidate 2, the $D = 0$ seasonal ARMA. Keep Module 4’s $(1,1,2)$ block, add a seasonal AR and a seasonal MA term, and do not seasonally difference:

round(coef(fit_d0_nsa), 4)

    ar1     ma1     ma2    sar1    sma1 
 0.7857 -0.7509  0.1479  0.9939 -0.7483

lb_both(fit_d0_nsa)

   h fitdf          Q         p
1 10     5  0.9431223 0.9670217
2 24     5 14.3839657 0.7608656

# The seasonal AR polynomial 1 - Phi_1 z, with z = L^12, has its root at z = 1 / Phi_1
round(abs(1 / coef(fit_d0_nsa)["sar1"]), 4)

  sar1 
1.0061

This model passes Ljung-Box at both horizons ($p \approx 0.97$ and $0.76$), and, as the IC table below shows, it beats every other $D = 0$ candidate on all three criteria. It also has $\hat{\Phi}_1 \approx 0.994$. Read the seasonal AR polynomial $1 - \hat{\Phi}_1 L^{12}$ as a polynomial in $z = L^{12}$: its root sits at $z = 1/\hat{\Phi}_1 \approx 1.006$ — within 0.01 of the unit circle. That is a seasonal unit root in disguise: asked to model the seasonal band without differencing it, the estimator put the seasonal AR coefficient as close to one as it could. The data are asking for $D = 1$ in the only way a $D = 0$ model can ask.

Candidate 3, auto.arima(). The intern’s answer, at the time these notes were rendered:

fit_auto_nsa

Series: unratensa 
ARIMA(2,0,2)(2,1,2)[12] 

Coefficients:
         ar1      ar2      ma1     ma2     sar1    sar2     sma1     sma2
      1.8299  -0.8374  -0.8276  0.1594  -0.6649  0.1434  -0.0375  -0.6223
s.e.  0.0567   0.0557   0.0650  0.0425   0.2179  0.0583   0.2162   0.1563

sigma^2 = 0.04164:  log likelihood = 117.99
AIC=-217.98   AICc=-217.73   BIC=-176.92

lb_both(fit_auto_nsa)

   h fitdf        Q         p
1 10     8  2.64863 0.2659851
2 24     8 15.82172 0.4654712

# Smallest modulus among the roots of the non-seasonal AR polynomial (> 1 is stationary)
ar_nsa <- coef(fit_auto_nsa)[grep("^ar", names(coef(fit_auto_nsa)))]
round(min(Mod(polyroot(c(1, -ar_nsa)))), 3)

[1] 1.093

It chose $D = 1$ — and $d = 0$: auto.arima() seasonally differences first and then finds no further unit root in the seasonally differenced series. Eight coefficients, a non-seasonal AR root modulus around $1.09$ (outside the unit circle, so stationary after the seasonal difference, but not by a wide margin), and clean diagnostics at both horizons. A very good intern; read its work. Note that it now lives in a third differencing class, $(d, D) = (0, 1)$, which matters in a moment.

Candidate 4, the winner. Take what candidate 1 got right (the seasonal band, $D = 1$ with one seasonal MA term) and what candidates 2 and 3 got right (a non-seasonal block that can carry the short-run dynamics):

round(coef(fit_win_nsa), 3)

   ar1    ma1    ma2   sma1 
 0.783 -0.752  0.158 -0.753

lb_both(fit_win_nsa)

   h fitdf         Q         p
1 10     4  2.387679 0.8808203
2 24     4 20.546739 0.4242287

checkresiduals(fit_win_nsa)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,2)(0,1,1)[12]
Q* = 20.547, df = 20, p-value = 0.4242

Model df: 4.   Total lags used: 24

Every coefficient is interior ($\hat{\Theta}_1 \approx -0.75$), Ljung-Box passes at both horizons ($p \approx 0.88$ and $0.42$), and the residual ACF is inside the band at lags 12, 24, and 36 and at the low lags. And — a check we did not plan — X-13ARIMA-SEATS’s own automatic model selection, run in Section 5.8.4 with no guidance from us, chose exactly these orders for this series.

Information criteria, one table per differencing class. Information criteria compare models of the same data. A model with $D = 1$ evaluates its likelihood on the seasonally differenced series, which is a different and shorter series than the one a $D = 0$ model uses, and a model with $d = 0$ is different again. So the four candidates fall into three classes and three tables, and rows are comparable only within a table (the Technical Note in Section 5.8.5 develops the rule):

ic_row <- function(name, fit) data.frame(model = name, nobs = fit$nobs,
                                         AIC = fit$aic, AICc = fit$aicc, BIC = fit$bic)

# Class (d, D) = (1, 1): the airline model and the winner
rbind(ic_row("airline (0,1,1)(0,1,1)[12]", fit_air_nsa),
      ic_row("winner  (1,1,2)(0,1,1)[12]", fit_win_nsa))

                       model nobs       AIC      AICc       BIC
1 airline (0,1,1)(0,1,1)[12]  707 -155.7062 -155.6720 -142.0231
2 winner  (1,1,2)(0,1,1)[12]  707 -206.2327 -206.1471 -183.4275

# Class (d, D) = (1, 0): Module 4's block with the seasonal band turned on in stages
fit_base_nsa <- Arima(unratensa, order = c(1, 1, 2))
fit_sar_nsa  <- Arima(unratensa, order = c(1, 1, 2),
                      seasonal = list(order = c(1, 0, 0), period = 12))
fit_sma_nsa  <- Arima(unratensa, order = c(1, 1, 2),
                      seasonal = list(order = c(0, 0, 1), period = 12))
d0_table <- rbind(ic_row("(1,1,2) non-seasonal [L4 block]", fit_base_nsa),
                  ic_row("(1,1,2)(1,0,0)[12]", fit_sar_nsa),
                  ic_row("(1,1,2)(0,0,1)[12]", fit_sma_nsa),
                  ic_row("(1,1,2)(1,0,1)[12]", fit_d0_nsa))
d0_table[order(d0_table$AIC), ]

                            model nobs        AIC       AICc        BIC
4              (1,1,2)(1,0,1)[12]  719 -189.77765 -189.65967 -162.31048
2              (1,1,2)(1,0,0)[12]  719   39.31666   39.40081   62.20597
3              (1,1,2)(0,0,1)[12]  719  512.59669  512.68085  535.48600
1 (1,1,2) non-seasonal [L4 block]  719  844.67357  844.72959  862.98501

# Class (d, D) = (0, 1): auto.arima() sits alone
ic_row("auto.arima", fit_auto_nsa)

       model nobs       AIC      AICc       BIC
1 auto.arima  708 -217.9839 -217.7261 -176.9219

In the $(1, 1)$ class the winner beats the airline model by about 50 AIC points and 41 BIC points — the whole gap is the non-seasonal block. In the $(1, 0)$ class the seasonal ARMA(1,1) band wins on every criterion, by hundreds of points over Module 4’s non-seasonal block, which is quantitative evidence that the seasonal band is doing real work — and its winner is the model with $\hat{\Phi}_1 = 0.994$. auto.arima()’s AICc is not comparable to either table. So the tables cannot pick between candidates 2, 3, and 4; diagnostics and the coefficients do, and they pick candidate 4: it passes both horizons, every coefficient is interior, and it does not hide a unit root inside an AR polynomial.

Put the three residual ACFs side by side, because the pictures carry the argument:

(ggAcf(residuals(fit_air_nsa), lag.max = 36) + ggtitle("Airline: low lags fail") + theme_bw()) |
  (ggAcf(residuals(fit_d0_nsa), lag.max = 36) + ggtitle("D = 0: clean, but Phi = 0.994") + theme_bw()) |
  (ggAcf(residuals(fit_win_nsa), lag.max = 36) + ggtitle("Winner (1,1,2)(0,1,1)[12]") + theme_bw())

$D$ is a commitment, and the data push back from both sides. This is the paragraph to remember from path 2. Section 5.8.5 shows what happens when you impose $D = 1$ on a series that has already been seasonally adjusted: the seasonal MA coefficient is driven to $\hat{\Theta}_1 \approx -1$, the over-differencing signature from Module 3, because the model is trying to undo a difference the data did not need. Here, on the raw series, you have just seen the mirror image: leave $D = 0$ on a series that does have a seasonal unit root and the seasonal AR coefficient is driven to $\hat{\Phi}_1 \approx +1$, because the model is trying to manufacture a difference the data did need. One seasonal difference too many pins $\hat{\Theta}_1$ at $-1$; one too few pins $\hat{\Phi}_1$ at $+1$. Either boundary estimate is the data telling you that the differencing decision, not the ARMA orders, is what to revisit — and neither shows up in an IC table, because IC cannot compare across differencing orders. Read the coefficients.

5.8.4 Path 3: Seasonal Adjustment, at Recognition Level

Path 3 answers a different need. A forecaster wants to model the seasonal band; a statistical agency, a journalist, or a policymaker wants it gone, so that this month’s number can be compared with last month’s without the January jump getting in the way. That is seasonal adjustment, and you need to recognize it for one reason above all: most of the headline macro series you will ever download — the unemployment rate, payrolls, retail sales, GDP, CPI — are its outputs.

What it does. Seasonal adjustment starts from the unobserved-components tradition of Section 5.6.4: the observed series is treated as the sum (or product) of three latent pieces,

\[y_t = \text{trend-cycle}_t + \text{seasonal}_t + \text{irregular}_t,\]

estimates the seasonal piece with a model, subtracts it, and publishes the rest. The published seasonally adjusted series is trend-cycle plus irregular. The seasonal piece is what the procedure decided was seasonal — under its definition, with its identifying restriction.

Who does it. Statistical agencies — the Bureau of Labor Statistics, the Census Bureau, the Bureau of Economic Analysis, and their counterparts abroad — run a program called X-13ARIMA-SEATS, maintained by the Census Bureau. It fits a SARIMA model (a regARIMA model, since it also carries regressors for trading-day, leap-year, and holiday effects and for outliers), uses that model to extend the series at both ends, and then extracts the seasonal component either by moving-average filters (X-11) or by a model-based decomposition (SEATS). Every option has a default; agencies tune them per series. That is all you need to know about the internals in this course.

A note on the four-step rhythm. Path 3 has no simulate step, and that is deliberate: seasonal adjustment enters this course at recognition level — you should recognize the procedure and its outputs when you meet them, not build one — so we go from concept straight to one real-data call.

The one call. The R package seasonal wraps X-13ARIMA-SEATS, and its default call is the program’s default specification — automatic model selection, automatic transformation choice, automatic outlier and calendar-effect detection. Run it once on the raw series, with no options:

install.packages(c("seasonal", "x13binary"))   # x13binary ships the Census Bureau program

library(seasonal)
fit_x13 <- seas(unratensa)
summary(fit_x13)


Call:
seas(x = unratensa)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
Leap Year         -0.0869031  0.0366671  -2.370 0.017785 *  
Weekday           -0.0001555  0.0014314  -0.109 0.913489    
LS1975.Jan         1.2033565  0.1819286   6.614 3.73e-11 ***
AR-Nonseasonal-01  0.8264050  0.0553491  14.931  < 2e-16 ***
MA-Nonseasonal-01  0.8126629  0.0635875  12.780  < 2e-16 ***
MA-Nonseasonal-02 -0.1490694  0.0397547  -3.750 0.000177 ***
MA-Seasonal-12     0.7485219  0.0238592  31.372  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SEATS adj.  ARIMA: (1 1 2)(0 1 1)  Obs.: 720  Transform: none
AICc: -247.1, BIC: -210.8  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 18.53   Shapiro (normality): 0.9938 **

Read the printout for what it chose, not for the numbers. The automatic model selection landed on ARIMA$(1,1,2)(0,1,1)$ — the same orders as our path-2 winner, reached with no guidance from us — with a SEATS decomposition and no log transform (the additive form; the automatic test rejected the log). Two calendar regressors survived the AIC tests, leap year and weekday (the latter estimated at essentially zero), and the outlier search inserted one level shift at January 1975. Nothing was specified by hand. This is the tool a statistical agency would start from.

The reveal. Take the adjusted series the call produced — final(fit_x13) — and lay it over the published UNRATE we loaded in Section 5.1:

sa_x13 <- final(fit_x13)

cmp <- data.frame(date = dates,
                  raw       = as.numeric(unratensa),
                  published = as.numeric(unrate),
                  x13       = as.numeric(sa_x13))

ggplot(cmp, aes(date)) +
  geom_line(aes(y = raw, colour = "UNRATENSA (raw)"), linewidth = 0.4, alpha = 0.5) +
  geom_line(aes(y = published, colour = "UNRATE (published, seasonally adjusted)"), linewidth = 1) +
  geom_line(aes(y = x13, colour = "seas(unratensa): one default X-13 call"),
            linewidth = 0.7, linetype = "22") +
  scale_colour_manual(values = c("UNRATENSA (raw)" = "grey55",
                                 "UNRATE (published, seasonally adjusted)" = "darkorange",
                                 "seas(unratensa): one default X-13 call" = "navy"), name = NULL) +
  ggtitle("Path 3: one X-13 call on the raw series, against the published series") +
  ylab("Percent") + xlab(NULL) + theme_bw() + theme(legend.position = "bottom")

diff_x13 <- cmp$x13 - cmp$published
round(c(correlation   = cor(cmp$x13, cmp$published),
        mean_abs_diff = mean(abs(diff_x13)),
        max_abs_diff  = max(abs(diff_x13)),
        share_within_0.1pp = mean(abs(diff_x13) <= 0.1)), 4)

       correlation      mean_abs_diff       max_abs_diff share_within_0.1pp 
            0.9985             0.0647             0.3726             0.8028

format(cmp$date[which.max(abs(diff_x13))], "%Y-%m")

[1] "1976-02"

The two adjusted lines are one line. Correlation 0.9985; a mean absolute difference of about 0.065 percentage points; four months in five within 0.1 of the published value; a maximum gap of 0.37 points in February 1976, next to the level shift the program inserted. The differences are what you would expect from a default call against a production run — BLS adjusts concurrently, with its own model span, its own regressors, and a revision history the default call knows nothing about — so the honest word is that one X-13 call reproduces the published series, not that it equals it.

Now say the sentence out loud. The UNRATE series you have used since Module 2 is the output of this procedure. Every ACF, every ADF test, every ARMA fit, every information criterion in Modules 2 through 4 was computed on trend-cycle-plus-irregular from an X-13 run at the Bureau of Labor Statistics, not on what the survey measured. That is not a complaint — it is why the series is usable — but it changes what you are modeling.

The cost. Adjustment is a model. It fits a SARIMA, imposes an identifying restriction, and removes what that restriction calls seasonal. Like any model, it can be wrong relative to the data, and the course’s rule from Section 5.2.5 applies to it exactly as it applies to your own fits: the adjustment procedure can have serial correlation of its own. When that serial correlation shows up at the seasonal lags of the published series, the course’s name for it is residual seasonality — seasonal autocorrelation left in a series that has already been adjusted. Three consequences follow. The published series carries the procedure’s choices as well as the data’s dynamics, and a model fitted to it is partly modeling the procedure. The field’s tests and its adjustment filters answer different questions — a diagnostic from one of the traditions in Section 5.6.4 can flag seasonality that a filter built in another tradition was never designed to remove, and the two do not speak the same language. And you cannot diagnose residual seasonality without saying what seasonality is — a test for “seasonality left over” needs a definition of seasonality to test against, which is exactly what Section 5.6.4 said the field has never agreed on. The Deeper Dive in Section 5.9 is one way of drawing that line so that it can be tested. First, though, the epilogue.

5.8.5 Epilogue — Residual Seasonality on the Published UNRATE

We now return to the series Module 4 handed us — the published, seasonally adjusted UNRATE — and to the lag-12 spike in the residuals of its ARMA(1,2). Everything in this subsection is the adjusted-series analysis that closes the module and hands Module 6 its model; read it in the vocabulary Sections 5.2.5 and 5.8.4 just built. 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.

Propose Candidates on the Adjusted Series

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. We include the airline model as a canonical benchmark, while recognizing that its $D = 1$ is a substantive seasonal-unit-root choice rather than a default. The other hand-fit candidates keep the non-seasonal ARMA(1,2) block that won Module 4 and turn on the seasonal band: the full band (candidate 4) is the general model, and candidates 2 and 3 are its one-term trims, which the general-to-specific check below tests. 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.

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, allowmean = FALSE, allowdrift = FALSE, approximation = FALSE)

fit_auto     # what did the algorithm pick?

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

Coefficients:
         ar1      ma1     ma2    sar1     sar2     sma1
      0.9057  -0.9369  0.2187  0.4784  -0.0756  -0.7194
s.e.  0.0267   0.0441  0.0396  0.0867   0.0495   0.0794

sigma^2 = 0.02476:  log likelihood = 310.83
AIC=-607.67   AICc=-607.51   BIC=-575.62

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 $-0.9999928$, effectively at the non-invertibility boundary. You have seen this red flag before — it is the Module 3 over-differencing signature (an MA coefficient driven near $-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$?”). For UNRATE, this is strong evidence that the imposed seasonal difference is one difference too many at the seasonal frequency, and it should send us back to the $D = 1$ decision rather than forward to another seasonal difference. Compare Section 5.8.3, where the same model on the raw series gave an interior $\hat{\Theta}_1 \approx -0.78$: the seasonal difference that was warranted on UNRATENSA has already been taken, in effect, by the adjustment procedure.

Why reconsider $D = 1$? Read the FRED page for UNRATE: the series is published seasonally adjusted. That documentation makes an additional seasonal difference something to justify, not assume. It does not identify the exact source of the remaining residual dependence at lags 12, 24, and 36. That pattern could reflect remaining seasonal dynamics, features of the source adjustment and revision process, calendar effects, structural change, or misspecification elsewhere in the baseline model; the residual ACF alone cannot separate those explanations. What the fitted comparison does establish is narrower and sufficient for this workflow: imposing $D = 1$ produces a boundary seasonal-MA estimate, while $D = 0$ models with seasonal ARMA terms handle the observed dependence without that red flag. The airline model is a canonical benchmark for raw seasonal data such as AirPassengers, not a rule for every monthly series.

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
5                auto.arima -607.6658 -607.5083 -575.6208
4        (1,1,2)(1,0,1)[12] -607.4004 -607.2825 -579.9333
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 quantitative evidence that terms at the seasonal frequency are 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: within this candidate set, both a seasonal AR and a seasonal MA term are needed to capture the observed dependence. That statement describes the fitted dynamics; it does not identify why the seasonally adjusted source series retains them. Third, auto.arima()’s pick, ARIMA$(1,1,2)(2,0,1)_{12}$, beats candidate 4 on AIC and AICc by about a quarter of a point and loses on BIC by about 4, with one more parameter: Module 4’s AIC-vs-BIC split again. We take the BIC choice, candidate 4, for parsimony, and then check it.

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.

Diagnostics on the Winner

The IC winner is ARIMA$(1,1,2)(1,0,1)_{12}$ — the BIC winner, and the AIC winner among the hand-built candidates. 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
# (Q here is the seasonal MA order, written Q_s in the prose)
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_s = 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 tested one-term seasonal reductions fail Ljung-Box at $h = 24$. Among this candidate set of seasonal reductions, the full model is therefore the smallest specification that still passes. This demonstrates the general-to-specific pruning step; it does not claim that every possible restriction has been tested or that Mizon’s procedure has been executed to completion over the entire model space.

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.

What the Leftover Was: Residual Seasonality

Now name it. The Module 4 ARMA(1,2) on the published series showed serial correlation at the seasonal lags — a verdict on that fit. The series it was fitted to is the output of a seasonal-adjustment procedure (Section 5.8.4), and under the course’s terminology, seasonal autocorrelation that survives in a published adjusted series is residual seasonality: the adjustment procedure’s own serial correlation, showing up in the data it published. The seasonal ARMA terms of ARIMA$(1,1,2)(1,0,1)_{12}$ are what absorbed it. That is the reinterpretation of the cliffhanger: the lag-12 spike is best read not as a property of “the unemployment rate” but as the trace of one model inside the residuals of another (with the limits above).

One more thing the numbers say, worded carefully as what they show. The lag-12 residual autocorrelation of the non-seasonal fit was negative ($\hat{\rho}_{12} = -0.145$), and so were its echoes at lags 24 and 36 ($-0.125$ and $-0.130$). Under a non-seasonal ARMA, negative autocorrelation at the seasonal lags is consistent with the published series carrying too little variation at the seasonal frequencies rather than too much — a seasonal deficit, which is the signature one would expect from an adjustment that removed slightly more than the seasonal component. So the residual seasonality in this series reads as a dip, not a peak, and the seasonal ARMA terms of the winner are modeling that dip. The Deeper Dive (Section 5.9.6, not assessable) measures the shape directly in the frequency domain and finds the same thing. Two honest limits on this reading: it is consistent with an adjustment-side explanation, not a claim about how the Bureau of Labor Statistics runs its procedure; and the numbers do not establish why the dip is there — the “Fit and Compare” discussion above already listed the other explanations the residual ACF cannot rule out. You do not need the frequency domain to see the direction of the leftover; the sign of the residual ACF does that.

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. Here the model passes the manually specified Box.test() check at $h = 10$ and the monthly seasonal check at $h = 24$, both with the appropriate fitdf; the default checkresiduals() call automatically uses $h = 24$ for this series. 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. That metadata warns that imposing $D = 1$ needs fresh evidence; it does not by itself predict a boundary estimate. The fitted $\hat{\Theta}_1 \approx -1$ is the evidence here that the imposed seasonal difference is being undone and is likely one difference too many. Raw seasonal series (UNRATENSA, retail sales not seasonally adjusted, AirPassengers) are where the full airline-model machinery, $D = 1$ included, earns its keep — and Section 5.8.3 showed that even there, the airline model is a benchmark to improve on, not an answer.

5.9 Deeper Dive — A Definition You Can Test

Deeper Dive · not assessable

Everything in Section 5.9 is optional. It is never required for credit and never appears on a problem set or exam. It exists because Section 5.8.4 ended on a question — how do you diagnose residual seasonality without a definition of seasonality? — and the instructor’s research is one answer. It is distilled from the working paper The Seasons They Are A-Changin’: A Century of Definitions and a Way Forward by Carter Bryson and Gary Cornwall (U.S. Bureau of Economic Analysis), and it runs on the authors’ R package freqseas, which is public. Every number below is computed when these notes are rendered.

5.9.1 Two Lossless Representations of One Series

Every finite series $\{x_t\}_{t=1}^{T}$ has two equivalent descriptions connected by the discrete Fourier transform (DFT):

\[d_x(\omega_j) = \sum_{t=1}^{T} x_t\, e^{-i\omega_j t}, \qquad x_t = \frac{1}{T}\sum_{j=0}^{T-1} d_x(\omega_j)\, e^{i\omega_j t}, \qquad \omega_j = \frac{2\pi j}{T}.\]

The time domain orders the observations by calendar position — the language of every ACF and every lag polynomial in this course. The frequency domain reallocates the same variance by rate of repetition: how much of the series’ movement repeats every two months, every six months, every year, every decade. Nothing is gained or lost by changing coordinates; the second formula recovers the series exactly from the first. Think of the series as a mixture and the DFT as the centrifuge that separates it into its pure tones.

For monthly data, the annual cycle lives at $\omega = 2\pi/12 = \pi/6$ radians per observation, and its harmonics — the cycles that repeat an integer number of times per year — at $\pi/3$ (semiannual), $\pi/2$ (quarterly), $2\pi/3$, $5\pi/6$, and $\pi$ (every two months). A seasonal series is one whose variance piles up at those six frequencies. That is the spectral tradition from Section 5.6.4, and the reason it is attractive is that it is model-free: deterministic seasonality, stochastic seasonality, and seasonal unit roots all show up the same way, as concentrated variance at the seasonal frequencies, whatever generated them.

The sample estimate of how variance is spread across frequencies is the periodogram,

\[I_T(\omega_j) = \frac{1}{T}\left|\sum_{t=1}^{T} x_t\, e^{-i\omega_j t}\right|^2,\]

one ordinate per Fourier frequency. A small helper computes it from base R’s spec.pgram() with no tapering or smoothing (the canonical copy is in helpers/seasonality.R):

periodogram_df <- function(x) {
  # Raw periodogram of a series on the Fourier grid of the padded series
  # (spec.pgram's default fast = TRUE pads T to a highly composite length;
  # 719 becomes 720): omega in radians
  # per observation on (0, pi], I the ordinate. No taper, no smoothing, no
  # detrending; the mean is removed. Difference first if the series is I(1).
  sp <- stats::spec.pgram(as.numeric(x), taper = 0, detrend = FALSE, demean = TRUE, plot = FALSE)
  data.frame(omega = 2 * pi * sp$freq, I = sp$spec)
}

5.9.2 The Same Seasonality in Both Domains

Here is the first difference of UNRATENSA — whose ACF spiked at lags 12, 24, 36 in Section 5.8.1 — seen from the frequency side, with the six seasonal harmonics marked:

pg_nsa <- periodogram_df(d_unratensa)
harmonics <- data.frame(omega = (1:6) * pi / 6,
                        label = c("pi/6", "pi/3", "pi/2", "2pi/3", "5pi/6", "pi"))

ggplot(pg_nsa, aes(omega, I)) +
  geom_segment(aes(xend = omega, yend = 0), colour = "navy") +
  geom_vline(data = harmonics, aes(xintercept = omega), colour = "darkorange", linetype = "dashed") +
  scale_x_continuous(breaks = harmonics$omega, labels = harmonics$label) +
  ggtitle("Periodogram of diff(UNRATENSA): variance by frequency, seasonal harmonics dashed") +
  xlab("omega (radians per month)") + ylab("periodogram ordinate") + theme_bw()

The same seasonality that produced the time-domain spikes produces frequency-domain spikes, sitting on the dashed lines. Look at which line is tallest: on this series the dominant seasonal energy is at $\pi/3$, the semiannual harmonic — the twice-a-year rhythm of a January jump and a June jump — not at the annual fundamental $\pi/6$. The lag-12 spike cannot tell you that; a spike at lag 12 is a spike at lag 12 whether it comes from a once-a-year or a twice-a-year cycle. The answer is spread across the other lags; the periodogram shows it at a glance.

5.9.3 The Definition: Draw the Line, Then Compare Peaks

Section 5.6.6 said a definition is a line you draw. Here is the line, drawn before any data are examined.

Choose the cycles of interest. For monthly data with an annual period, the seasonal frequencies are the six harmonics $\Omega = \{\pi/6, \pi/3, \pi/2, 2\pi/3, 5\pi/6, \pi\}$.
Partition the frequency axis $(0, \pi]$ into $M$ equal-width bins $B_1, \ldots, B_M$. A bin is seasonal if it contains one of the harmonics; call the set of seasonal bins $\mathcal{M}_1$ and the rest $\mathcal{M}_0$. $M$ is a design resolution chosen in advance; the package picks a boundary-safe value for you.
Compare the tallest peaks. For a spectral density $f$, let $b_m(f)$ be the highest value of $f$ inside bin $m$. Define

\[\Delta(f) = \max_{m \in \mathcal{M}_1} b_m(f) \;-\; \max_{m \in \mathcal{M}_0} b_m(f).\]

Definition. A series is seasonal, relative to the chosen cycles and partition, if and only if $\Delta > 0$.

In plain language: the chosen cycles are demonstrably larger than every other cycle. The tallest peak among the seasonal bins beats the tallest peak among the non-seasonal bins. It is a relative statement — a comparison of two mountain ranges against a common sea level — so it captures a narrow spike at $\pi/6$ and a diffuse hump around $\pi/3$ alike, and it makes no commitment to what generated the peak. Whatever lands in a seasonal bin counts, trading-day leakage included; that is the “no clean isn’t” of Section 5.6.6 made operational.

One line on notation, because it matters. The course has used $\Delta$ since Module 1 for the difference operator, $\Delta y_t = y_t - y_{t-1}$. The $\Delta$ in this section is a different object: a population quantity (a difference of two peak heights), and its sample version $\hat\Delta$ below is a test statistic. The hat and the absence of an operand are the tell — $\hat\Delta$ never acts on a $y_t$. The notation dictionary records both.

5.9.4 The Test: A Null Distribution in Closed Form

Replace the unknown spectral density with the periodogram, take the same two maxima, and you have the sample statistic:

\[\hat\Delta = \max_{m \in \mathcal{M}_1} b_m(I_T) \;-\; \max_{m \in \mathcal{M}_0} b_m(I_T).\]

The reason this is a test and not just a descriptive number is a classical fact about the periodogram. Under the null that the series is white noise with innovation variance $\sigma^2$, the ordinates $I_T(\omega_j)/\sigma^2$ at distinct Fourier frequencies are approximately independent $\text{Exp}(1)$ draws. Extreme-value theory then does all the work: the maximum of $N$ independent exponentials is approximately Gumbel with location $\sigma^2 \log N$ and scale $\sigma^2$, and the difference of two independent Gumbels with a common scale is Logistic. With $N_1$ ordinates in the seasonal bins and $N_0$ in the non-seasonal bins,

\[\hat\Delta \;\big|\; H_0 \;\Rightarrow\; \text{Logistic}\!\left(\sigma^2 \log\frac{N_1}{N_0},\; \sigma^2\right).\]

(The working paper writes $\tau$ for this scale parameter; the course keeps $\tau$ for the Dickey-Fuller statistic of Module 2 and $\sigma^2$ for an innovation variance, so the null is written here in course notation.) No simulation, no bandwidth, no kernel: the critical value is a quantile of a named distribution whose only data-dependent input is the partition.

Two preparatory steps make the null honest on real data, and they are the same two steps a careful diagnostic always takes:

Prewhiten. A series with ordinary (non-seasonal) autocorrelation — a stationary AR(1), say — has a periodogram that is coloured: high at low frequencies for positive $\phi_1$. That colour inflates the non-seasonal maximum and shrinks $\hat\Delta$, so the raw statistic has almost no power against real seasonality. The fix is to fit a non-seasonal ARIMA chosen by BIC first (differencing once if the series is $I(1)$, then an AR whose order BIC picks) and run the test on its residuals. In the package this whitening step estimates the AR from the non-seasonal ordinates only, so that seasonal power cannot leak into the whitener.
Standardize the whitened residuals to unit variance, so that $\sigma^2 = 1$ and the null collapses to a single known Logistic$(\log(N_1/N_0), 1)$.

That is the whole test. Install the package once, then run it on the raw series:

install.packages("remotes")
remotes::install_github("grycrnwll/freqseas")

library(freqseas)

test_nsa <- seas_test(unratensa)   # defaults: annual period, M chosen by the package,
                                   # alpha = 0.05, difference if needed, AR order by BIC
summary(test_nsa)                  # verdict, statistic and critical value, ordinate counts,

freqseas seasonality test -- summary

Detection
  decision    : seasonal
  p-value     : 0  (alpha = 0.05)
  statistic   : 122.2  (critical = 2.191)
  ordinates   : N1 = 113 seasonal, N0 = 240 nonseasonal

Specification
  spec        : band
  shoulder p  : 0.00606
  phase R     : 0.99

Bins
  M = 19 (suggest_M, offset-free)
  rho_M = 0.263, offset_u = 0

Whitener
  d = 1 (lag-1 ACF 0.96 > 0.95 (auto: differenced once))
  AR order p = 2
  AR coefficients: 0.0922, 0.1878
  seasonal alignment: trimmed 9 leading residuals (first 12 observations unadjusted)
  BIC path (M0-Whittle):
    0 = -1002   1 = -1000   2 = -1006   3 = -1002

Donor pool size: 216 (of 708 whitened observations)

Per-harmonic table
# A tibble: 6 × 4
  omega label  excess elevated
  <dbl> <chr>   <dbl> <lgl>   
1 0.524 pi/6     8.90 TRUE    
2 1.05  pi/3   122.   TRUE    
3 1.57  pi/2    73.1  TRUE    
4 2.09  2*pi/3   5.67 TRUE    
5 2.62  5*pi/6  50.4  TRUE    
6 3.14  pi       6.79 TRUE

                                   # the whitener it built, and the per-harmonic excess table

The verdict is not close: $\hat\Delta \approx 122$ against a 5% critical value of about 2.19, and a $p$-value that is zero to machine precision. The summary also records what the test did on the way — the series was differenced once because its lag-1 autocorrelation (0.96) exceeded the package’s 0.95 threshold, an AR(2) whitener was chosen by BIC, $M = 19$ bins were used, and the seasonal bins hold $N_1 = 113$ ordinates against $N_0 = 240$ non-seasonal ones. The per-harmonic table at the bottom says where the seasonality lives, and it repeats what the raw periodogram showed: every harmonic is elevated, and the excess at $\pi/3$ (about 122) and $\pi/2$ (about 73) dwarfs the annual fundamental $\pi/6$ (about 9).

The package draws the periodogram the test actually looked at — prewhitened and standardized, with the seasonal bins shaded and the extreme-value threshold marked:

autoplot(test_nsa) +
  ggtitle("What the test sees: prewhitened, standardized periodogram of UNRATENSA, seasonal bins shaded")

5.9.5 Adjustment as the Logical Conjugate

Section 5.8.4 said adjustment is a model, and that the field’s tests and the field’s adjustment filters answer different questions — the test asks about one population object, the filter removes another, and the two do not speak the same language. The design principle of the working paper is to close that gap: the adjustment should remove exactly the excess the test detects, so that re-running the test on the adjusted series fails to reject by construction. Test and adjustment share a partition and a target. They are logical conjugates.

The operator that delivers it is stochastic spectral imputation (SSI), in three moves:

Setup. Take the DFT of the whitened series. The ordinates in the seasonal bins are the targets; the quiet ordinates below a chosen quantile of the non-seasonal floor are the donors.
Impute. For each target, draw a donor at random and reset the target’s magnitude to the donor’s level, keeping the target’s own phase. This is the stochastic step — the seasonal bins are filled with draws from the non-seasonal floor rather than set to zero, so the adjusted series keeps a realistic noise level at those frequencies instead of a spectral hole.
Finalize. Enforce conjugate symmetry, invert the DFT, recolour through the whitener, and restore the levels.

The call below is the classic SSI path — every seasonal-bin ordinate imputed, each keeping its own phase as in step 2 (phase_rule = "zero" is the package’s name for this classic rule, not an instruction to zero phases) — with a fixed seed so the draw is reproducible:

adj_ssi <- seas_adjust(unratensa, phase_rule = "zero", target_set = "bins", seed = 6376)
adj_ssi

SSI seasonal adjustment (classic)
 detection p = 0 | spec: band (shoulder p = 0.00606, phase R = 0.99)
 phase rule: zero  [declared identification]
 imputation: bootstrap, B = 1, targets = bins (113), seed = 6376
 M = 19 (suggest_M, offset-free) | whitener: d = 1, AR(2) on M0 ordinates | seasonal alignment trimmed 9 leading residuals (first 12 obs. unadjusted)
 post-adjustment detection p = 0.877

ssi_nsa  <- adjusted(adj_ssi)     # the adjusted series (first 12 months left as raw, see below)
test_ssi <- seas_test(ssi_nsa)    # re-run the same test on the output
test_ssi

freqseas seasonality test
 decision: not seasonal (detection p = 0.877, alpha = 0.05)
 spec: line (shoulder p = 0.872, phase R = 0.26)
 M = 19 (suggest_M, offset-free) | whitener: d = 1, AR(2) on M0 ordinates | seasonal alignment trimmed 9 leading residuals (first 12 obs. unadjusted)

unlist(test_ssi$evt)              # statistic, p-value, critical value, N1, N0

  statistic           p    critical          N1          N0 
 -2.7166387   0.8768991   2.1911879 113.0000000 240.0000000

The post-adjustment test fails to reject — $\hat\Delta \approx -2.7$, $p \approx 0.88$ — which is the conjugate property doing what it promises. (The adjustment leaves the first twelve months unadjusted: with one difference and an AR(2) whitener, the seasonal alignment of the imputed series starts after the whitener’s burn-in, and the package says so in its printout.)

One honest sentence about the package, because you will be tempted to call it with no options: seas_adjust()’s default branch is a different, experimental object — on a series the test classifies as a “band” (a broad seasonal hump rather than a line), the default divides by an estimated gain and imputes nothing, and on UNRATENSA its own post-test rejects. The “by construction” claim attaches to the classic call above, which is the path the working paper describes. Verify it yourself:

adj_default <- seas_adjust(unratensa, seed = 6376)   # package defaults: the band branch
signif(c(post_p_default = glance(adj_default)$post_p,
         post_p_classic = glance(adj_ssi)$post_p), 3)

post_p_default post_p_classic 
      2.71e-10       8.77e-01

Lay the SSI output against the published series and the raw series:

cmp$ssi <- as.numeric(ssi_nsa)

ggplot(cmp, aes(date)) +
  geom_line(aes(y = raw, colour = "UNRATENSA (raw)"), linewidth = 0.4, alpha = 0.5) +
  geom_line(aes(y = published, colour = "UNRATE (published)"), linewidth = 1) +
  geom_line(aes(y = ssi, colour = "classic SSI (freqseas)"), linewidth = 0.7, linetype = "42") +
  scale_colour_manual(values = c("UNRATENSA (raw)" = "grey55",
                                 "UNRATE (published)" = "darkorange",
                                 "classic SSI (freqseas)" = "darkgreen"), name = NULL) +
  ggtitle("Stochastic spectral imputation of UNRATENSA against the published series") +
  ylab("Percent") + xlab(NULL) + theme_bw() + theme(legend.position = "bottom")

keep <- seq_along(cmp$date) > 12   # exclude the unadjusted head
round(c(corr_with_published = cor(cmp$ssi[keep], cmp$published[keep]),
        mean_abs_diff       = mean(abs(cmp$ssi[keep] - cmp$published[keep]))), 3)

corr_with_published       mean_abs_diff 
              0.997               0.298

SSI tracks the published series but is rougher than X-13’s default output — a mean absolute gap of about 0.3 percentage points against X-13’s 0.065. That is by design: SSI removes seasonal excess down to the noise floor and stops; it does not try to replicate a production adjustment, and it does not smooth. The working paper argues that this is a feature — X-13 tends to remove more variance than the seasonal component contains and to dig a dip at the seasonal frequencies — but that argument is beyond these notes.

5.9.6 The Bridge Back to the Epilogue

Finally, the question Section 5.8.4 could not answer without a definition: does the published UNRATE still carry seasonality? Run the same test on it, and on the X-13 output from Section 5.8.4 for comparison:

test_pub <- seas_test(unrate)
test_pub

freqseas seasonality test
 decision: not seasonal (detection p = 0.233, alpha = 0.05)
 spec: line (shoulder p = 0.225, phase R = 0.42)
 M = 19 (suggest_M, offset-free) | whitener: d = 1, AR(3) on M0 ordinates | seasonal alignment trimmed 8 leading residuals (first 12 obs. unadjusted)

unlist(test_pub$evt)              # statistic, p-value, critical value, N1, N0

  statistic           p    critical          N1          N0 
  0.4361415   0.2333676   2.1911879 113.0000000 240.0000000

tidy(test_pub)                    # per-harmonic excess

# A tibble: 6 × 4
  omega label  excess elevated
  <dbl> <chr>   <dbl> <lgl>   
1 0.524 pi/6    -4.61 FALSE   
2 1.05  pi/3    -3.85 FALSE   
3 1.57  pi/2    -4.45 FALSE   
4 2.09  2*pi/3  -4.58 FALSE   
5 2.62  5*pi/6  -4.28 FALSE   
6 3.14  pi      -4.42 FALSE

test_x13 <- seas_test(sa_x13)
test_x13

freqseas seasonality test
 decision: not seasonal (detection p = 0.651, alpha = 0.05)
 spec: line (shoulder p = 0.641, phase R = 0.27)
 M = 19 (suggest_M, offset-free) | whitener: d = 1, AR(3) on M0 ordinates | seasonal alignment trimmed 8 leading residuals (first 12 obs. unadjusted)

unlist(test_x13$evt)

  statistic           p    critical          N1          N0 
 -1.3779603   0.6512888   2.1911879 113.0000000 240.0000000

Both fail to reject: $\hat\Delta \approx 0.44$ with $p \approx 0.23$ on the published series, $\hat\Delta \approx -1.4$ with $p \approx 0.65$ on the default X-13 output. Under this definition, neither adjusted series has seasonal excess. And the per-harmonic table for the published series already hints at the epilogue’s reading: every harmonic’s measured excess is negative — the seasonal-bin ordinates sit below the donor level.

Is it a dip? A peak-dominance test is one-sided: it was built to detect excess, and it correctly reports none. To ask the opposite question — whether the seasonal frequencies carry less power than the bands beside them — we can measure the shape directly on the same object the test looked at. Take the test’s fixed partition, and for each seasonal bin divide the mean standardized periodogram ordinate of the whitened residuals by the mean over the bin’s two neighbouring non-seasonal bins. A ratio below one is a dip; above one is an excess. The canonical copy of this helper lives in helpers/seasonality.R:

seasonal_bin_ratios <- function(test) {
  # Seasonal-bin dip check on a freqseas seas_test() object: for each seasonal
  # bin of the test's fixed partition, the mean standardized periodogram
  # ordinate of the whitened residuals divided by the mean over the bin's
  # neighbouring non-seasonal bins. Below 1 = dip, above 1 = excess. The
  # Nyquist ordinate is excluded, as in the package's detection statistic.
  breaks <- test$partition$breaks
  M1     <- test$partition$M1
  n_bins <- length(breaks) - 1L
  n      <- test$pgram$n
  j      <- seq_len(floor(n / 2))
  omega  <- 2 * pi * j / n                 # Fourier frequencies, radians per obs
  I      <- test$pgram$pgram_std[j + 1L]   # element i is Fourier index j = i - 1
  keep   <- omega > 0 & abs(omega - pi) > 1e-9
  omega  <- omega[keep]; I <- I[keep]
  bin    <- findInterval(omega, breaks, left.open = TRUE, rightmost.closed = TRUE)

  neighbours <- function(m) setdiff(intersect(c(m - 1L, m + 1L), seq_len(n_bins)), M1)
  labels <- if (length(M1) == 6) c("pi/6", "pi/3", "pi/2", "2pi/3", "5pi/6", "pi") else paste0("bin ", M1)
  harm   <- if (length(M1) == 6) (1:6) * pi / 6 else rep(NA_real_, length(M1))

  rows <- lapply(seq_along(M1), function(k) {
    m <- M1[k]; nb <- neighbours(m); nb_mean <- mean(I[bin %in% nb])
    # the single ordinate sitting on the harmonic itself (NA if none is on the
    # Fourier grid, as at the Nyquist frequency, whose ordinate is excluded)
    ih <- which.min(abs(omega - harm[k]))
    on_grid <- !is.na(harm[k]) && abs(omega[ih] - harm[k]) <= 0.51 * pi / length(omega)
    data.frame(harmonic = labels[k], bin = m,
               bin_mean = mean(I[bin == m]), neighbour_mean = nb_mean,
               ratio = mean(I[bin == m]) / nb_mean,
               harmonic_ratio = if (on_grid) I[ih] / nb_mean else NA_real_)
  })
  out <- do.call(rbind, rows)
  nb_all <- unique(unlist(lapply(M1, neighbours)))
  attr(out, "overall") <- mean(I[bin %in% M1]) / mean(I[bin %in% nb_all])
  out
}

tests3 <- list("UNRATENSA (raw)" = test_nsa,
               "UNRATE (published)" = test_pub,
               "X-13 default output" = test_x13)
dips <- lapply(tests3, seasonal_bin_ratios)

dip_table <- do.call(rbind, lapply(names(dips), function(nm) {
  r <- dips[[nm]]
  data.frame(series = nm, t(setNames(round(r$ratio, 2), r$harmonic)),
             overall = round(attr(r, "overall"), 2),
             bins_below_1 = sum(r$ratio < 1), check.names = FALSE)
}))
dip_table

               series pi/6  pi/3  pi/2 2pi/3 5pi/6   pi overall bins_below_1
1     UNRATENSA (raw) 7.01 58.78 24.93  4.69 19.13 1.78   18.39            0
2  UNRATE (published) 0.84  0.74  0.68  0.87  0.53 0.84    0.75            6
3 X-13 default output 1.41  0.86  0.55  0.88  0.69 0.74    0.82            5

# The ordinate on each harmonic itself, over the same neighbour mean
harm_table <- do.call(rbind, lapply(names(dips), function(nm) {
  r <- dips[[nm]]
  data.frame(series = nm, t(setNames(round(r$harmonic_ratio, 2), r$harmonic)), check.names = FALSE)
}))
harm_table

               series  pi/6    pi/3   pi/2 2pi/3  5pi/6 pi
1     UNRATENSA (raw) 95.84 1069.05 419.36 37.37 285.25 NA
2  UNRATE (published)  0.07    0.63   0.14  0.06   0.38 NA
3 X-13 default output  0.03    0.10   0.29  0.13   0.13 NA

# The time-domain link: lag-12 autocorrelation of the same whitened residuals
sapply(tests3, function(t) round(acf(t$e, lag.max = 12, plot = FALSE)$acf[13], 3))

    UNRATENSA (raw)  UNRATE (published) X-13 default output 
              0.837              -0.113              -0.059

The shape is unambiguous. On the published series every one of the six seasonal bins sits below its neighbours — about three-quarters of the neighbouring power overall (ratio $\approx 0.75$) — and the ordinates on the harmonics themselves are lower still (second table). The default X-13 output shows the same shape, slightly shallower (overall $\approx 0.82$, five of six bins below one). Its one exception, the $\pi/6$ bin, is instructive: the bin’s mean is pulled above one by mass beside the harmonic, while the ordinate at $\pi/6$ is only 0.03 of its neighbours. The raw series shows the opposite, an excess of about eighteen times by bin and far more at the harmonic ordinates. That is a dip at the seasonal frequencies in the published series — less seasonal-frequency power than the bands on either side — which is the shape a seasonal filter that removes more than the seasonal component would leave.

The last line closes the loop with Module 4. The lag-12 sample autocorrelation $\hat{\rho}_{12}$ of a series is a spectrum-weighted average of $\cos(12\omega)$, and $\cos(12\omega) = +1$ exactly at the seasonal harmonics, so a deficit of power there pulls $\hat{\rho}_{12}$ negative. The whitened residuals of the published series have $\hat{\rho}_{12} \approx -0.11$; those of the raw series, $+0.84$; and the Module 4 ARMA(1,2) residuals on the published series had $\hat{\rho}_{12} = -0.145$ — a non-seasonal ARMA cannot manufacture a trough at those frequencies, so the trough shows up in its residuals. The seasonal-bin dip and the negative Module 4 spikes are the same fact seen from two sides; the test’s fail-to-reject only says no excess was detected. What none of this establishes is why the dip is there; that is the epilogue’s honest limit, and the reason the seasonal ARMA terms in Module 6’s model are there.

References for Sections 5.6.4 and 5.9. Bell, W. R., & Hillmer, S. C. (1984), “Issues involved with the seasonal adjustment of economic time series,” Journal of Business & Economic Statistics. Bryson, C., & Cornwall, G., The Seasons They Are A-Changin’: A Century of Definitions and a Way Forward, working paper, U.S. Bureau of Economic Analysis. Canova, F., & Hansen, B. E. (1995), “Are seasonal patterns constant over time? A test for seasonal stability,” Journal of Business & Economic Statistics. Falkner, H. D. (1924), “The measurement of seasonal variation,” Journal of the American Statistical Association. Granger, C. W. J. (1978), “Seasonality: causation, interpretation, and implications,” in Zellner (ed.). Hillmer, S. C., & Tiao, G. C. (1982), “An ARIMA-model-based approach to seasonal adjustment,” Journal of the American Statistical Association. Hylleberg, S., Engle, R. F., Granger, C. W. J., & Yoo, B. S. (1990), “Seasonal integration and cointegration,” Journal of Econometrics. Kallek, S. (1978), “An overview of the objectives and framework of seasonal adjustment,” in Zellner (ed.). Nerlove, M. (1964), “Spectral analysis of seasonal adjustment procedures,” Econometrica. Zellner, A. (ed.) (1978), Seasonal Analysis of Economic Time Series, U.S. Department of Commerce, Bureau of the Census.

Common Pitfalls and Misconceptions

“My data has serial correlation, so I need to fix the data.” Backwards, and the vocabulary says why. The data has autocorrelation — every time series does; that is what the ACF measures. Serial correlation is a verdict on your fitted model’s residuals: a property of the 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_s$ (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 a series is already seasonally adjusted — and most headline FRED macro series are; Section 5.8.4 showed that one X-13 call on the raw unemployment rate reproduces the published UNRATE — do not assume $D = 1$. The seasonal difference the raw series needed has, in effect, already been taken by the adjustment procedure. Require evidence of a remaining seasonal unit root; otherwise seasonal differencing risks over-differencing. Treat a seasonal MA estimate stuck near $\pm 1$ as a red flag that sends you back to the differencing decision, not as an ordinary interior 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.
“The seasonally adjusted series is the data.” It is not. It is a model’s output — trend-cycle plus irregular from an X-13ARIMA-SEATS run at the agency, with the agency’s regressors, outliers, and identifying restriction baked in. What the survey measured is the not seasonally adjusted series. Know which one you are modeling, and remember that the procedure can leave residual seasonality behind, because it is a model and models can be wrong.
Reading the differencing decision from only one side. $D = 1$ on an already adjusted series drives $\hat{\Theta}_1$ toward $-1$ (one seasonal difference too many); $D = 0$ on a raw series with a seasonal unit root drives $\hat{\Phi}_1$ toward $+1$ (one too few). Both are boundary estimates, both mean “revisit $D$,” and neither shows up in an IC table, because IC cannot compare across differencing orders. Look at the coefficients, not just the $p$-values.

Connection to Enders

Residual diagnostics and Box-Jenkins model adequacy checking: Enders Chapter 2, pp. 72–79
The Ljung-Box Q statistic and fitted-model degrees-of-freedom correction: Enders Chapter 2, p. 68
Seasonal processes, multiplicative seasonality, and seasonal differencing: Enders Chapter 2, pp. 96–102
The airline model: Enders Chapter 2, p. 102
General-to-specific vs specific-to-general: Mizon (1995), “Progressive Modelling of Macroeconomic Time Series: The LSE Methodology”
Seasonal adjustment at recognition level (decomposition, X-11 and SEATS): Hyndman & Athanasopoulos, Chapter 3 (time series decomposition)
The periodogram and spectral analysis (for the Deeper Dive): Shumway & Stoffer, Chapter 4; Hamilton, Chapter 6

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 in their ARIMA chapter (fpp2 Chapter 8; fpp3 Chapter 9, which replaces checkresiduals() with ljung_box()), including residual checks and the seasonal-differencing decision; note that their $(P,D,Q)_m$ notation uses $m$ where we use $s$ for the seasonal period, and that their decomposition chapter treats seasonal adjustment as a descriptive tool where the course treats it as a model with its own diagnostics. Shumway & Stoffer measure frequency in cycles per observation on $(0, 1/2]$ where the Deeper Dive uses radians on $(0, \pi]$ — multiply by $2\pi$ to convert. 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?
Three paths on UNRATENSA. Load the cached data/UNRATENSA.csv over the 1960–2019 window, as in Section 5.8.1.
1. Path 1. Regress the first difference on eleven month dummies plus an intercept. Plot the residual ACF to lag 36 and report $\hat{\rho}_e$ at lags 12, 24, and 36 against the band; run Ljung-Box at $h = 24$ with fitdf = 0. In two sentences, say what the dummies absorbed and what they left, and why the low lags also breach the band.
2. Path 2. Fit the airline model, ARIMA$(1,1,2)(1,0,1)_{12}$, auto.arima(unratensa), and ARIMA$(1,1,2)(0,1,1)_{12}$. Also fit ARIMA$(1,1,2)$, $(1,1,2)(1,0,0)_{12}$ and $(1,1,2)(0,0,1)_{12}$, so your $(d, D) = (1, 0)$ table has the four rows of Section 5.8.3. Report the airline model’s Ljung-Box verdict at both horizons and say where its rejection comes from. Report $\hat{\Phi}_1$ from the $D = 0$ fit and say, in one sentence, what a value of 0.99 means for the differencing decision. Build the IC table for the two $(d, D) = (1, 1)$ rows and, separately, for the $D = 0$ rows, stating why they cannot share a table and where auto.arima()’s pick belongs. Run checkresiduals() on your winner.
3. Path 3. Run seasonal::seas(unratensa) once, with no options. Overlay final() of the result on the published UNRATE and report the correlation and the mean absolute difference. One sentence: which series have you been modeling since Module 2?
4. The epilogue. Fit the Module 4 non-seasonal ARMA(1,2) to diff(unrate) — the published series. In one sentence, say what its residual ACF at lag 12 shows, and which of the three terms from Section 5.2.5 is the right name for it.
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. Three terms, kept apart: autocorrelation belongs to a series (the data has it); serial correlation is the verdict on a fitted model; residual seasonality is the adjustment procedure’s serial correlation, left in a published series.
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.
Nobody fully agrees what seasonality is, but you can see it. Four traditions — calendar means, seasonal unit roots, unobserved components, spectral — each consistent, none nesting the others. A definition is a line you draw; whatever lands inside is seasonal by construction. The one safe non-example: monthly sampling alone makes nothing seasonal.
Seasonality comes in two flavors, with a boundary inside one of them. Deterministic (fixed calendar effects, seasonal dummies or Fourier) and stochastic (drifting cycles, seasonal ARMA and/or seasonal differencing), with the seasonal unit root as the persistence boundary inside the stochastic case. Macro monthly series are usually stochastic — but check whether the series was already seasonally adjusted at the source before reaching for $D = 1$.
Three paths, and when each applies. Deterministic → seasonal dummies; stochastic → SARIMA; “I need it gone, not modeled” → seasonal adjustment. On the raw unemployment rate, dummies absorbed the fixed part and left a drifting part; SARIMA’s winner was ARIMA$(1,1,2)(0,1,1)_{12}$, chosen by diagnostics and interior coefficients because IC cannot compare across differencing orders; and one X-13 call reproduced the published series.
$D$ is a commitment, and the data push back from both sides. One seasonal difference too many pins $\hat{\Theta}_1$ at $-1$; one too few pins $\hat{\Phi}_1$ at $+1$. Either boundary estimate says “revisit $D$.”
Seasonal adjustment, at recognition level. Trend-cycle + seasonal + irregular; remove the seasonal; publish the rest. Statistical agencies run X-13ARIMA-SEATS; most headline FRED macro series are its outputs, including the UNRATE you have modeled since Module 2. An adjusted series is a model’s output, not the data, and the procedure can leave residual seasonality behind.
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 two-parameter benchmark for raw monthly seasonal series — a benchmark to improve on, not an answer.
The Module 6 handoff is unchanged. ARIMA$(1,1,2)(1,0,1)_{12}$ on the published, seasonally adjusted UNRATE, diagnostic-clean at both horizons, with its seasonal ARMA terms absorbing the residual seasonality the adjustment left behind.
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: split $\rightarrow$ fit on the training window $\rightarrow$ forecast $\rightarrow$ evaluate on the test window. This module was fit + check; Modules 6 and 7 add the rest.

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.

#	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

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

Residual pattern	First candidate diagnosis — not a verdict
Mostly within the \(\pm 2/\sqrt{T}\) bands, with no systematic pattern	No visible evidence of remaining serial correlation at the inspected lags
Spike at lag 1 only	May suggest omitted short-run AR or MA structure
Geometric decay from lag 1	Suggests omitted low-order AR dynamics
Apparent cutoff after a few lags	Suggests omitted low-order MA dynamics
Spikes at lags \(s, 2s, 3s\) (say \(s = 12\))	Suggests omitted dependence at the seasonal frequency — the Module 5 transition
Large correlations that persist rather than decay	May indicate insufficient differencing or another source of non-stationarity; reassess the transformation