Module 9 Lab: Granger Causality in Dynamic Regression

Econ 6376 — Applied Time Series Econometrics · AI Learning Pack

How to use this lab

Everything runs as shipped, offline. Your work is to predict, run, modify and explain. Record predictions before the corresponding chunks, in your private LEARNING_LOG.md, an R comment, or on paper. With an AI tutor, the prediction ask takes one turn; “just run it” releases it. Printed subprompts are optional self-study opportunities: they do not re-arm the tutor’s prediction ask within one work unit, such as a whole lab section and its variations. No AI is needed to use the lab.

We follow the deck’s sequence: the predictive question → the joint test → simulations with a known mechanism → claims and unemployment → uses and limits. The notes are canonical for facts, notation and the Core boundary. Sections 2–7 are the main path; sections 8–9 are optional.

Where are we on the mixing board? No new slider. We keep Module 8’s ARDL board — outcome lags on, observed-predictor lags on, explicit MA terms and deterministic trend off — and ask whether switching off the lagged- $x$ block costs us predictive information. The current- $x$ term $\beta_0x_t$ is left out of both regressions in the past-only comparison.

1. Setup

Open this file from the pack folder and run chunks in order. In RStudio, choose Session → Set Working Directory → To Source File Location before running individual chunks; getwd() should show ai_packs/module_09 and dir.exists("data") should be TRUE. Quarto’s file render uses the lab’s directory automatically. Installed R packages are ggplot2, urca, knitr and rmarkdown; rendering also needs Quarto. Render this file by path, not the whole folder.

The canonical course helpers granger_nested_test() and granger_bg_test() live in R/granger_utils.R. We build the F statistic by hand first, then load the helpers in section 4 and make them agree with our arithmetic.

2. The predictive question

Keep three questions apart:

Question	What an answer would tell us
Do $x$ and $y$ move together?	Association.
Does the past of $x$ improve prediction of $y$ beyond the history already used?	Additional predictive information — the Granger question.
Would intervening on $x$ change $y$ ?	A causal intervention effect, under an identification argument.

The past of $x$ Granger-causes $y$ if it contains information useful for predicting $y$ beyond the past of $y$ and the other information specified in the comparison. Three words do the work: incremental (with versus without $x$ ’s history), directional ( $x$ helping predict $y$ is a different question from $y$ helping predict $x$ ), and conditional (say what else is known). In shorthand the null is $H_0: x\not\rightarrow y$ and the alternative $H_1: x\rightarrow y$ ; the arrows stand for a joint restriction on a whole included lag block, conditional on the stated own lags and other information.

Hamilton’s thesis is the module’s thesis: these are tests of whether one series helps forecast another, not of whether it causes it; they bear on the causal question only together with further assumptions (Hamilton §11.2; see the notes for his oil-price example). Claims can be an early signal of a labor-demand deterioration without being a lever: changing the number in the claims database would not change employment.

Within-period effects are a different question

Add a same-period term to the working ARDL:

$y_t=\alpha+\sum_{j=1}^{p}\phi_j\,y_{t-j} +\underbrace{\beta_0x_t}_{\text{same period}} +\underbrace{\sum_{\ell=1}^{q}\beta_\ell\,x_{t-\ell}}_{\text{past }x\text{: the Granger question}}+\epsilon_t.$

The Granger null concerns only $\beta_1,\ldots,\beta_q$ . A same-period link is what Granger (1969) called instantaneous causality. Our past-only comparison leaves $x_t$ out of both regressions; it does not test $\beta_0$ .

Take the simplest case: $x_t$ independent over time and $y_t=0.8x_t+\epsilon_t$ . Changing $x_t$ changes $y_t$ in this DGP.

Predict: in a past-only comparison ( $y_{t-1}$ versus $y_{t-1}$ plus $x_{t-1}$ ), will the test (a) reject decisively because $x$ clearly affects $y$ , (b) typically not reject because yesterday’s $x$ carries nothing about today’s $y$ , or (c) fail to run?

set.seed(6376093)
n_wp <- 600
x_all <- rnorm(n_wp + 1L)
y_all <- 0.8 * x_all + rnorm(n_wp + 1L)
wp <- data.frame(y = y_all[-1], x = x_all[-1],
                 y_l1 = y_all[-(n_wp + 1L)], x_l1 = x_all[-(n_wp + 1L)])
wp_r <- lm(y ~ y_l1, data = wp, na.action = na.fail)
wp_u <- lm(y ~ y_l1 + x_l1, data = wp, na.action = na.fail)
wp_past_p <- anova(wp_r, wp_u)[2, "Pr(>F)"]
wp_beta0 <- unname(coef(lm(y ~ y_l1 + x + x_l1, data = wp))["x"])
# Residuals from the two past-only equations: same-period dependence.
wp_rx <- lm(x ~ x_l1 + y_l1, data = wp, na.action = na.fail)
wp_cor <- cor(residuals(wp_u), residuals(wp_rx))
knitr::kable(data.frame(
  Quantity = c("Past-only Granger p-value (x lag in y equation)",
               "Estimated beta_0 when current x is included",
               "Correlation of the two equations' residuals"),
  Value = c(format.pval(wp_past_p, digits = 3), round(wp_beta0, 3),
            round(wp_cor, 3))))

Quantity	Value
Past-only Granger p-value (x lag in y equation)	0.18
Estimated beta_0 when current x is included	0.842
Correlation of the two equations’ residuals	0.647

The past-only p-value is 0.18: lagged $x$ carries no information here, so the test has nothing to find. Yet $\hat\beta_0=$ 0.842 when current $x$ enters, close to the DGP value 0.8, and the two equations’ residuals correlate at 0.647. That residual correlation is symmetric — it is the same number whichever variable you call the outcome — so by itself it cannot say which moved first within the period.

Modify: put current x in both regressions and re-test x_l1. That asks whether earlier $x$ adds information given current $x$ — a different question again. Putting x inside the tested block instead would make a joint test of instantaneous and lagged terms. Neither is our past-only comparison.

Expectations can reverse the apparent story

In Hamilton’s present-value example (§11.2), stock prices are set by expected future dividends and do not affect dividends. Yet past prices help predict dividends while past dividends do not help predict prices: the price already contains investors’ news. The predictive direction runs opposite to the structural one. Forward-looking behavior is a central reason to keep “helps forecast” and “causes” apart.

3. Turn the question into a joint test

The unrestricted and restricted regressions are

$y_t=\alpha+\sum_{j=1}^{p}\phi_jy_{t-j}+\sum_{\ell=1}^{q}\beta_\ell x_{t-\ell}+\epsilon_t \quad\text{and}\quad y_t=\alpha+\sum_{j=1}^{p}\phi_jy_{t-j}+\epsilon_t,$

with $H_0:\beta_1=\cdots=\beta_q=0$ against “at least one $\beta_\ell\ne0$ .” Here $q$ is the largest predictor lag, not an MA order. Re-estimate everything in the restricted model; do not freeze the unrestricted estimates. The null is the whole block: one insignificant $\hat\beta_2$ is not this test, and neither is $\beta_1+\beta_2+\beta_3=0$ , because offsetting lag effects can still change predictions.

With $RSS_R$ and $RSS_U$ from the two fits on the same $n$ observations, and $k_U$ counting every unrestricted coefficient including the intercept,

$F=\frac{(RSS_R-RSS_U)/q}{RSS_U/(n-k_U)}.$

The reference $F(q,n-k_U)$ is an approximation under maintained dynamic- regression conditions, not an exact finite-sample law.

Predict: $n=120$ ; the unrestricted model has an intercept, two own lags and three $x$ lags; $RSS_R=108$ , $RSS_U=90$ . Count $k_U$ and $q$ first. Is $F$ closer to (a) 0.2, (b) 7.6, or (c) 18?

F_example <- ((108 - 90) / 3) / (90 / (120 - 6))
data.frame(F = F_example, numerator_df = 3, denominator_df = 114,
           p_value = pf(F_example, 3, 114, lower.tail = FALSE))

    F numerator_df denominator_df      p_value
1 7.6            3            114 0.0001118079

Now write the statistic as a function of two fitted models. It is deliberately naive: it trusts its inputs.

f_by_hand <- function(restricted, unrestricted) {
  rss_r <- sum(residuals(restricted)^2)
  rss_u <- sum(residuals(unrestricted)^2)
  n <- nobs(unrestricted)
  k_u <- length(coef(unrestricted))
  q <- k_u - length(coef(restricted))
  F_stat <- ((rss_r - rss_u) / q) / (rss_u / (n - k_u))
  data.frame(n_R = nobs(restricted), n_U = n, q = q, k_U = k_u,
             RSS_R = rss_r, RSS_U = rss_u, F = F_stat,
             p_value = pf(F_stat, q, n - k_u, lower.tail = FALSE))
}

On a common sample the unrestricted model can reproduce the restricted fit by setting its extra coefficients to zero, so $RSS_U\le RSS_R$ . A negative RSS reduction means the models or rows are wrong, not that you found something.

Reverse the direction with another equation. To ask whether $y$ helps predict $x$ , make $x_t$ the outcome with its own lags and test the added $y$ -lag block. That is a new regression and a new null.

An in-sample test of a predictive idea

The F-test judges a predictive idea by in-sample fit. Ashley, Granger and Schmalensee (1980) argued that credible Granger testing should lean on out-of-sample forecasting. The trade-off: in-sample fit can reward overfitting, while a pre-chosen holdout spends data and power. Ashley and Tsang (2014) cross-sample validation splits at $\tau$ , fits each side and predicts the other, so every prediction uses coefficients from data not used to make it. Cornwall, Mills, Sauley and Weng (2019) hold out one block at a time and compare the restricted and unrestricted models’ held-out losses.

Predict: in-sample, $RSS_U\le RSS_R$ always on a common sample. On held-out data, can the unrestricted model predict worse? Hold your answer for section 6.

4. Simulate a known one-way relationship

$x_t=0.55x_{t-1}+\epsilon_{x,t},\qquad y_t=0.4y_{t-1}+0.6x_{t-1}+\epsilon_t,$

with mutually independent standard normal innovations. The $y$ equation has lagged $x$ ; the $x$ equation has no lagged $y$ . The for-loop makes the timing literal.

Predict: which direction has population predictive content? Must every finite sample produce exactly one rejection and one nonrejection — (a) yes, (b) no, sampling variation can produce either error?

set.seed(6376091)
n <- 600
burn <- 200
total <- n + burn + 1L
epsilon_x <- rnorm(total)
epsilon <- rnorm(total)
x <- y <- numeric(total)
for (t in 2:total) {
  x[t] <- 0.55 * x[t - 1] + epsilon_x[t]
  y[t] <- 0.4 * y[t - 1] + 0.6 * x[t - 1] + epsilon[t]
}
keep <- (burn + 2L):total
sim <- data.frame(t = seq_len(n), y = y[keep], x = x[keep],
                  y_l1 = y[keep - 1L], x_l1 = x[keep - 1L])
sim_y_r <- lm(y ~ y_l1, data = sim, na.action = na.fail)
sim_y_u <- lm(y ~ y_l1 + x_l1, data = sim, na.action = na.fail)
sim_x_r <- lm(x ~ x_l1, data = sim, na.action = na.fail)
sim_x_u <- lm(x ~ x_l1 + y_l1, data = sim, na.action = na.fail)
knitr::kable(rbind(
  cbind(Direction = "x history -> y", f_by_hand(sim_y_r, sim_y_u)),
  cbind(Direction = "y history -> x", f_by_hand(sim_x_r, sim_x_u))),
  digits = 4)

Direction	n_R	n_U	q	k_U	RSS_R	RSS_U	F	p_value
x history -> y	600	600	1	3	853.8729	595.2698	259.3547	0.0000
y history -> x	600	600	1	3	571.4650	570.4480	1.0643	0.3026

Lag columns come from the pre-sample burn-in, so both regressions use the same 600 rows. Now confirm the arithmetic against R’s own nested-model F test and the course helper:

source("R/granger_utils.R")
anova(sim_y_r, sim_y_u, test = "F")

Analysis of Variance Table

Model 1: y ~ y_l1
Model 2: y ~ y_l1 + x_l1
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    598 853.87                                  
2    597 595.27  1     258.6 259.35 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sim_helper <- granger_nested_test(sim_y_r, sim_y_u)
c(by_hand = f_by_hand(sim_y_r, sim_y_u)$F,
  anova = anova(sim_y_r, sim_y_u)[2, "F"], helper = sim_helper$F)

 by_hand    anova   helper 
259.3547 259.3547 259.3547

The cross-lag estimate is $\hat\beta_1=$ 0.600 against the DGP value 0.6 — agreement after rounding, not exact recovery. The reverse-direction p-value is 0.303. A p-value is evidence from one sample, not the definition of the mechanism.

ggplot(data.frame(t = sim$t, e = residuals(sim_y_u)), aes(t, e)) +
  geom_line(colour = "#315375") +
  labs(x = "Observation", y = "Fitted residual",
       title = "Known one-way DGP: residuals from the unrestricted outcome equation")

Line plot of fitted residuals from the known one-way simulation, indexed by observation.

Break the comparison: unequal rows

Build lag columns by padding inside the retained sample instead, and let lm() drop incomplete rows silently. The restricted model uses one own lag; the unrestricted model adds three $x$ lags.

Predict: will f_by_hand() (a) refuse, (b) quietly return a number, or (c) return exactly the same F as a correct common-sample comparison?

trap <- data.frame(y = sim$y, x = sim$x)
for (j in 1:3) {
  trap[[paste0("y_l", j)]] <- c(rep(NA_real_, j), head(trap$y, -j))
  trap[[paste0("x_l", j)]] <- c(rep(NA_real_, j), head(trap$x, -j))
}
bad_r <- lm(y ~ y_l1, data = trap)                       # drops 1 row
bad_u <- lm(y ~ y_l1 + x_l1 + x_l2 + x_l3, data = trap)  # drops 3 rows
naive <- f_by_hand(bad_r, bad_u)
helper_says <- tryCatch(granger_nested_test(bad_r, bad_u)$F,
                        error = function(e) paste("refused:", conditionMessage(e)))
common <- trap[complete.cases(trap), ]
good_r <- lm(y ~ y_l1, data = common, na.action = na.fail)
good_u <- lm(y ~ y_l1 + x_l1 + x_l2 + x_l3, data = common, na.action = na.fail)
knitr::kable(rbind(cbind(Comparison = "Unequal rows (naive)", naive),
                   cbind(Comparison = "Common rows", f_by_hand(good_r, good_u))),
             digits = 4)

Comparison	n_R	n_U	q	k_U	RSS_R	RSS_U	F	p_value
Unequal rows (naive)	599	597	3	5	846.0794	584.137	88.4895	0
Common rows	597	597	3	5	845.9761	584.137	88.4546	0

helper_says

[1] "refused: identical(rownames(r), rownames(u)) is not TRUE"

The naive function returned a number from 599 restricted rows and 597 unrestricted rows: the RSS difference mixes a model change with a data change. Here only 2 rows differ, so the naive F (88.49) happens to sit near the common-rows F (88.45) — closeness in one draw is luck, not validity. With twelve own lags in one model and three in the other, or a missing month in the middle, the gap need not be small. The helper checks matching outcomes, rows and nested design columns and refuses. The repair is to build the complete lagged sample before fitting either model and use na.action = na.fail so nothing is dropped silently.

Modify and reconcile: keep the seed and innovations, replace 0.6 by zero in the loop, regenerate $y$ and re-test both directions. The null is now true in both. A p-value below 0.05 would be a false rejection, not a hidden effect. Then shorten the retained sample and explain why weaker evidence does not mean a weaker DGP coefficient.

5. A useful signal need not be a lever

Unobserved news $z_t$ arrives independently each period. The recorded signal is noisy news, and $y$ responds to yesterday’s news:

$x_t=z_t+\eta_t,\qquad y_t=0.4y_{t-1}+0.8z_{t-1}+\epsilon_t.$

$z_t$ , $\eta_t$ and $\epsilon_t$ are mutually independent standard normal draws. There is no $x$ term in the structural equation for $y$ .

Predict: (1) Does $x_{t-1}$ help predict $y_t$ when $z$ is omitted — (a) yes, (b) no? (2) Once $z_{t-1}$ is included, is the sample test (a) guaranteed not to reject, or (b) able to reject even though the population cross-lag coefficient is zero?

set.seed(6376092)
z <- rnorm(total)
eta <- rnorm(total)
epsilon <- rnorm(total)
x <- z + eta
y <- numeric(total)
for (t in 2:total) {
  y[t] <- 0.4 * y[t - 1] + 0.8 * z[t - 1] + epsilon[t]
}
driver <- data.frame(y = y[keep], y_l1 = y[keep - 1L],
                     x_l1 = x[keep - 1L], z_l1 = z[keep - 1L])
pair_r <- lm(y ~ y_l1, data = driver, na.action = na.fail)
pair_u <- lm(y ~ y_l1 + x_l1, data = driver, na.action = na.fail)
cond_r <- lm(y ~ y_l1 + z_l1, data = driver, na.action = na.fail)
cond_u <- lm(y ~ y_l1 + z_l1 + x_l1, data = driver, na.action = na.fail)
driver_results <- rbind(
  cbind(Information = "Own history only", granger_nested_test(pair_r, pair_u)),
  cbind(Information = "Own history plus z", granger_nested_test(cond_r, cond_u)))
knitr::kable(transform(driver_results[, c("Information", "n", "q", "F", "p_value")],
  p_value = format.pval(p_value, digits = 3, eps = 1e-6)), digits = 3)

Information	n	q	F	p_value
Own history only	600	1	100.286	<1e-06
Own history plus z	600	1	5.875	0.0157

The first comparison finds genuine predictive information: the signal reveals part of yesterday’s news. It is not a false positive merely because $x$ does not drive $y$ . The second asks whether the signal helps once the actual news is known; its population coefficient is zero. In this draw its p-value is 0.01566 — a false rejection in this sample. We keep this draw on purpose. Adding the driver removes $x$ ’s population content; it does not guarantee a sample nonrejection.

Now impose a change on recorded $x$ , holding $z$ , the initial outcome and its innovations fixed. The $y$ equation contains no $x$ , so the $y$ path does not move. Informative in observational data; not a driver under intervention.

How often does a true null get rejected?

One false rejection is one draw. Repeat the conditional comparison, whose null is true by construction.

Predict: at the 5% level, across 200 independent samples, will the rejection rate be closer to (a) 0%, (b) 5%, or (c) 25%?

set.seed(6376095)
B <- 200
null_p <- numeric(B)
for (b in seq_len(B)) {
  zb <- rnorm(total); etab <- rnorm(total); epsb <- rnorm(total)
  xb <- zb + etab
  yb <- numeric(total)
  for (t in 2:total) yb[t] <- 0.4 * yb[t - 1] + 0.8 * zb[t - 1] + epsb[t]
  db <- data.frame(y = yb[keep], y_l1 = yb[keep - 1L],
                   x_l1 = xb[keep - 1L], z_l1 = zb[keep - 1L])
  null_p[b] <- granger_nested_test(lm(y ~ y_l1 + z_l1, data = db),
                                   lm(y ~ y_l1 + z_l1 + x_l1, data = db))$p_value
}
null_rate <- mean(null_p < 0.05)
c(rejection_rate = null_rate, MCSE = sqrt(null_rate * (1 - null_rate) / B))

rejection_rate           MCSE 
    0.04500000     0.01465862

4.5% of these null-true samples reject at 5%. That is the size working as designed, within Monte Carlo error. Our kept draw is one of the roughly one-in-twenty.

6. Application: claims and unemployment

$U_t$ is the monthly unemployment rate (percent). $C_t$ is the mean of the weekly initial-claims observations whose week-ending Saturdays fall in month $t$ — not a monthly total. Both are seasonally adjusted. As in Module 8,

$y_t=\Delta U_t,\qquad x_t=100\Delta\log C_t.$

$y$ is in percentage points; $x$ is approximately percentage growth of the monthly claims mean.

Dates on the data and dates on the clock

The forecast origin waits until both previous-month measures are released, before the target month’s unemployment release: usually a nowcast of the current calendar month, written with information $\mathcal I_t^{\mathrm{release}}$ . The calendar is messier than the index:

Frequency and aggregation. $U_t$ refers to the CPS reference week (the week containing the 12th) and is released early in month $t+1$ . $C_t$ averages weeks before and after that reference week, so the monthly “lag” $x_{t-1}$ is not a clean calendar lead: the late weeks inside $C_{t-1}$ fall within the period $y_t$ measures. Aggregation can create or hide predictive direction, and can push a one-week response into the same month.
Seasonal adjustment. BLS and the Department of Labor adjust the series separately; adjusted history is revised with later data, so a past-only comparison may not be strictly past-only.
Release timing. Weekly claims come out each Thursday. By the time $U_t$ is released, claims covering month $t$ ’s reference week have been public for about two weeks. Our past-only design leaves them out on purpose, so it understates claims’ practical nowcasting value. Using them would be the same-period question of section 2.

The caches are revised values collected on different dates: this is a revised-data pseudo-out-of-sample exercise, not a real-time vintage study.

Load and check the data

Predict: if one week’s claims were missing, should the monthly mean be (a) computed from the remaining weeks, or (b) treated as an error?

u <- read.csv(file.path(data_dir, "UNRATE.csv"))
c_week <- read.csv(file.path(data_dir, "ICSA.csv"))
u$observation_date <- as.Date(u$observation_date)
c_week$observation_date <- as.Date(c_week$observation_date)
u <- u[u$observation_date >= as.Date("1988-12-01") &
         u$observation_date <= as.Date("2019-12-01"), ]
c_week <- c_week[order(c_week$observation_date), ]
stopifnot(!anyNA(u), !anyNA(c_week), all(c_week$ICSA > 0),
  identical(u$observation_date,
            seq(as.Date("1988-12-01"), as.Date("2019-12-01"), by = "month")),
  identical(c_week$observation_date,
            seq(as.Date("1989-01-07"), as.Date("2019-12-28"), by = "week")),
  all(as.POSIXlt(c_week$observation_date)$wday == 6))
u$y <- c(NA_real_, diff(u$UNRATE))
u$U_previous <- c(NA_real_, head(u$UNRATE, -1L))
u <- u[u$observation_date >= as.Date("1989-01-01"), ]
c_week$month <- as.Date(format(c_week$observation_date, "%Y-%m-01"))
stopifnot(all(table(c_week$month) %in% c(4L, 5L)))
c_month <- aggregate(ICSA ~ month, data = c_week, FUN = mean)
app <- merge(u, c_month, by.x = "observation_date", by.y = "month", all = TRUE)
names(app)[match(c("observation_date", "UNRATE", "ICSA"), names(app))] <-
  c("date", "U", "C")
stopifnot(nrow(app) == 372L, !anyNA(app))
app$x <- c(NA_real_, 100 * diff(log(app$C)))
for (j in 1:6) {
  app[[paste0("y_l", j)]] <- c(rep(NA_real_, j), head(app$y, -j))
  app[[paste0("x_l", j)]] <- c(rep(NA_real_, j), head(app$x, -j))
}
app <- app[app$date >= as.Date("1990-01-01"), ]
stopifnot(nrow(app) == 360L, !anyNA(app))
train <- app[app$date <= as.Date("2017-12-01"), ]
test <- app[app$date >= as.Date("2018-01-01"), ]
stopifnot(nrow(train) == 336L, nrow(test) == 24L)

Training is January 1990–December 2017 (336 rows); January 2018–December 2019 (24 rows) is reserved for evaluation. 1989 supplies differences and lags, so every specification below uses the same training outcomes. The window avoids the pandemic, which also limits what the results say about later periods.

paths <- rbind(
  data.frame(date = train$date, value = train$y,
             Series = "Unemployment change (percentage points)"),
  data.frame(date = train$date, value = train$x,
             Series = "Claims: 100 x log change"))
ggplot(paths, aes(date, value)) + geom_line(colour = "#315375") +
  facet_wrap(~ Series, ncol = 1, scales = "free_y") +
  labs(x = NULL, y = NULL, title = "Transformed training data, 1990-2017")

Two training-sample time-series panels: unemployment changes in percentage points and 100 times log changes in the monthly mean of weekly claims, from 1990 to 2017.

The transformations are a stationary working specification. Reuse Module 2: ADF has a unit-root null; KPSS has a level-stationarity null. The ADF lag choice below does not choose the ARDL’s lag orders.

stationarity <- data.frame()
for (variable in c("y", "x")) {
  adf <- urca::ur.df(train[[variable]], type = "drift", lags = 12, selectlags = "BIC")
  kpss <- urca::ur.kpss(train[[variable]], type = "mu", lags = "short")
  stationarity <- rbind(stationarity, data.frame(
    Variable = variable, ADF = as.numeric(adf@teststat[1, "tau2"]),
    ADF_5pct = as.numeric(adf@cval["tau2", "5pct"]),
    KPSS = as.numeric(kpss@teststat), KPSS_5pct = as.numeric(kpss@cval[1, "5pct"])))
}
knitr::kable(stationarity, digits = 4)

Variable	ADF	ADF_5pct	KPSS	KPSS_5pct
y	-3.6919	-2.87	0.2535	0.463
x	-12.7947	-2.87	0.1263	0.463

Both ADF statistics fall below their 5% critical values and neither KPSS statistic exceeds its own. That supports the working specification; it does not prove a time-invariant DGP.

Fix the comparison, then diagnose it

Main specification in each direction: six own lags and three lags of the other variable. For unemployment this carries forward the specification reached in Module 8’s exploratory training work; it is not a fresh prespecified confirmatory test. The reverse equation reuses the lag lengths for comparability.

Forward null: conditional on six lags of unemployment changes, the three claims-growth lag coefficients are all zero. The reverse null exchanges roles.

Predict: for the forward equation, count $n$ , $q$ , $k_U$ and the residual degrees of freedom before running. Then: if the nominal F rejects but the Breusch–Godfrey diagnostic also rejects, do you (a) report a Granger-causality verdict, (b) report the nominal result and the limitation, or (c) add lags until the diagnostic passes?

y_r <- lm(reformulate(paste0("y_l", 1:6), response = "y"),
          data = train, na.action = na.fail)
y_u <- lm(reformulate(c(paste0("y_l", 1:6), paste0("x_l", 1:3)), response = "y"),
          data = train, na.action = na.fail)
x_r <- lm(reformulate(paste0("x_l", 1:6), response = "x"),
          data = train, na.action = na.fail)
x_u <- lm(reformulate(c(paste0("x_l", 1:6), paste0("y_l", 1:3)), response = "x"),
          data = train, na.action = na.fail)
y_test <- granger_nested_test(y_r, y_u)
x_test <- granger_nested_test(x_r, x_u)
y_bg <- granger_bg_test(y_u, order = 12)
x_bg <- granger_bg_test(x_u, order = 12)
main_results <- data.frame(
  Direction = c("Claims growth -> unemployment change", "Unemployment change -> claims growth"),
  n = c(y_test$n, x_test$n), q = c(y_test$q, x_test$q), k_U = c(y_test$k_U, x_test$k_U),
  F = c(y_test$F, x_test$F), Nominal_F_p = c(y_test$p_value, x_test$p_value),
  BG_12_p = c(y_bg$p_value, x_bg$p_value))
knitr::kable(transform(main_results,
  Nominal_F_p = format.pval(Nominal_F_p, digits = 3, eps = 1e-6),
  BG_12_p = format.pval(BG_12_p, digits = 3, eps = 1e-6)), digits = 3)

Direction	n	q	k_U	F	Nominal_F_p	BG_12_p
Claims growth -> unemployment change	336	3	10	21.058	<1e-06	0.0031
Unemployment change -> claims growth	336	3	10	1.804	0.146	0.2549

The Breusch–Godfrey (BG) diagnostic regresses the fitted residual on the original regressors and its own lags 1–12 (dropping the first twelve rows), and compares $n_aR_a^2$ with an asymptotic $\chi^2(12)$ . Its null is no residual serial correlation through lag 12. It is a specification warning, not an exogeneity test.

Forward: nominal F p-value 1.71e-12, BG(12) p-value 0.0031. Reverse: nominal p-value 0.146, BG(12) p-value 0.255. The forward equation fails the residual-dynamics check at 5%. Its tiny nominal p-value is therefore not a clean basis for a validated Granger-causality verdict. Report the arithmetic and the limitation, and stop the illustration here. The reverse nonrejection is “not detected with these lags,” not “unemployment changes carry no information about claims.”

residual_paths <- rbind(
  data.frame(date = train$date, e = residuals(y_u), Equation = "Unemployment-change outcome"),
  data.frame(date = train$date, e = residuals(x_u), Equation = "Claims-growth outcome"))
ggplot(residual_paths, aes(date, e)) + geom_line(colour = "#315375") +
  facet_wrap(~ Equation, ncol = 1, scales = "free_y") +
  labs(x = NULL, y = "Fitted residual", title = "Residual behavior is part of the evidence")

Separate panels of unrestricted-model fitted residuals over time for the unemployment-change and claims-growth outcome equations.

A limited sensitivity table, not a search for significance

Three specifications fixed in advance, $(p,q)=(3,1),(6,1),(6,3)$ , with $p$ the outcome’s own lags and $q$ the other variable’s lags, all on the same 336 rows.

Predict: does the reverse-direction decision at 5% (a) stay the same across all three rows, or (b) depend on the lag menu?

menu <- data.frame(p = c(3, 6, 6), q = c(1, 1, 3))
sensitivity <- data.frame()
for (outcome in c("y", "x")) {
  predictor <- if (outcome == "y") "x" else "y"
  for (i in seq_len(nrow(menu))) {
    own <- paste0(outcome, "_l", seq_len(menu$p[i]))
    other <- paste0(predictor, "_l", seq_len(menu$q[i]))
    r <- lm(reformulate(own, response = outcome), data = train, na.action = na.fail)
    u_fit <- lm(reformulate(c(own, other), response = outcome), data = train, na.action = na.fail)
    result <- granger_nested_test(r, u_fit)
    diagnostic <- granger_bg_test(u_fit, order = 12)
    sensitivity <- rbind(sensitivity, data.frame(
      Outcome = outcome, p = menu$p[i], q = menu$q[i], n = result$n,
      F = result$F, Nominal_F_p = result$p_value, BG_12_p = diagnostic$p_value))
  }
}
knitr::kable(transform(sensitivity,
  Nominal_F_p = format.pval(Nominal_F_p, digits = 3, eps = 1e-6),
  BG_12_p = format.pval(BG_12_p, digits = 3, eps = 1e-6)), digits = 3)

Outcome	p	q	n	F	Nominal_F_p	BG_12_p
y	3	1	336	7.542	0.006356	<1e-06
y	6	1	336	11.080	0.000972	<1e-06
y	6	3	336	21.058	< 1e-06	0.0031
x	3	1	336	7.922	0.005177	0.1363
x	6	1	336	4.835	0.028578	0.2282
x	6	3	336	1.804	0.146263	0.2549

rev_p <- sensitivity$Nominal_F_p[sensitivity$Outcome == "x"]
fwd_bg <- sensitivity$BG_12_p[sensitivity$Outcome == "y"]

Reverse direction, nominal p-values: 0.00518 at (3,1), 0.0286 at (6,1), and 0.146 at the main (6,3). Whether unemployment history appears to help predict claims depends on the finite lag menu. In the forward direction every row rejects, but the BG(12) p-values are 1.09e-13, 3.5e-07 and 0.0031: the shorter specifications fail the residual check far more decisively.

Hamilton (§11.2) reports the same lag-length sensitivity in the money–income literature, and Enders (Ch. 5 §8) frames the trade-off: too few lags misspecify, too many waste degrees of freedom. The remedy is to fix lags in advance by a stated rule on a common sample — not to report the row you like. Stop here. If residual dependence remains, revisit the economic specification as a separate modeling exercise.

Compare forecasts separately

Freeze the training coefficients of the restricted (6 own lags) and unrestricted (6 own, 3 claims lags) unemployment models. At each 2018–2019 origin, newly observed previous-month values supply the lags: a sequence of one-step predictions, not a 24-step path from December 2017.

Predict: recall section 3’s question. Must the claims model win this 24-month holdout because $RSS_U\le RSS_R$ in training — (a) yes, (b) no?

predictions <- data.frame(
  Own_history = as.numeric(predict(y_r, newdata = test)),
  Add_claims = as.numeric(predict(y_u, newdata = test)),
  No_change = rep(0, nrow(test)))
scores <- data.frame()
forecast_paths <- data.frame()
for (model in names(predictions)) {
  error <- test$y - predictions[[model]]
  level_prediction <- test$U_previous + predictions[[model]]
  stopifnot(isTRUE(all.equal(error, test$U - level_prediction)))
  scores <- rbind(scores, data.frame(Model = model, n = length(error),
    RMSE = sqrt(mean(error^2)), MAE = mean(abs(error))))
  forecast_paths <- rbind(forecast_paths,
    data.frame(date = test$date, U = level_prediction, Model = model))
}
knitr::kable(scores, digits = 4)

Model	n	RMSE	MAE
Own_history	24	0.1129	0.0929
Add_claims	24	0.1057	0.0831
No_change	24	0.1099	0.0875

ggplot(forecast_paths, aes(date, U, colour = Model, linetype = Model)) +
  geom_line() + geom_line(data = test, aes(date, U), inherit.aes = FALSE,
                          linewidth = 0.9, colour = "black") +
  labs(x = NULL, y = "Unemployment rate (percent)",
       title = "One-step index predictions at post-release origins, 2018-2019",
       subtitle = "Black line: observed unemployment; coefficients fixed after training") +
  scale_x_date(breaks = as.Date(c("2018-01-01", "2018-07-01", "2019-01-01", "2019-07-01")),
               date_labels = "%b %Y") +
  theme(legend.position = "bottom")

Observed unemployment in black and three one-step prediction series from 2018 to 2019: own-history, added-claims, and no-change forecasts, using fixed training coefficients.

The level prediction is $\tilde U_{t\mid\mathcal I_t^{\mathrm{release}}}=U_{t-1}+\tilde y_{t\mid\mathcal I_t^{\mathrm{release}}}$ , so change and level errors are identical, and RMSE and MAE are in percentage points. The added-claims RMSE is 0.1057 against 0.1129 for own history and 0.1099 for no change: a modest realized improvement on 24 particular months. That is descriptive. It does not repair the training diagnostics, identify an intervention effect, or promise a future gain, and we attach no ordinary Diebold–Mariano p-value to these nested estimated models.

Keep three conclusions separate: the nominal restriction tests (arithmetic plus limitations), the diagnostics and sensitivity (what the regressions leave and how much results lean on lags), and the held-out scores (realized accuracy on revised data at stated origins).

Write your report. Using the notes’ reporting template (§9.10), draft a short paragraph for the forward direction in your own words: transformed outcome and units, frequency and dates, lags and information, the excluded block, $n$ , F and degrees of freedom, nominal p-value, diagnostics, a bounded conclusion, and the separate holdout result. “We computed the test but cannot give it the usual interpretation under this specification” is a legitimate finding.

7. Uses and limits

Integrated levels. Our clean baseline is a stationary dynamic regression. With integrated levels, lag-exclusion statistics can have nonstandard reference distributions, so a routine p-value does not settle anything. “Difference everything” is not a universal remedy either: two wandering levels may share a stationary long-run relationship whose adjustment carries predictive information. Recognizing this is Core; the machinery is Module 10’s.

Economics and beyond. Sims (1972) asked whether money helps predict income and the reverse. Friedman and Kuttner’s evidence, discussed by Enders, shows the money–income result changing across sample periods under a fixed specification. In brain-signal research the same ARDL comparison asks whether region A’s past helps predict region B; slow, indirect fMRI measurement complicates reading a predictive direction as neural influence. In each case, state the series, the outcome, the other histories, the lags and the period.

Research connection. Cornwall, Mills, Sauley and Weng (2019) find that some in-sample Phillips-curve rejections coexist with weak held-out predictive evidence. Ask what was predicted, against which benchmark, on which observations, with what loss and what information.

8. Looking Ahead — optional, outside Core

Starting from an ARDL $(1,1)$ with current $x_t$ , $y_t=\alpha+\phi_1y_{t-1}+\beta_0x_t+\beta_1x_{t-1}+\epsilon_t$ , subtracting $y_{t-1}$ and collecting terms gives, for $\phi_1\ne1$ ,

$\Delta y_t=\lambda\bigl(y_{t-1}-c-\kappa x_{t-1}\bigr)+\beta_0\Delta x_t+\epsilon_t,$

with $\lambda=\phi_1-1$ , $c=\alpha/(1-\phi_1)$ and $\kappa=(\beta_0+\beta_1)/(1-\phi_1)$ . This is algebra, not proof of cointegration; Module 10 develops the evidence. Ask the tutor to build any preview on request, labelled as such.

9. Explorations — optional, guided to open

These are outside the assessable surface: starting points, not required questions or official answers.

Size and sample length. Rerun section 5’s null-rejection loop with the retained sample cut to 60 and 150 observations. Predict first: does the 5% rejection rate drift? Report Monte Carlo standard errors.
Power curve. In section 4, replace 0.6 by 0.05, 0.1 and 0.2 and repeat each over a few hundred seeded samples at two sample sizes. Plot rejection rate against the coefficient. Which matters more here, coefficient or $n$ ?
Aggregation moves the lag. Simulate a weekly process in which $x$ affects $y$ with a one-week lag. Average both to four-week “months.” Run the past-only monthly test and a monthly regression that includes current $x$ . Where did the effect go? Connect your answer to instantaneous causality.
A persistent hidden driver. Make $z_t$ an AR(1) in section 5. Does including only $z_{t-1}$ still remove $x$ ’s predictive content with one lag? State the population target before running.
Levels versus changes. Simulate two independent random walks, run a Granger F-test in levels with own lags, and record the 5% rejection rate over many seeds. Then repeat in differences. Label this a Module 10 preview.
Cross-sample validation, one split. On the training data, pick one $\tau$ , fit the restricted and unrestricted unemployment models on each side and predict the other side. Sum the squared prediction errors and form the Ashley–Tsang pseudo-F. Can it be negative? Try several $\tau$ and explain why a single split is arbitrary.

10. Run the referee and reconcile

Read CHECKS.md for the scope of these checks. They verify shipped examples and invariants; they cannot certify a causal interpretation or every modification. Running them in a private environment leaves your objects unchanged.

source("R/checks.R", local = new.env())

Module 9 checks
---------------
PASS 1. Complete Saturday calendar, 4-5 weeks per month, mean weekly aggregation, dated lags, 336 training and 24 evaluation rows  [weeks=1617; months merged=372; training=336; evaluation=24]
PASS 2. Manual F = anova F = helper F on common rows; q=3, k_U=10, df=326; RSS_U <= RSS_R; unequal rows refused  [manual=21.058129; anova=21.058129; helper=21.058129]
PASS 3. Known one-way DGP: lab draw's cross-lag estimate rounds to 0.600 and rejects; large sample near 0.6  [lab estimate=0.6001; large-n estimate=0.6040]
PASS 4. Common driver: pairwise signal rejects; kept conditional draw is a false rejection (p near 0.0157); population targets 0.4 and 0  [conditional p=0.01566; large-n pairwise=0.4018; conditional=0.0139]
PASS 5. Application: forward F and BG(12) failure, reverse nonrejection, reverse decision depends on lags  [forward F=21.0581; BG p=0.00310; reverse p=0.1463]
PASS 6. Frozen one-step predictions use only lagged inputs; change and level errors match; holdout RMSEs reproduce  [RMSE=0.112942/0.105747/0.109924]
---------------
6/6 checks passed.

Finish your private learning-log entry. Record your revised reasoning and remaining uncertainty in your own words. If your prediction was wrong, say which information set, lag choice, or timing assumption changed your mind.