Module 5 Lab: Model Diagnostics and Seasonality

Econ 6376 — Applied Time Series Econometrics · AI Learning Pack

How to use this lab

Everything here runs as shipped. No API key, no internet, no blanks to fill in. The point is not to get the code working. The point is to find out whether you can look at what a model left behind and say what it missed.

The one rule from Module 1 still holds:

Write your prediction down before you run the chunk.

Not in your head. In LEARNING_LOG.md, in a comment, on paper. Being wrong on paper is the single most useful thing that can happen to you today. Most prompts below ask you to sketch a residual ACF before you see it, and a sketch you can check is worth ten you can’t.

If you are working with an AI tutor, it will ask for that prediction once, wait one turn, and run the code either way. It is making room for the “huh, I didn’t expect that” moment, because that is the one you remember in November.

Where you are on the mixing board: Module 3 lit the AR and MA blocks, and Module 4 chose their orders. One band has been silent the whole time: the seasonal band, the same AR and MA machinery running at lag 12 instead of lag 1. Today you turn it on, and after today every slider of the univariate mean equation is on. Later modules add covariates (Module 8) and give the variance its own equation (Module 14).

Two halves, one picture, or as the notes put it, the same lesson in a trench coat. First you learn to check a model. Then you meet the thing Module 4’s residuals kept pointing at.

---

1. Setup

What do you need on the table? Two packages, the course helpers, and two cached series: the unemployment rate you know, and its raw twin.

library(ggplot2)
library(forecast)

source("R/arma_simulator.R")     # the Module 3 simulator
source("R/sarima_simulator.R")   # its seasonal big sibling
source("R/diagnostics.R")        # lb_both(): Ljung-Box at h = 10 and h = 24
source("R/seasonality.R")        # seasonal_dummies(), periodogram_df()

theme_set(theme_bw(base_size = 13))

# Both caches run 1948 onward; Modules 4 and 5 use 1960-01 to 2019-12.
read_window <- function(path, column) {
  raw <- read.csv(path, stringsAsFactors = FALSE)
  raw <- raw[raw$observation_date >= "1960-01-01" &
             raw$observation_date <= "2019-12-01", ]
  ts(raw[[column]], start = c(1960, 1), frequency = 12)
}
unrate      <- read_window("data/UNRATE.csv", "UNRATE")        # published, adjusted
unratensa   <- read_window("data/UNRATENSA.csv", "UNRATENSA")  # raw, not adjusted
d_unrate    <- diff(unrate)
d_unratensa <- diff(unratensa)

# Two small plotting helpers used all the way down: a sample ACF as a data
# frame, and a correlogram with the +/- 1.96/sqrt(T) reference band.
acf_df <- function(x, lag_max = 36, label = "") {
  r <- stats::acf(as.numeric(x), lag.max = lag_max, plot = FALSE)$acf[-1]
  data.frame(lag = seq_len(lag_max), rho = as.numeric(r), fit = label)
}
acf_plot <- function(d, T_n, title = NULL, subtitle = NULL, ncol = 2) {
  band <- 1.96 / sqrt(T_n)
  d$fit <- factor(d$fit, levels = unique(d$fit))
  g <- ggplot(d, aes(lag, rho)) +
    geom_hline(yintercept = c(-band, band), linetype = "dashed", colour = "grey55") +
    geom_hline(yintercept = 0, colour = "grey40", linewidth = 0.3) +
    geom_segment(aes(xend = lag, y = 0, yend = rho), colour = "steelblue4",
                 linewidth = 1) +
    labs(title = title, subtitle = subtitle, x = "lag k", y = expression(hat(rho)[k]))
  if (length(levels(d$fit)) > 1) g <- g + facet_wrap(~ fit, ncol = ncol)
  g
}

c(levels = length(unrate), changes = length(d_unrate))

 levels changes 
    720     719 

So what: you now hold the rate you know and a raw twin you have never looked at, 720 months each, and the gap between them is where this module ends up.

---

2. Last time ended on a spike

Module 4 finished with a model it liked and a picture it didn’t. An ARMA(1,2) on the monthly change in unemployment won every information criterion. Then you plotted its residual ACF, and something poked through the band at lag 12. Reproduce it: same window, same orders, same zero-mean policy.

Predict first

Before you draw it: write one sentence saying what you would need to see in this residual ACF to sign off on the model. Which lags, how big, relative to what?

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)
T_e    <- length(e_hat)

callback_acf <- acf_df(e_hat, 36)
acf_plot(callback_acf, T_e, title = "Residual ACF: ARMA(1,2) on diff(UNRATE), 1960-2019")

round(c(lag12 = callback_acf$rho[12], lag24 = callback_acf$rho[24],
        lag36 = callback_acf$rho[36], band = 1.96 / sqrt(T_e)), 3)

 lag12  lag24  lag36   band 
-0.145 -0.125 -0.130  0.073 

The low lags are quiet; Module 4 did its job there. But at lag 12 the residual autocorrelation is $-0.145$ against a band of $\pm 0.073$ ($T = 719$), and it comes back at lag 24 ($-0.125$) and lag 36 ($-0.130$). Every twelve months, same sign each time. Module 4 called this “residual annual-lag dependence” and wisely refused to act on it. Today it gets a proper name, serial correlation, and by the end of the lab a second, more specific one.

Modify and re-run

Module 4 fit this model with include.mean = FALSE, and so does this line. Delete that argument so R estimates a mean, and re-run. Before you do: will the three spikes move in the second decimal, the third, or not at all? Write one line on why a model’s mean would, or would not, reach its residual autocorrelations.

So what: the model won the scoring contest and still left a pattern behind, and winning among candidates is not the same as being right.

---

3. What residuals should look like

Picture a defensive coordinator in the film room on Monday. The play you called was $\hat y_t$, the field gave you $y_t$, and the gap is the residual $e_t$. Gaps always exist. What you hunt for is a pattern, the same unblocked blitz every week, because a pattern means your playbook is missing a page.

That is Module 1’s notation at work. You never see the innovation $\epsilon_t$; you see every residual $e_t = y_t - \hat y_t$. If your model is right, $e_t$ should behave like $\epsilon_t$: plausibly white noise. Four checks, ranked: serial correlation (today’s subject), heteroskedasticity (eyeball it now, model it in Module 14), zero mean, approximate normality.

Predict first

Which decade of the Module 4 residuals do you bet is the loudest, and how much louder than the quietest? Name the decade and a ratio.

resid_frame <- data.frame(date = as.numeric(time(e_hat)), e = as.numeric(e_hat))
ggplot(resid_frame, aes(date, e)) +
  geom_hline(yintercept = 0, colour = "grey50") +
  geom_line(colour = "steelblue4", linewidth = 0.35) +
  labs(title = "Module 4 residuals over time", x = NULL, y = expression(e[t]))

resid_frame$decade <- paste0(floor(resid_frame$date / 10) * 10, "s")
decade_sd <- aggregate(e ~ decade, data = resid_frame, FUN = sd)
decade_sd$e <- round(decade_sd$e, 3)
names(decade_sd)[2] <- "sd_of_residual"
decade_sd

  decade sd_of_residual
1  1960s          0.168
2  1970s          0.189
3  1980s          0.195
4  1990s          0.136
5  2000s          0.139
6  2010s          0.148

round(c(mean_residual = mean(e_hat)), 4)

mean_residual 
       -9e-04 

The mean is essentially zero, but the spread is not constant: the 1980s are about 1.4 times as loud as the 1990s. Keep that picture for Module 14. For the first check, you need vocabulary.

Three words the course keeps apart

Term	Meaning	Belongs to
Autocorrelation	A property of any series, measured by its ACF	the data
Serial correlation	The verdict that a fitted model’s residuals $e_t$ show autocorrelation	the model, relative to the data
Residual seasonality	Seasonal autocorrelation left in a published seasonally adjusted series	the adjustment procedure, itself a model

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.” Say “serial correlation” out loud and somebody always hears “cereal killer.” Neither one lives in your data: one lives in your model, and the other lives in the breakfast aisle.

The textbooks and most applied work use the first two words interchangeably; you will see that in the wild. The course’s line is deliberate.

Reading residuals with the Module 3 table

In Module 3 the identification table told you what to fit. Pointed at residuals, it proposes what your fit missed:

Residual ACF	First candidate diagnosis
Inside the band, no pattern	Nothing visible left
Spike at lag 1 only	An omitted short-run AR or MA term
Geometric decay from lag 1	Omitted low-order AR dynamics
Cutoff after a few lags	Omitted low-order MA dynamics
Spikes at $s, 2s, 3s$ (12, 24, 36)	Omitted seasonal dependence
Large, persistent correlations	Not enough differencing

Look back at section 2. Which row are you in?

So what: diagnostics are just the Module 3 table, pointed at the part of the data your model gave up on.

---

4. Misspecify and watch

What does a model that is one block short actually leave behind? Build data where you know the truth, fit the wrong thing on purpose, and look. Three DGPs, $T = 500$ each: AR(2) with $\phi = (0.5, 0.3)$, MA(2) with $\theta = (0.5, 0.4)$, and ARMA(2,2) with both. Each fit reports its residual ACF and a Ljung-Box $p$-value; until section 5, read a small $p$ as “structure left over”. (Same experiment as the class widget, now in R.)

set.seed(6376); y_arma22 <- arma_simulator(n = 500, phi = c(0.5, 0.3), theta = c(0.5, 0.4))
set.seed(6377); y_ma2    <- arma_simulator(n = 500, theta = c(0.5, 0.4))
set.seed(6378); y_ar2    <- arma_simulator(n = 500, phi = c(0.5, 0.3))

# Fit one ARMA(p, q), keep its residual ACF and both Ljung-Box verdicts.
try_fit <- function(y, p, q, label) {
  fit <- tryCatch(Arima(y, order = c(p, 0, q)), error = function(e) NULL)
  if (is.null(fit)) return(NULL)
  list(fit = fit, acf = acf_df(residuals(fit), 24, label), lb = lb_both(fit))
}
lb_line <- function(res) {
  data.frame(fit = res$acf$fit[1],
             p_h10 = signif(res$lb$p[1], 3), p_h24 = signif(res$lb$p[2], 3),
             biggest_rho = round(res$acf$rho[which.max(abs(res$acf$rho))], 3),
             at_lag = which.max(abs(res$acf$rho)))
}
fmt_p <- function(p) if (p < 0.001) "below 0.001" else as.character(signif(p, 2))

Beat 1: the ARMA(2,2), fit as if it were an AR(2)

Predict first

The AR block is right; the MA block is missing. You have the section 3 table, so you can probably guess where the leftover lives. The harder call is how much is left. How big is the biggest residual bar: about 0.1, 0.3 or 0.5? And does Ljung-Box reject at both horizons, or only at one?

b1 <- try_fit(y_arma22, 2, 0, "ARMA(2,2) data, AR(2) fit")
acf_plot(b1$acf, 500, title = "Beat 1")

lb_line(b1)

                        fit    p_h10   p_h24 biggest_rho at_lag
1 ARMA(2,2) data, AR(2) fit 1.38e-09 4.9e-07       0.266      2

Low lags: 0.27 at lag 2, and Ljung-Box rejects at both horizons. That is the missing MA block’s fingerprint, smeared into your residuals.

Beat 2: add the missing block

Predict first

Now fit ARMA(2,2) to the same data. Name the size of the largest residual autocorrelation you expect, to one decimal place.

b2 <- try_fit(y_arma22, 2, 2, "ARMA(2,2) data, ARMA(2,2) fit")
acf_plot(b2$acf, 500, title = "Beat 2")

lb_line(b2)

                            fit p_h10 p_h24 biggest_rho at_lag
1 ARMA(2,2) data, ARMA(2,2) fit 0.726 0.698      -0.076     18

Same data, one more block, and nothing larger than 0.08 remains ($p$ = 0.7 at $h = 24$). The pattern was never in the data; it was in the gap between the data and your model.

Beat 3: an MA(2), faked by longer and longer ARs

Beats 3-5 go one step past the notes’ Core: practice, not assessable.

An invertible MA(2) is an AR($\infty$) (Module 3). Can a long enough AR fake it?

Predict first

Fit AR($p$) for $p = 0, 1, 2, 3$ to the MA(2) data. Sketch how the biggest residual spike changes as $p$ grows. Then commit: does AR(3) pass at $h = 24$?

b3 <- lapply(0:3, function(p) try_fit(y_ma2, p, 0, paste0("MA(2) data, AR(", p, ") fit")))
acf_plot(do.call(rbind, lapply(b3, `[[`, "acf")), 500, title = "Beat 3")

do.call(rbind, lapply(b3, lb_line))

                    fit    p_h10   p_h24 biggest_rho at_lag
1 MA(2) data, AR(0) fit 0.000000 0.00000       0.528      1
2 MA(2) data, AR(1) fit 0.000138 0.00104      -0.198      3
3 MA(2) data, AR(2) fit 0.000158 0.00149      -0.207      3
4 MA(2) data, AR(3) fit 0.077500 0.03250       0.128     11

b3_right <- try_fit(y_ma2, 0, 2, "MA(2) data, MA(2) fit")
lb_line(b3_right)

                    fit p_h10 p_h24 biggest_rho at_lag
1 MA(2) data, MA(2) fit 0.692 0.143       0.118     11

The leftover shrinks, though not in a straight line: the biggest spike is 0.53 with no AR terms, about 0.2 with one or two, 0.13 with three, and the $h = 24$ $p$-value climbs from below 0.001 to 0.032. At this draw even AR(3) still rejects, if only just. Duality works, at the price of ever more lags; the MA(2) gets there directly ($p$ = 0.14).

Beat 4: the mirror image

Predict first

Fit an MA(3) to the AR(2) data. Three MA terms against a two-term AR: pass or fail at $h = 24$? One sentence on why.

b4 <- try_fit(y_ar2, 0, 3, "AR(2) data, MA(3) fit")
acf_plot(b4$acf, 500, title = "Beat 4")

rbind(lb_line(b4), lb_line(try_fit(y_ar2, 2, 0, "AR(2) data, AR(2) fit")))

                    fit p_h10 p_h24 biggest_rho at_lag
1 AR(2) data, MA(3) fit 0.000 0.000       0.271      4
2 AR(2) data, AR(2) fit 0.561 0.643      -0.072     12

An AR(2) is an MA($\infty$), and MA(3) chops it off after three terms. The residuals keep what was chopped ($p$ below 0.001), while the AR(2) fit on the same data is clean.

Beat 5: too many knobs

Predict first

Over-fit the ARMA(2,2) data with ARMA(3,3). Will the residuals pass? And will the coefficients look like $\phi = (0.5, 0.3)$, $\theta = (0.5, 0.4)$ plus two zeros?

b5 <- try_fit(y_arma22, 3, 3, "ARMA(2,2) data, ARMA(3,3) fit")
if (is.null(b5)) {
  cat("ARMA(3,3) failed to converge at this draw; try another seed.\n")
} else {
  print(lb_line(b5))
  print(round(coef(b5$fit), 3))
}

                            fit p_h10 p_h24 biggest_rho at_lag
1 ARMA(2,2) data, ARMA(3,3) fit 0.591 0.599        0.08     11
      ar1       ar2       ar3       ma1       ma2       ma3 intercept 
    0.943     0.145    -0.223     0.071     0.066    -0.151    -0.211 

The residuals pass. Now read the coefficient line: would you defend it to anyone? Extra AR and MA roots can nearly cancel, so the estimates wander. Diagnostics cannot see over-fitting. Only parsimony, Module 4’s information criteria, can.

Modify and re-run

In the misspec-data chunk, make the MA(2) stronger: theta = c(0.8, 0.6), then re-run beat 3. Before you do: will AR(3) get closer to passing or further away? Why? (How slowly do the AR($\infty$) weights of a stronger MA die out?)

Three lessons. The tools work: you did not need the DGP to see the leftover. The fit will not volunteer the problem: every model printed coefficients, and none printed “short”. And the structure came from the model: your data never moved between beats 1 and 2. A stain that survives the wash tells you the washing machine failed; you fix the machine, not the stain. Serial correlation is a model report card, not a disease the data has.

So what: when residuals carry structure, change the model, not the data, and remember that a clean report card can still belong to a bloated model.

---

5. Ljung-Box, built once by hand

How do you turn 24 bars into one verdict? With a 95% band you expect about one bar in 24 outside by pure chance, and which one it is changes the story you tell yourself. You want one number that asks whether the first $h$ residual autocorrelations are jointly consistent with zero:

$$Q(h) = T(T+2)\sum_{k=1}^{h}\frac{\hat\rho_k^{\,2}}{T-k} \;\stackrel{a}{\sim}\; \chi^2_{\,h-\text{fitdf}} .$$

The $\hat\rho_k$ are your residuals’ autocorrelations, and fitdf counts the ARMA parameters you estimated: $p + q$, or $p + q + P + Q_s$ once seasonal terms arrive, or zero for a raw series. You spent those degrees of freedom fitting; the test charges you for them. Computing it yourself once is worth more than calling Box.test() ten times, so take beat 2’s clean ARMA(2,2), fitdf $= 4$.

Predict first

Under the null, $Q(24)$ with fitdf $= 4$ averages about its degrees of freedom, 20. Name a number for this fit’s $Q(24)$ before you compute it.

e_b2  <- residuals(b2$fit)
T_b2  <- length(e_b2)
h     <- 24
rho_k <- stats::acf(e_b2, lag.max = h, plot = FALSE)$acf[-1]

Q_hand <- T_b2 * (T_b2 + 2) * sum(rho_k^2 / (T_b2 - seq_len(h)))
p_hand <- pchisq(Q_hand, df = h - 4, lower.tail = FALSE)

Q_R <- Box.test(e_b2, lag = h, type = "Ljung-Box", fitdf = 4)
c(Q_by_hand = Q_hand, Q_Box.test = unname(Q_R$statistic),
  p_by_hand = p_hand, p_Box.test = Q_R$p.value, df = unname(Q_R$parameter))

 Q_by_hand Q_Box.test  p_by_hand p_Box.test         df 
16.2970169 16.2970169  0.6980353  0.6980353 20.0000000 

Same number both ways: $Q(24) =$ 16.3 on 20 degrees of freedom. Notice which way the good news runs. In regression you wanted to reject; here the null is “nothing left”, so a large $p$-value is good news.

Two habits, fixed before you see output. Pick $h$ in advance: $h = 10$ for non-seasonal residuals, $h = 24$ for monthly data, and on monthly data report both. And never shop across $h$ until something rejects: that is $p$-hacking with a chi-square (exploration 1 measures it).

Modify and re-run

Re-run with fitdf = 0 in Box.test() and df = h in pchisq(). Before you do: will the $p$-value go up or down? Why? What would that do to your verdict on a borderline model?

So what: Ljung-Box is just the residual ACF squared and summed; the only decisions it needs from you are $h$ and fitdf, made before you look.

---

6. The verdict on Module 4’s model

So: does the ARMA(1,2) from section 2 pass? Monthly data means both horizons, with fitdf $= p + q = 3$.

Predict first

Write “pass” or “fail” for $h = 10$ and for $h = 24$ before running. If you think they disagree, say which way and why.

lb_both(fit_l4)

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

They disagree, and the disagreement is the lesson. At $h = 10$, $Q = 3.76$ on $\chi^2_7$, $p = 0.807$: no evidence of serial correlation over ten lags. At $h = 24$, $Q = 46.1$ on $\chi^2_{21}$, $p = 0.0012$: widen the window past lag 12 and the test sees what your eye saw. A test that never looks at lag 12 cannot reject because of lag 12. And the table? No decay from lag 1, no cutoff, no lag-1 spike. Just bars at 12, 24 and 36, the same sign each time.

Spikes that recur at regular intervals. What are those?

So what: if you had only checked the short horizon, you would have signed off on this model.

---

7. Can the diagnostics catch a seasonal miss on purpose?

Would you trust a smoke detector you had never seen go off? Plant something and watch: a series whose only structure is seasonal, this month equal to 0.7 times the same month last year plus a shock,

$$y_t = \Phi_1 y_{t-12} + \epsilon_t, \qquad \Phi_1 = 0.7, \quad T = 500 .$$

Capital $\Phi$ marks a seasonal coefficient; this is SARIMA$(0,0,0)(1,0,0)_{12}$, and the loop is the equation, line for line.

set.seed(2005)
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)

Now fit what someone fresh out of Module 4 would fit: a non-seasonal ARMA(1,1).

Predict first

Where will the ARMA(1,1)’s residual ACF leave the band? And which horizons will reject: $h = 10$, $h = 24$, both, or neither?

fit_wrong <- Arima(y_sar, order = c(1, 0, 1))
round(coef(fit_wrong), 3)

      ar1       ma1 intercept 
   -0.249     0.372     0.124 

acf_plot(acf_df(residuals(fit_wrong), 36), 500, title = "Residual ACF, ARMA(1,1) fit")

lb_both(fit_wrong)

   h fitdf          Q         p
1 10     2   2.741567 0.9495146
2 24     2 377.587352 0.0000000

The coefficients look unremarkable, and at $h = 10$ Ljung-Box passes ($p$ = 0.95). Only $h = 24$, which can see lag 12, rejects ($Q$ = 378). Section 6 again, on data where you know the answer. Now give the model the right class:

fit_right <- Arima(y_sar, order = c(0, 0, 0),
                   seasonal = list(order = c(1, 0, 0), period = 12))
round(coef(fit_right), 3)

     sar1 intercept 
    0.703     0.100 

lb_both(fit_right)

   h fitdf         Q         p
1 10     1  4.513781 0.8744702
2 24     1 14.714333 0.9045789

$\hat\Phi_1 =$ 0.703, close to the 0.7 you planted, and both horizons pass. Same data; the serial correlation lived in the model.

How good is this detector?

This subsection goes one step past the notes’ Core: practice, not assessable.

How weak can the seasonal structure get before Ljung-Box stops noticing? Below, sarima_simulator() draws data with $\Phi_1 \in \{0, 0.15, 0.3, 0.5\}$, fits the wrong ARMA(1,1), and counts rejections at $h = 24$. A few seconds.

Predict first

At $\Phi_1 = 0$ the ARMA(1,1) is actually fine. What rejection rate should you see there? And at what $\Phi_1$ does the test catch the miss half the time? Name a number.

set.seed(505)
Phi_grid <- c(0, 0.15, 0.3, 0.5)
n_reps   <- 60

power <- do.call(rbind, lapply(Phi_grid, function(P) {
  rejects <- replicate(n_reps, {
    y   <- sarima_simulator(n = 500, Phi = P, s = 12)
    fit <- tryCatch(Arima(y, order = c(1, 0, 1)), error = function(e) NULL)
    if (is.null(fit)) NA else
      Box.test(residuals(fit), lag = 24, type = "Ljung-Box", fitdf = 2)$p.value < 0.05
  })
  data.frame(Phi_1 = P, replications = n_reps, failed_fits = sum(is.na(rejects)),
             rejection_rate = round(mean(rejects, na.rm = TRUE), 3))
}))
power

  Phi_1 replications failed_fits rejection_rate
1  0.00           60           0          0.017
2  0.15           60           0          0.350
3  0.30           60           0          0.983
4  0.50           60           1          1.000

The $\Phi_1 = 0$ row is the test’s size: with nothing to find it should reject about 5% of the time (here 0.017, from only 60 draws). By $\Phi_1 = 0.3$ it catches the miss 98% of the time, so half-power sits somewhere between 0.15 and 0.3.

Modify and re-run

Add $0.2$ and $0.25$ to Phi_grid and raise n_reps to 150 (slower). Before you do: will the size row move toward 0.05 or away? Why do more replications help the size row more than the $\Phi_1 = 0.5$ row?

So what: the diagnostic catches a moderately strong seasonal miss reliably, but only at a horizon long enough to reach the seasonal lag.

---

8. You can see it

Here is the same unemployment rate before anyone touched it. FRED publishes a second series, UNRATENSA: the civilian unemployment rate not seasonally adjusted, what the survey measures each month. You have never looked at it.

Predict first

The Module 4 residuals had a lag-12 autocorrelation of $-0.145$. Name a number for the lag-12 autocorrelation of the raw series’ monthly changes: sign and size.

raw_frame <- data.frame(date = as.numeric(time(unratensa)), rate = as.numeric(unratensa))
ggplot(raw_frame, aes(date, rate)) +
  geom_line(colour = "grey25", linewidth = 0.35) +
  labs(title = "UNRATENSA: the unemployment rate, not seasonally adjusted, 1960-2019",
       x = NULL, y = "Percent")

raw_acf <- acf_df(d_unratensa, 36)
acf_plot(raw_acf, length(d_unratensa), title = "ACF of diff(UNRATENSA), no model fitted")

round(c(lag12 = raw_acf$rho[12], lag24 = raw_acf$rho[24], lag36 = raw_acf$rho[36]), 3)

lag12 lag24 lag36 
0.804 0.777 0.742 

# The same lag on the published, adjusted series
round(acf_df(d_unrate, 12)$rho[12], 3)

[1] -0.089

You can see the annual hiring and layoff cycle without any statistics. In the monthly changes, lags 12, 24 and 36 sit at 0.804, 0.777 and 0.742: ten times the band, barely decaying. On the published series the same lag is $-0.089$. No model has been fitted, so this is autocorrelation, a property of the raw series; section 2’s spike was serial correlation, a verdict on a fit. The pictures rhyme and mean different things.

Is every monthly series seasonal?

Predict first

Take a plain random walk and label it with frequency = 12 so R thinks it is monthly. Sketch the ACF of its monthly changes at lags 11, 12 and 13. Does lag 12 stand out?

set.seed(6376)
rw_monthly <- ts(cumsum(rnorm(720)), start = c(1960, 1), frequency = 12)
rw_acf <- acf_df(diff(rw_monthly), 36)
acf_plot(rw_acf, 719, title = "ACF of the monthly changes of a random walk labeled monthly")

round(rw_acf$rho[c(11, 12, 13)], 3)

[1]  0.027 -0.035 -0.011

Nothing: lag 12 is -0.035, inside the band. Being monthly gives a series a seasonal frequency where structure could live; it does not put structure there. That is the one safe non-example of seasonality.

Nobody agrees what it is, but you can see it

The field has argued for a century, in four traditions: calendar means, seasonal unit roots, unobserved components, spectral peaks. None nests the others. A 1978 Census–NBER volume compared the search to hunting a black cat in a dark room, and added that it could be worse: the cat might not be there. The course’s position: we do not fully agree on what seasonality is, but you can see it. You just did.

Seasonality is deterministic (a fixed calendar pattern) or stochastic (a pattern that drifts), and at the stochastic edge sits the seasonal unit root, the pattern that never mean-reverts, just as the random walk was the AR(1) at $\phi = 1$. There is no clean list of what seasonality isn’t: a definition is a line you draw, and whatever lands inside is seasonal by construction.

The triage. Fixed pattern: path 1, seasonal dummies. Drifting pattern you want to model and forecast: path 2, SARIMA. “I need it gone, not modeled”: path 3, seasonal adjustment. A journalist comparing this month with last, and a forecaster building next year’s fan chart: which path does each need?

So what: the raw series is loudly seasonal before any model touches it, and the next three sections try each path on it.

---

9. Path 1: seasonal dummies

What is the cheapest possible theory of the January jump? That January is simply January, the same bump every year: eleven month dummies plus an intercept, fit to the monthly change.

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

month_means <- coef(fit_dum)[1] + c(0, coef(fit_dum)[-1])
names(month_means) <- month.abb
round(month_means, 2)

  Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec 
 0.90 -0.07 -0.26 -0.52 -0.11  0.64 -0.16 -0.26 -0.18 -0.14  0.11  0.02 

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 = "steelblue4") + geom_hline(yintercept = 0) +
  labs(title = "Path 1: mean monthly change in the raw rate, by month, 1960-2019",
       x = NULL, y = "Percentage points")

The dummies soak up a lot: $R^2 = 0.721$ on a differenced series, with January averaging $+0.90$ points, June $+0.64$ and April $-0.52$.

Predict first

With 72% of the variance explained, will the residual ACF at lag 12 clear the band? If not, guess its size.

e_dum   <- residuals(fit_dum)
dum_acf <- acf_df(e_dum, 36)
acf_plot(dum_acf, length(e_dum), title = "Residual ACF after the month dummies")

round(c(lag12 = dum_acf$rho[12], lag24 = dum_acf$rho[24], lag36 = dum_acf$rho[36]), 3)

lag12 lag24 lag36 
0.350 0.298 0.223 

which(abs(dum_acf$rho) > 1.96 / sqrt(length(e_dum)))    # lags outside the band

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

Box.test(e_dum, lag = 24, type = "Ljung-Box", fitdf = 0)  # no ARMA terms: fitdf = 0

    Box-Ljung test

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

It does not clear. Lags 12, 24 and 36 sit at 0.350, 0.298 and 0.223, three to five times the band, and Ljung-Box (fitdf $= 0$: no ARMA terms) rejects flat out. Some of that is non-seasonal (lags 2, 3, 5, 7 and 11 also breach, because dummies carry none of Module 4’s short-run dynamics), so read the seasonal lags themselves. The fixed part is gone; a changing part remains.

Drifting, or just a noisy fixed pattern? Let the month effects differ between 1960–89 and 1990–2019 and see what you get.

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

anova(lm(d_unratensa ~ month + half), lm(d_unratensa ~ month * half))

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

round(tapply(d_unratensa, list(month, half), mean)[c("Jun", "Jul"), ], 2)

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

The June jump halves, from $+0.84$ to $+0.43$ points, and July flips sign, from $-0.33$ to $+0.02$. The $F(11, 695) = 12.5$ is nominal: it assumes uncorrelated errors, which the residual ACF above just ruled out, so its $p$-value is far too small. The size of those shifts is the evidence. The seasonal pattern of 2019 is not the pattern of 1965, and fixed dummies cannot say that; neither can the sine–cosine pairs of forecast::fourier(), because the problem is fixedness, not parsimony.

So what: dummies absorb the fixed part of a seasonal pattern and leave the part that changed over the sample in your residuals, where the diagnostics find it.

---

10. Path 2: SARIMA

Think about the tide under the weather. Day to day you notice the weather, the non-seasonal band. Underneath runs a slower rhythm that repeats every year. Model the weather, ignore the tide, and you keep getting blindsided by something perfectly predictable. SARIMA puts a second row of faders on your console:

$$\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 .$$

Lowercase is the band you know, orders $(p, d, q)$. Uppercase is the seasonal band at lags $s, 2s, \ldots$, orders $(P, D, Q)$, $s = 12$ for monthly data, same sign conventions. Standing alone in a sentence, the seasonal MA order is $Q_s$, so it cannot be confused with the Ljung-Box $Q$.

What happens when the bands multiply?

The notes work this expansion in a Deeper Dive (section 5.7); here it takes four lines of R.

The bands multiply. Take ARIMA$(1,0,0)(1,0,0)_{12}$, whose AR side is $(1 - \phi_1 L)(1 - \Phi_1 L^{12})$.

Predict first

With $\phi_1 = 0.5$ and $\Phi_1 = 0.7$, which lags get a nonzero coefficient after you multiply out? List them with their values.

polymul <- function(a, b) {                        # coefficient vectors, constant first
  out <- rep(0, length(a) + length(b) - 1)
  for (i in seq_along(a)) out[i:(i + length(b) - 1)] <- out[i:(i + length(b) - 1)] + a[i] * b
  out
}
nonseasonal <- c(1, -0.5)                          # 1 - 0.5 L
seasonal    <- c(1, rep(0, 11), -0.7)              # 1 - 0.7 L^12
product     <- polymul(nonseasonal, seasonal)
lags <- seq_along(product) - 1
data.frame(lag = lags[product != 0 & lags > 0], coefficient = product[product != 0 & lags > 0])

  lag coefficient
1   1       -0.50
2  12       -0.70
3  13        0.35

Lags 1, 12 and 13. The lag-13 term, $\phi_1\Phi_1 = 0.35$, is not a free parameter; the product forces it. Two parameters buy three lags of dynamics, which is why seasonal models stay small (and exactly what sarima_simulator() does before its loop).

Decaying or flat?

Compare a stationary seasonal AR ($\Phi_1 = 0.8$) with a seasonal random walk, $y_t = y_{t-12} + \epsilon_t$.

Predict first

Sketch both ACFs at lags 12, 24, 36. Which decays and which stays nearly flat? For the one that decays, guess the lag-24 value.

set.seed(6376)
y_sar08 <- sarima_simulator(n = 480, Phi = 0.8, s = 12)

eps_rw <- rnorm(480 + 12); y_srw <- numeric(480 + 12)
for (t in 13:(480 + 12)) y_srw[t] <- y_srw[t - 12] + eps_rw[t]
y_srw <- ts(y_srw[-(1:12)], frequency = 12)

gallery <- rbind(acf_df(y_sar08, 40, "Seasonal AR, Phi_1 = 0.8"),
                 acf_df(y_srw,   40, "Seasonal random walk, Phi_1 = 1"))
acf_plot(gallery, 480, title = "Two seasonal fingerprints")

gallery[gallery$lag %in% c(12, 24, 36), ]

   lag       rho                             fit
12  12 0.7308801        Seasonal AR, Phi_1 = 0.8
24  24 0.5138877        Seasonal AR, Phi_1 = 0.8
36  36 0.3388403        Seasonal AR, Phi_1 = 0.8
52  12 0.9418471 Seasonal random walk, Phi_1 = 1
64  24 0.8874295 Seasonal random walk, Phi_1 = 1
76  36 0.8356964 Seasonal random walk, Phi_1 = 1

The seasonal random walk is Sheldon from Module 1, remembering every January exactly, forever. Decaying spikes mean a stationary seasonal ARMA; spikes that stay high mean a seasonal unit root, fixed by the seasonal difference $\Delta_{12} y_t = y_t - y_{t-12}$. The raw unemployment spikes in section 8 went 0.804, 0.777, 0.742. Which picture is that?

Four candidates on the raw series

#	Model	Why it is here
1	$(0,1,1)(0,1,1)_{12}$	The airline model, Box and Jenkins’ two-parameter benchmark
2	$(1,1,2)(1,0,1)_{12}$	Module 4’s block plus a seasonal ARMA band, no seasonal difference
3	auto.arima()	Let the algorithm choose everything
4	$(1,1,2)(0,1,1)_{12}$	The airline’s seasonal band on Module 4’s block

Predict first

Which of the four fail Ljung-Box, at either horizon? And for candidate 2, with no seasonal difference, guess $\hat\Phi_1$ to two decimals.

seasonal_fit <- function(y, order, seasonal) {
  Arima(y, order = order, seasonal = list(order = seasonal, period = 12))
}
cand_1 <- seasonal_fit(unratensa, c(0, 1, 1), c(0, 1, 1))
cand_2 <- seasonal_fit(unratensa, c(1, 1, 2), c(1, 0, 1))
cand_3 <- auto.arima(unratensa)                    # about five seconds
cand_4 <- seasonal_fit(unratensa, c(1, 1, 2), c(0, 1, 1))

orders_of <- function(fit) {                       # fit$arma is p, q, P, Q, s, d, D
  a <- fit$arma
  sprintf("(%d,%d,%d)(%d,%d,%d)[%d]", a[1], a[6], a[2], a[3], a[7], a[4], a[5])
}
candidates <- list(`1` = cand_1, `2` = cand_2, `3` = cand_3, `4` = cand_4)
do.call(rbind, lapply(names(candidates), function(k) {
  fit <- candidates[[k]]; lb <- lb_both(fit)
  data.frame(candidate = k, model = orders_of(fit),
             fitdf = lb$fitdf[1], p_h10 = signif(lb$p[1], 3), p_h24 = signif(lb$p[2], 3))
}))

  candidate              model fitdf    p_h10    p_h24
1         1 (0,1,1)(0,1,1)[12]     2 1.54e-13 2.48e-11
2         2 (1,1,2)(1,0,1)[12]     5 9.67e-01 7.61e-01
3         3 (2,0,2)(2,1,2)[12]     8 2.66e-01 4.65e-01
4         4 (1,1,2)(0,1,1)[12]     4 8.81e-01 4.24e-01

round(coef(cand_1), 3)

   ma1   sma1 
 0.066 -0.777 

round(acf_df(residuals(cand_1), 5)$rho[2:5], 3)    # candidate 1's residuals, lags 2-5

[1] 0.214 0.139 0.119 0.112

round(coef(cand_2), 4)

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

round(coef(cand_4), 3)

   ar1    ma1    ma2   sma1 
 0.783 -0.752  0.158 -0.753 

• Candidate 1, the airline model, fails both horizons. Its seasonal half is fine ($\hat\Theta_1 = -0.78$, interior), but its non-seasonal half is too thin: $\hat\theta_1 \approx 0.07$, and Module 4’s short-run dynamics sit in its residuals at lags 2 to 5 (0.21, 0.14, 0.12, 0.11). A benchmark, not an answer.
• Candidate 2 passes, and that is the trap. $\hat\Phi_1 = 0.994$: asked to model the seasonal band without differencing it, the estimator pushed the seasonal AR coefficient as close to 1 as it could. That is the signature of a seasonal unit root, the data asking for $D = 1$ the only way a $D = 0$ model can.
• Candidate 3, auto.arima()’s pick, has eight coefficients, $D = 1$, $d = 0$, clean diagnostics. It is a very good intern; read its work before you sign it.
• Candidate 4 passes ($p = 0.88$ and $0.42$) with every coefficient interior ($\hat\Theta_1 = -0.75$).

Information criteria compare only within a differencing class

Can AIC referee? Only partly: a $D = 1$ likelihood is computed on the seasonally differenced series, a shorter, different series than a $D = 0$ model sees, so rows compare only when they share $(d, D)$.

ic_row <- function(name, fit) data.frame(model = name, d_D = paste0("(", fit$arma[6], ",", fit$arma[7], ")"),
                                         nobs = fit$nobs, AIC = round(fit$aic, 1), BIC = round(fit$bic, 1))
ic_all <- rbind(ic_row("1: airline", cand_1), ic_row("4: (1,1,2)(0,1,1)", cand_4),
                ic_row("2: (1,1,2)(1,0,1)", cand_2), ic_row("3: auto.arima", cand_3))
split(ic_all, ic_all$d_D)       # one table per differencing class (d, D)

$`(0,1)`
          model   d_D nobs  AIC    BIC
4 3: auto.arima (0,1)  708 -218 -176.9

$`(1,0)`
              model   d_D nobs    AIC    BIC
3 2: (1,1,2)(1,0,1) (1,0)  719 -189.8 -162.3

$`(1,1)`
              model   d_D nobs    AIC    BIC
1        1: airline (1,1)  707 -155.7 -142.0
2 4: (1,1,2)(0,1,1) (1,1)  707 -206.2 -183.4

In the $(1, 1)$ class, candidate 4 beats the airline model by about 50 AIC points ($-206.2$ against $-155.7$): the whole gap is the non-seasonal block. Candidates 2 and 3 sit alone in their classes, so no table picks among 2, 3 and 4. The coefficients do: candidate 2 carries the signature of a seasonal unit root in its AR polynomial, and candidate 4 does the job with four coefficients instead of eight.

What if the series was already adjusted?

The published UNRATE has already been through seasonal adjustment. What happens if you give it the airline model’s seasonal difference anyway?

Predict first

Fit the airline model to the published UNRATE, which has already been through seasonal adjustment. Guess $\hat\Theta_1$: near $-0.8$, like the raw series, or somewhere else?

airline_published <- seasonal_fit(unrate, c(0, 1, 1), c(0, 1, 1))
coef(airline_published)

        ma1        sma1 
 0.06768543 -0.99999280 

$\hat\Theta_1 = -0.99999$, pinned against the invertibility boundary: the Module 3 over-differencing signature, a model trying to undo a difference the data did not need. $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 (candidate 2) pushes $\hat\Phi_1$ to within 0.01 of $+1$. Either way, revisit $D$, and read the coefficients, since no IC table spans differencing classes. Module 4’s guardrail still stands: a residual spike at lag 12 alone never licenses $D = 1$.

So what: SARIMA models the change that dummies could not, and boundary estimates, not the IC table, tell you when the differencing is wrong.

---

11. Path 3: seasonal adjustment

A journalist wants to know whether unemployment rose last month. The raw series jumps every January and June whatever the economy does, so someone has to take the seasonal pattern out first. That is seasonal adjustment, and you need it at recognition level: what it does, and how to recognize its outputs.

It treats the series as trend-cycle + seasonal + irregular. X-13ARIMA-SEATS, the Census Bureau’s program, fits a SARIMA with calendar and outlier terms, pulls out the seasonal piece, and hands back the rest. The Bureau of Labor Statistics uses it differently from what you are about to do: it never adjusts the rate itself. It adjusts unemployment and the labor force, group by group (age and sex), with the program’s X-11 filters rather than the SEATS default you will run, and then computes the rate from the adjusted totals.

Predict first

You are about to run one X-13 call on the raw series with every setting at its default, and lay the result over the published UNRATE. How close will it get? Name the mean absolute gap in percentage points, and say whether you expect the two lines to look different on a plot.

The call needs the optional seasonal and x13binary packages; if you lack them, the chunk says so and the comparison still runs from a saved copy of the same call’s output.

x13_ready <- requireNamespace("seasonal", quietly = TRUE) &&
  requireNamespace("x13binary", quietly = TRUE) &&
  !nzchar(Sys.getenv("ECON6376_SKIP_X13"))

fit_x13 <- NULL
if (x13_ready) {
  fit_x13 <- tryCatch(seasonal::seas(unratensa), error = function(e) NULL)
}

if (is.null(fit_x13)) {
  cat("X-13 is not available in this R session, so the live call was skipped.\n",
      "To run it yourself: install.packages(c(\"seasonal\", \"x13binary\"))\n",
      "The comparison below uses data/x13_vs_published.csv, the saved output of\n",
      "this same default call.\n", sep = "")
  saved <- read.csv("data/x13_stats.csv", stringsAsFactors = FALSE)
  cat("The saved run's automatic model choice:", saved$x13_arima_model, "\n")
} else {
  cat("X-13's automatic model choice:", fit_x13$model$arima$model, "\n")
}

X-13's automatic model choice: (1 1 2)(0 1 1) 

Now the reveal, read from the saved file so everyone sees the same thing; if your live call ran, the last line checks it against the file.

x13 <- read.csv("data/x13_vs_published.csv", stringsAsFactors = FALSE)
x13$date <- as.Date(x13$date)

gap <- x13$X13_default - x13$UNRATE_published
round(c(correlation        = cor(x13$X13_default, x13$UNRATE_published),
        mean_abs_gap_pp    = mean(abs(gap)),
        max_abs_gap_pp     = max(abs(gap)),
        share_within_0.1pp = mean(abs(gap) <= 0.1)), 4)

       correlation    mean_abs_gap_pp     max_abs_gap_pp share_within_0.1pp 
            0.9985             0.0647             0.3726             0.8028 

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

[1] "1976-02"

ggplot(x13, aes(date)) +
  geom_line(aes(y = UNRATENSA, colour = "raw (UNRATENSA)"), linewidth = 0.35, alpha = 0.6) +
  geom_line(aes(y = UNRATE_published, colour = "published (UNRATE)"), linewidth = 0.9) +
  geom_line(aes(y = X13_default, colour = "one default X-13 call"), linewidth = 0.55,
            linetype = "22") +
  scale_colour_manual(values = c("raw (UNRATENSA)" = "grey60",
                                 "published (UNRATE)" = "darkorange",
                                 "one default X-13 call" = "navy"), name = NULL) +
  labs(title = "One X-13 call on the raw series, against the published series",
       x = NULL, y = "Percent") +
  theme(legend.position = "bottom")

if (!is.null(fit_x13)) {
  live_gap <- as.numeric(seasonal::final(fit_x13)) - x13$X13_default
  cat("Largest gap between your live call and the saved file:",
      signif(max(abs(live_gap)), 3), "percentage points\n")
}

Largest gap between your live call and the saved file: 0 percentage points

On this plot the two adjusted lines are one line. Correlation 0.9985, a mean absolute gap of about 0.065 points, four months in five within 0.1, and the biggest gap 0.37 points, in February 1976. X-13’s automatic choice was $(1,1,2)(0,1,1)$: your candidate 4, reached with no help from you. Close, not equal: one default call closely approximates the published series, which BLS builds from adjusted pieces, with a different filter, its own production setup and revisions (and, for the older years, earlier programs). So why does one call on the rate itself get this close? Ask your tutor to help you think it through, then check your answer against the gaps: where is the largest one, and what was the economy doing then?

Now say it out loud. The UNRATE you have modeled since Module 2 is the output of a seasonal-adjustment procedure like this one. Every ACF, unit-root test and ARMA fit in Modules 2 through 4 ran on trend-cycle plus irregular, not on what the survey measured. That is why the series is usable, and it changes what you were modeling. The cost: adjustment is a model, so it can leave serial correlation of its own behind. At the seasonal lags of a published series, that is called residual seasonality.

So what: “seasonally adjusted” means “passed through a model”, and a model’s leftovers can end up in your data.

---

12. Epilogue: residual seasonality

Why would a Wall Street Journal headline in 2015 read “The Phrase of the Day is ‘Residual Seasonality’”? The Federal Reserve Board, the Cleveland Fed and the UK’s Office for National Statistics have all published on it too. You can now see why a model’s leftovers make the news.

Back to your published series. Start general, as Mizon’s general-to-specific approach says: Module 4’s block plus a full seasonal ARMA band, and no seasonal difference, since this series has had its seasonal pattern removed once already. That is ARIMA$(1,1,2)(1,0,1)_{12}$, with fitdf $= p + q + P + Q_s = 5$.

Predict first

Will the lag-12 spike survive this model? Guess the new residual autocorrelation at lag 12, and whether Ljung-Box passes at $h = 24$.

fit_full <- seasonal_fit(unrate, c(1, 1, 2), c(1, 0, 1))
round(coef(fit_full), 3)

   ar1    ma1    ma2   sar1   sma1 
 0.907 -0.935  0.213  0.525 -0.791 

lb_both(fit_full)

   h fitdf         Q         p
1 10     5  4.000874 0.5492902
2 24     5 21.175527 0.3271836

after_acf <- acf_df(residuals(fit_full), 36)
round(c(lag12 = after_acf$rho[12], lag24 = after_acf$rho[24], lag36 = after_acf$rho[36]), 3)

 lag12  lag24  lag36 
 0.028 -0.024 -0.065 

before_after <- rbind(transform(callback_acf, fit = "Before: ARMA(1,2) on diff(UNRATE)"),
                      transform(after_acf,    fit = "After: ARIMA(1,1,2)(1,0,1)[12]"))
acf_plot(before_after, T_e, title = "Same series, two models")

Gone. The residual autocorrelations at 12, 24 and 36 are now 0.028, $-0.024$ and $-0.065$, all inside the band, and Ljung-Box passes at both horizons ($p = 0.55$ and $0.33$). The whole module sits between your two pictures: the left one asked the question, the right one answers it. But general-to-specific does not stop at “it passes”. Can you make it smaller?

Predict first

Drop the seasonal MA, giving $(1,1,2)(1,0,0)_{12}$, or drop the seasonal AR, giving $(1,1,2)(0,0,1)_{12}$. Does either trim still pass at $h = 24$?

trim_sar <- seasonal_fit(unrate, c(1, 1, 2), c(1, 0, 0))
trim_sma <- seasonal_fit(unrate, c(1, 1, 2), c(0, 0, 1))
rbind(cbind(model = "(1,1,2)(1,0,0)[12]", lb_both(trim_sar)),
      cbind(model = "(1,1,2)(0,0,1)[12]", lb_both(trim_sma)))

               model  h fitdf         Q           p
1 (1,1,2)(1,0,0)[12] 10     4  3.578364 0.733516157
2 (1,1,2)(1,0,0)[12] 24     4 43.168710 0.001941911
3 (1,1,2)(0,0,1)[12] 10     4  3.876668 0.693362172
4 (1,1,2)(0,0,1)[12] 24     4 39.710456 0.005433612

Neither survives: both fail at $h = 24$ ($p = 0.0019$ and $0.0054$). Among these reductions the full seasonal band is the smallest model that passes: one honest step of Mizon’s procedure, not a search of every possible model.

What the leftover was. The Module 4 fit showed serial correlation at the seasonal lags. But you fitted it to the output of an adjustment procedure, and seasonal autocorrelation surviving in a published adjusted series is residual seasonality: the procedure’s own serial correlation, showing up in the data it published. The seasonal ARMA terms absorbed it. The sign says something too: the spikes were negative, consistent with the published series carrying slightly too little variation at the seasonal frequencies. A dip, not a peak, as if the adjustment removed a bit more than the seasonal component.

Hold that reading loosely

The dip reading is consistent with an adjustment-side explanation. It is not a claim about how BLS runs its procedure, and the numbers do not say why the dip is there.

Modify and re-run

Fit candidate 4’s orders to the published series: seasonal_fit(unrate, c(1, 1, 2), c(0, 1, 1)). Before you run it: will $\hat\Theta_1$ land in the interior, as on the raw series, or near $-1$, like the airline model on the published series? Why?

So what: the lag-12 spike is best read not as a fact about unemployment alone but as the trace of one model inside the residuals of another.

---

13. Explorations

Optional, and outside the assessable surface

Nothing in this section is required, and none of it appears on problem sets or exams. Assessments draw only on the Core tier of the module notes. These are here because the interesting questions usually start where the syllabus stops.

If one of these grabs you, chase it. If a different question grabs you, chase that one instead — that is a better outcome than any of the six below. Your AI tutor’s job in this section is to help you turn a hunch into an experiment you can actually run, then connect what you find back to a Core concept.

1. Can you p-hack a Ljung-Box test? (guided)

How bad is shopping across $h$, really? Take residuals from a correct model, pure white noise with fitdf $= 0$. Test each series once at $h = 10$, then at every $h$ from 1 to 36, keeping the smallest $p$.

set.seed(6376)
shop <- replicate(300, {
  e <- rnorm(500)
  p <- sapply(1:36, function(h) Box.test(e, lag = h, type = "Ljung-Box")$p.value)
  c(fixed_h10 = p[10] < 0.05, best_of_36 = min(p) < 0.05)
})
rowMeans(shop)

 fixed_h10 best_of_36 
0.07666667 0.25666667 

At a fixed $h = 10$ this run rejects 7.7% of 300 honest series. Keep the best of 36 and it rejects 26%. Why isn’t that 36 times 5%? (Are $Q(10)$ and $Q(11)$ independent?)

2. Why does the airline model pass on AirPassengers and fail on UNRATENSA? (guided)

Box and Jenkins built the airline model on monthly airline passengers, 1949–1960, where it famously passes. On your raw unemployment rate it failed at the low lags.

fit_ap <- seasonal_fit(log(AirPassengers), c(0, 1, 1), c(0, 1, 1))
round(coef(fit_ap), 3)

   ma1   sma1 
-0.402 -0.557 

lb_both(fit_ap)

   h fitdf         Q         p
1 10     2  8.809729 0.3586005
2 24     2 26.445847 0.2330325

low_lags <- round(rbind(AirPassengers = acf_df(residuals(fit_ap), 5)$rho[2:5],
                        UNRATENSA     = acf_df(residuals(cand_1), 5)$rho[2:5]), 3)
colnames(low_lags) <- paste0("lag", 2:5)
low_lags

               lag2   lag3   lag4  lag5
AirPassengers 0.024 -0.125 -0.111 0.061
UNRATENSA     0.214  0.139  0.119 0.112

What is different about the two series’ short-run dynamics, and what kind of series would you guess the airline model fits well?

3. Where does the seasonality live, by rate of repetition? (guided; Deeper Dive · not assessable)

A monthly pattern can repeat once a year, twice a year, up to six times: the harmonics $\omega = k\pi/6$. The periodogram spreads variance across those rates.

pg <- periodogram_df(d_unratensa)
ggplot(pg, aes(omega, I)) +
  geom_vline(xintercept = (1:6) * pi / 6, colour = "firebrick", linetype = "dotted") +
  geom_segment(aes(xend = omega, y = 0, yend = I), colour = "steelblue4") +
  labs(title = "Periodogram of diff(UNRATENSA); dotted lines at the seasonal harmonics",
       x = expression(omega ~ "(radians per month)"), y = "ordinate")

Is the tallest harmonic the once-a-year one at $\pi/6$, or something else? What does that say about the shape of the pattern? Try d_unrate too. If you want the formal version, section 5.9 of the notes runs a frequency-domain seasonality test with the freqseas package; the starter below is off by default.

# Optional: remotes::install_github("grycrnwll/freqseas")
freqseas::seas_test(unratensa)
freqseas::seas_test(unrate)

4. Can you invent a DGP that fools the diagnostics? (open, off-syllabus)

Build a series with real structure, fit a sensible model, and get Ljung-Box to pass at both $h = 10$ and $h = 24$. Dependence that only appears at lag 30? Or volatility clustering, where shocks are uncorrelated but not independent (Module 1’s subtlety, and a preview of ARCH in Module 14)? What would you have to look at to catch your own trick?

5. What does a seasonal series you care about look like? (open, off-syllabus)

Retail sales, housing starts, electricity output, the price of eggs. Find the not seasonally adjusted version on FRED, run the triage, and try the three paths. If X-13 is installed, compare seas() on the raw series to the agency’s adjusted version: how close does one default call get this time?

Use data/UNRATENSA.csv as a template for any CSV you download by hand from FRED; no key needed. With a free FRED API key, R/fetch_fred.R pulls live, reading the key from your environment so it never lands in a file you share.

# Optional. Requires the fredr package and a FRED_KEY in your environment.
source("R/fetch_fred.R")
retail_nsa <- fetch_fred("RSXFSN", start_date = "1992-01-01", frequency = 12)

6. Is the seasonal pattern still changing? (open)

Section 9 caught June and July drifting between halves. Slow drift, or a break? Estimate the month effects in rolling 10-year windows and plot June’s path. Or extend the window past 2019 as a labeled sensitivity analysis: does 2020 break candidate 4, or the epilogue model?

---

14. Checks

R/checks.R is the referee: seven assertions about things this module claims are true, in base R so you can read every line. It does not know what you concluded and it cannot grade you. It just refuses to let an incorrect claim about the code stand, whether that claim came from you, your AI tutor, or a typo. If a check fails, CHECKS.md says what it asserts and what a failure usually means.

source("R/checks.R")

Module 5 checks
---------------
PASS Module 4 residuals spike at 12/24/36; Ljung-Box passes h = 10, fails h = 24  [rho = -0.1449 / -0.1252 / -0.1299; p(h=10) = 0.807, p(h=24) = 0.0012]
PASS Hand-built Q(24) equals Box.test(), on 24 - fitdf degrees of freedom  [hand = 46.115372, Box.test = 46.115372, df = 21]
PASS Seasonal DGP: ARMA(1,1) rejected at h = 24; SARIMA recovers Phi_1 and passes  [wrong-class p = 0; Phi_1-hat = 0.703, right-class p = 0.905]
PASS Dummies: 11 columns, R^2 near 0.721, residual lag-12 near 0.350, June falls and July rises between halves  [cols = 11, R^2 = 0.7213, rho_12 = 0.350, June shift = -0.41, July shift = +0.35]
PASS D too few pins Phi_1 near +1; D too many pins Theta_1 near -1; the winner is interior and clean  [Phi_1 (D=0, raw) = 0.9939; Theta_1 (airline, published) = -0.999993; winner Theta_1 = -0.753, p = 0.881 / 0.424]
PASS ARIMA(1,1,2)(1,0,1)_12 on published UNRATE passes both horizons, lag 12 inside the band  [fitdf = 5, p = 0.549 / 0.327, rho_12 = 0.028]
PASS UNRATE/UNRATENSA caches and the saved X-13 comparison are intact  [window = 720 + 720 months; X-13 file rows = 720, summary matches = TRUE]
---------------
7/7 checks passed.

---

Where you are now

You turned on the seasonal band, the last sliders of the univariate mean equation, and found:

• A model can win every information criterion and still leave structure in its residuals, and residuals that don’t behave like innovations mean the model is wrong.
• Autocorrelation belongs to a series; serial correlation is a verdict on a fit; residual seasonality is an adjustment procedure’s serial correlation.
• Ljung-Box needs $h$ and fitdf chosen first; on monthly data, check $h = 10$ and $h = 24$. It catches under-fitting and misses over-fitting.
• Seasonality can be fixed or drifting, three paths treat it, and $D$ is a commitment the coefficients will argue with.
• The UNRATE you have modeled since Module 2 is a model’s output.

You now hold a model that passes every check you know how to run. The open question you carry into Module 6 is the one a model is finally for: what does it say about next month, and how far should you trust it?

Before you close this file: open LEARNING_LOG.md and write down, in your own words, one prediction you got wrong today and what you now think instead. Your words, not your tutor’s — that is the entire value of the exercise.