Module 8 Lab: Autoregressive Distributed Lag Models

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: pictures → exogeneity → controlled experiment → dynamic responses → forecasting calendar → real data. The notes are canonical for facts, notation and the Core boundary. Sections 2–7 are the main path; sections 8–9 are optional. The model expands the mixing board with observed predictors, retaining outcome lags and switching off explicit MA terms in the baseline. A fitted value is ŷt\hat y_t; a forecast is ỹt+h∣t\tilde y_{t+h\mid t}.

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_08 and dir.exists("data") should be TRUE. Quarto’s file render uses the lab’s directory automatically. Installed R packages are ggplot2, forecast, urca, knitr and rmarkdown; rendering also needs Quarto. Render this file by path, not the whole folder.

2. Draw the information before estimating the equation

An ARMA remembers outcomes and innovations. An ARDL brings an observed second series into the room. Predict: if current xx is unavailable, which route must a forecast replace or omit: (a) outcome memory, (b) current-xx input, or (c) every lagged-xx input? Draw the restriction before reading on.

ARMA: memory of outcomes and innovations Past innovation epsilon t minus 1 enters past outcome y t minus 1 and current outcome y t. Past outcome enters current outcome. Current innovation epsilon t enters current outcome. ε(t−1)ε(t)y(t−1)y(t) MA memoryAR memory
Two kinds of memory in an ARMA model. Innovations are unobserved DGP quantities, not estimated residuals.
ARDL: observed predictor history and contemporaneous information Lagged x enters current x along a dashed contextual arrow, and enters current y. Current x, lagged y, and current innovation also enter current y. x(t−1)x(t)y(t−1)y(t)ε(t) x dynamics: contextcurrent xlagged xAR memoryPrevious periodCurrent period
An ARDL outcome equation. The dashed arrow supplies context about x's persistence; this one equation does not model x's own evolution.

These are posited dependencies, not identified causal effects. Unrolling time keeps yt−1y_{t-1} and yty_t separate and lets feedback across dates remain acyclic.

The ARDL(1,1)(1,1) equation is

yt=α+ϕ1yt−1+β0xt+β1xt−1+ϵt.y_t=\alpha+\phi_1y_{t-1}+\beta_0x_t+\beta_1x_{t-1}+\epsilon_t.

Generalize to pp outcome lags and predictor lags 0,…,q0,\ldots,q only after you can point to each term. Here qq is maximum predictor lag; ARMA uses qq for MA order. Setting β0=0\beta_0=0 removes current xx, not its history. Setting all βℓ=0\beta_\ell=0 gives an AR. All ϕj=0\phi_j=0 gives a distributed-lag regression.

With (α,ϕ1,β0,β1)=(0.2,0.6,0.4,0.2)(\alpha,\phi_1,\beta_0,\beta_1)=(0.2,0.6,0.4,0.2) and (yt−1,xt,xt−1)=(2,1,0.5)(y_{t-1},x_t,x_{t-1})=(2,1,0.5), compute the conditional mean:

0.2 + 0.6 * 2 + 0.4 * 1 + 0.2 * 0.5
[1] 1.9

Calling this a conditional mean requires the disturbance to have mean zero given those inputs. A realization also contains ϵt\epsilon_t; an estimated equation replaces the coefficients with hats and the disturbance with ete_t.

A software translation to check. Regression with AR(1) disturbances, yt=α+βxt+uty_t=\alpha+\beta x_t+u_t, ut=ϕ1ut−1+ϵtu_t=\phi_1u_{t-1}+\epsilon_t, implies

yt=(1−ϕ1)α+ϕ1yt−1+βxt−ϕ1βxt−1+ϵt.y_t=(1-\phi_1)\alpha+\phi_1 y_{t-1}+\beta x_t -\phi_1\beta x_{t-1}+\epsilon_t.

It is a restricted ARDL. xreg does not by itself tell you which specification was fitted. Modify on paper: allow xt−2x_{t-2}; then set β0=0\beta_0=0. Explain which paths survive and which information a one-step forecast needs.

3. The Office fire drill: who else heard the news?

Think of two coworkers reacting to the same emergency in The Office fire drill. One moves and the other moves; both may be reacting to the alarm. The diagram makes that omitted route visible without needing a clip or a network connection. The reference is an analogy for common news, not a claim that the scene is a controlled experiment.

Two different concerns: common current news and lagged feedback Lagged y points to current x and current y. Lagged x points to current x and current y. Current x points to current y. Common current news points to both current x and current y, making a source of current disturbance dependence visible. x(t−1)x(t)y(t−1)y(t)Common current news Lagged feedback: y(t−1) → x(t)
Common news can connect current x with the outcome disturbance. Lagged feedback is a different route. The diagram depicts possible mechanisms, not identified effects.

Let ℱt−1\mathcal F_{t-1} contain the relevant histories. Our sufficient conditional-mean restriction is

E[ϵt∣ℱt−1,xt]=0.E[\epsilon_t\mid\mathcal F_{t-1},x_t]=0.

It is about the whole conditioning information set, not just a sequence of pairwise correlations. Alongside suitable sampling, stability, finite moments and full rank, it supports consistency. It does not promise exact finite-sample unbiasedness or exact dynamic-regression t/Ft/F inference.

Predict: which situation directly threatens the structural coefficient’s OLS interpretation: (a) xtx_t is correlated with yty_t, (b) current xx reveals part of the structural disturbance, or (c) xtx_t uses already-conditioned-on yt−1y_{t-1} and the current shocks are independent of each other and the past?

The current-shock and feedback stories differ. In the feedback case, the path ϵt→yt→xt+1\epsilon_t\to y_t\to x_{t+1} reveals why future xx can contain information about today’s disturbance while the sequential restriction still holds. Strict exogeneity over the entire xx path is stronger.

White noise need not be unrelated to current xx. Lagging a variable does not remove every persistent omitted driver. A clean residual ACF and robust standard errors do not certify structural exogeneity. Keep three questions separate: conditional prediction, response along a stated input path, and causal intervention effect.

4. Build both series, then break the moment condition

The pair of equations is

xt=0.5xt−1+cyt−1+ϵx,t,yt=0.2+0.6yt−1+0.4xt+0.2xt−1+ϵt.x_t=0.5x_{t-1}+c y_{t-1}+\epsilon_{x,t},\qquad y_t=0.2+0.6y_{t-1}+0.4x_t+0.2x_{t-1}+\epsilon_t.

The shock pairs are independent across dates, jointly Gaussian with unit variances and current correlation rr. Compare (c,r)=(0,0)(c,r)=(0,0), (0,0.7)(0,0.7), and (0.2,0)(0.2,0). These three configurations have stable joint dynamics. In a general feedback system, checking only |ϕ1|<1|\phi_1|<1 would be insufficient.

Predict before running this section: with more data, does the shared-shock estimate of the current-xx coefficient approach (a) 0.4, (b) 1.1, or (c) zero? What changes when the only new feature is lagged feedback?

ardl_simulator <- function(n = 400, alpha = 0.2, phi = 0.6,
                           beta = c(0.4, 0.2), rho_x = 0.5,
                           feedback = 0, shock_corr = 0, burn_in = 200) {
  stopifnot(n >= 2, n == as.integer(n), burn_in >= 0,
            burn_in == as.integer(burn_in), length(beta) == 2,
            all(is.finite(c(alpha, phi, beta, rho_x, feedback, shock_corr))),
            abs(shock_corr) < 1)
  total <- n + burn_in + 1L
  epsilon_x <- rnorm(total)
  epsilon <- shock_corr * epsilon_x +
    sqrt(1 - shock_corr^2) * rnorm(total)
  x <- y <- numeric(total)
  for (t in 2:total) {
    x[t] <- rho_x * x[t - 1] + feedback * y[t - 1] + epsilon_x[t]
    y[t] <- alpha + phi * y[t - 1] + beta[1] * x[t] +
      beta[2] * x[t - 1] + epsilon[t]
  }
  keep <- (burn_in + 2L):total
  data.frame(t = seq_len(n), y = y[keep], x = x[keep],
             epsilon = epsilon[keep], epsilon_x = epsilon_x[keep])
}

The for-loop is the DGP: compute xtx_t first, then yty_t. Burn-in reduces the effect of arbitrary initial values in a stable system; it cannot create stability. This function and the multiplier function below are built inline, with exact canonical copies in R/ardl_utils.R for reuse.

case_names <- c("Baseline", "Shared current shocks", "Lagged feedback only")
case_corr <- c(0, 0.7, 0)
case_feedback <- c(0, 0, 0.2)
case_data <- case_fits <- vector("list", 3)
single_estimates <- matrix(NA_real_, 3, 4,
                           dimnames = list(case_names,
                                           c("alpha", "phi_1", "beta_0", "beta_1")))
for (s in seq_along(case_names)) {
  set.seed(637608)  # Same underlying independent draws across the three cases.
  d <- ardl_simulator(shock_corr = case_corr[s], feedback = case_feedback[s])
  d$y_lag1 <- c(NA, head(d$y, -1))
  d$x_lag1 <- c(NA, head(d$x, -1))
  fit <- lm(y ~ y_lag1 + x + x_lag1, data = d)
  case_data[[s]] <- d
  case_fits[[s]] <- fit
  single_estimates[s, ] <- coef(fit)
}
knitr::kable(rbind("DGP coefficients" = c(0.2, 0.6, 0.4, 0.2),
                   single_estimates), digits = 3)
alpha phi_1 beta_0 beta_1
DGP coefficients 0.200 0.600 0.400 0.200
Baseline 0.221 0.555 0.352 0.276
Shared current shocks 0.217 0.565 1.065 -0.068
Lagged feedback only 0.265 0.552 0.349 0.259

Lags are explicit shifted columns. stats::lag() on a numeric vector is not a substitute for checking those columns against their dates. One sample is not the distribution; repeat the experiment serially with a fixed seed.

set.seed(637609)
B <- 200
mc <- array(NA_real_, c(B, 4, 3),
            dimnames = list(NULL, colnames(single_estimates), case_names))
for (b in seq_len(B)) {
  for (s in seq_along(case_names)) {
    d <- ardl_simulator(shock_corr = case_corr[s], feedback = case_feedback[s])
    d$y_lag1 <- c(NA, head(d$y, -1))
    d$x_lag1 <- c(NA, head(d$x, -1))
    mc[b, , s] <- coef(lm(y ~ y_lag1 + x + x_lag1, data = d))
  }
}
mc_table <- data.frame()
truth <- c(0.2, 0.6, 0.4, 0.2)
for (s in seq_along(case_names)) {
  for (j in seq_along(truth)) {
    mc_table <- rbind(mc_table, data.frame(
      case = case_names[s], coefficient = colnames(single_estimates)[j],
      DGP = truth[j], mean_estimate = mean(mc[, j, s]),
      simulation_sd = sd(mc[, j, s]),
      MCSE_of_mean = sd(mc[, j, s]) / sqrt(B)))
  }
}
knitr::kable(mc_table, digits = 3)
case coefficient DGP mean_estimate simulation_sd MCSE_of_mean
Baseline alpha 0.2 0.204 0.053 0.004
Baseline phi_1 0.6 0.599 0.033 0.002
Baseline beta_0 0.4 0.396 0.053 0.004
Baseline beta_1 0.2 0.206 0.055 0.004
Shared current shocks alpha 0.2 0.205 0.039 0.003
Shared current shocks phi_1 0.6 0.597 0.028 0.002
Shared current shocks beta_0 0.4 1.097 0.036 0.003
Shared current shocks beta_1 0.2 -0.144 0.055 0.004
Lagged feedback only alpha 0.2 0.206 0.058 0.004
Lagged feedback only phi_1 0.6 0.595 0.034 0.002
Lagged feedback only beta_0 0.4 0.403 0.042 0.003
Lagged feedback only beta_1 0.2 0.200 0.053 0.004 
mc_plot <- data.frame(estimate = as.vector(mc[, "beta_0", ]),
                      case = rep(case_names, each = B))
ggplot(mc_plot, aes(x = estimate, fill = case)) +
  geom_histogram(bins = 25, color = "white") +
  geom_vline(xintercept = 0.4, linetype = 2, linewidth = 0.8) +
  facet_wrap(~ case, ncol = 1, scales = "free_y") +
  labs(x = "Estimated contemporaneous coefficient", y = "Replications",
       title = "Same structural beta_0; different regressor–disturbance relations",
       subtitle = "Dashed line: structural DGP coefficient 0.4") +
  guides(fill = "none")

The Monte Carlo SD measures dispersion across samples. Its MCSE measures precision of the simulation mean, not uncertainty about the structural interpretation. Inspect the table before declaring which discrepancy is sampling error and which persists with sample size.

In the shared-shock case, E[ϵt∣ϵx,t]=0.7ϵx,tE[\epsilon_t\mid\epsilon_{x,t}]=0.7\epsilon_{x,t} and ϵx,t=xt−0.5xt−1\epsilon_{x,t}=x_t-0.5x_{t-1}. Substitution gives

E[yt∣ℱt−1,xt]=0.2+0.6yt−1+1.1xt−0.15xt−1.E[y_t\mid\mathcal F_{t-1},x_t] =0.2+0.6y_{t-1}+1.1x_t-0.15x_{t-1}.

The known structural coefficients are still 0.4 and 0.2. OLS is estimating a different conditional relationship. Here ϵt\epsilon_t is a structural disturbance, not an innovation relative to current xx and the history. The remaining conditional disturbance is ϵt−0.7ϵx,t\epsilon_t-0.7\epsilon_{x,t}.

forecast::ggAcf(residuals(case_fits[[2]]), lag.max = 20) +
  labs(title = "Residual ACF: the shared-shock equation can look well specified",
       subtitle = "A residual diagnostic cannot certify a structural interpretation")

Modify: raise n to 2,000 in ardl_simulator() calls and repeat the comparison, then restore shock_corr = 0 and use feedback = 0.2. Explain why one change sharpens a different target and the other preserves the structural target under the stated independent-shock assumptions.

5. Follow every path: pulse, step, and predictor innovation

Hold initial conditions and outcome innovations fixed. A pulse raises one xtx_t by one then returns xx to its baseline path. A step raises all future xx values by one. The direct lagged coefficient is only one of tomorrow’s routes:

m0=β0,m1=ϕ1β0+β1,mh=ϕ1mh−1(h≥2).m_0=\beta_0,\qquad m_1=\phi_1\beta_0+\beta_1,\qquad m_h=\phi_1m_{h-1}\quad(h\ge2).

Predict: with (ϕ1,β0,β1)=(0.6,0.4,0.2)(\phi_1,\beta_0,\beta_1)=(0.6,0.4,0.2), is tomorrow’s imposed-pulse response (a) 0.2, (b) 0.44, or (c) 0.84? Which number belongs to the permanent step instead?

ardl_multipliers <- function(phi, beta, horizon = 20) {
  stopifnot(is.numeric(phi), is.numeric(beta), length(beta) >= 1,
            all(is.finite(c(phi, beta))), horizon >= 0,
            horizon == as.integer(horizon))
  pulse <- numeric(horizon + 1L)
  for (h in 0:horizon) {
    direct <- if (h < length(beta)) beta[h + 1L] else 0
    inherited <- 0
    if (h > 0 && length(phi) > 0) {
      for (j in seq_len(min(length(phi), h))) {
        inherited <- inherited + phi[j] * pulse[h - j + 1L]
      }
    }
    pulse[h + 1L] <- direct + inherited
  }
  data.frame(h = 0:horizon, pulse = pulse, step = cumsum(pulse))
}
response <- ardl_multipliers(phi = 0.6, beta = c(0.4, 0.2), horizon = 20)
knitr::kable(head(response, 5), digits = 4)
h pulse step
0 0.4000 0.4000
1 0.4400 0.8400
2 0.2640 1.1040
3 0.1584 1.2624
4 0.0950 1.3574 
response_long <- rbind(
  data.frame(h = response$h, response = response$pulse, experiment = "Pulse"),
  data.frame(h = response$h, response = response$step, experiment = "Permanent step"))
ggplot(response_long, aes(h, response, color = experiment)) +
  geom_line(linewidth = 0.9) + geom_point(size = 1.7) +
  labs(x = "Periods since the imposed change", y = "Difference in y",
       color = "Input path", title = "One equation, two different input experiments") +
  theme(legend.position = "bottom")

The cumulative pulse response equals the step response in this linear, time-invariant equation. Under stable outcome dynamics its limit is

∑ℓ=0qβℓ1−∑j=1pϕj=1.5.\frac{\sum_{\ell=0}^{q}\beta_\ell}{1-\sum_{j=1}^{p}\phi_j}=1.5.

This is a change, not the mean of yy. Estimated counterparts need hats and uncertainty; a near-zero denominator magnifies coefficient uncertainty.

Break stability

Before running: does a finite algebraic ratio imply convergence when ϕ1=1.05\phi_1=1.05? Follow the response, not just the denominator.

stable <- ardl_multipliers(0.6, c(0.4, 0.2), horizon = 20)
unit_root <- ardl_multipliers(1, c(0.4, 0.2), horizon = 20)
unstable <- ardl_multipliers(1.05, c(0.4, 0.2), horizon = 20)
break_plot <- rbind(
  data.frame(h = stable$h, step = stable$step, phi = "phi = 0.6"),
  data.frame(h = unit_root$h, step = unit_root$step, phi = "phi = 1"),
  data.frame(h = unstable$h, step = unstable$step, phi = "phi = 1.05"))
ggplot(break_plot, aes(h, step, color = phi)) +
  geom_line(linewidth = 1) +
  labs(x = "Periods since permanent step", y = "Difference in y",
       color = NULL, title = "A formula needs its stability condition") +
  theme(legend.position = "bottom")

The unit-root step accumulates without bound; the explosive step also fails to converge. The finite ratio 0.6/(1−1.05)=−120.6/(1-1.05)=-12 is not an eventual response. The helper intentionally permits unstable inputs so this failure is visible.

Change the input experiment

Now use the baseline with independent current shocks, no lagged feedback, and fixed outcome innovations. An isolated innovation in the predictor equation generates xx differences 1,0.5,0.25,…1,0.5,0.25,\ldots, not a one-period pulse. Compute the outcome-path differences directly:

x_difference <- 0.5^(0:20)
y_difference <- numeric(length(x_difference))
for (k in seq_along(x_difference)) {
  inherited <- if (k == 1) 0 else 0.6 * y_difference[k - 1]
  lagged_input <- if (k == 1) 0 else 0.2 * x_difference[k - 1]
  y_difference[k] <- inherited + 0.4 * x_difference[k] + lagged_input
}
knitr::kable(data.frame(h = 0:3, imposed_pulse = response$pulse[1:4],
                        x_innovation_path = y_difference[1:4]), digits = 4)
h imposed_pulse x_innovation_path
0 0.4000 0.4000
1 0.4400 0.6400
2 0.2640 0.5840
3 0.1584 0.4504

At horizon one the innovation-induced response is 0.64, compared with 0.44 for an imposed pulse. With shared shocks, conditioning on xx news can change the expected outcome disturbance too. With feedback, follow both equations. None of these calculations alone identifies an intervention in observational data.

Modify: use ϕ1=−0.5\phi_1=-0.5 in the pulse experiment. Explain the cancellation at horizon one by following both paths. Then compare a stable value near one with the explosive value and say why a large multiplier needs uncertainty.

6. A predictor needs a calendar

A current predictor in an outcome equation becomes a future predictor in a forecast. Under the maintained model,

ỹt+1∣t=α̂+ϕ̂1yt+β̂0x̃t+1∣t+β̂1xt.\tilde y_{t+1\mid t}=\hat\alpha+\hat\phi_1y_t +\hat\beta_0\tilde x_{t+1\mid t}+\hat\beta_1x_t.

The predictor can be forecast, stipulated as a scenario, or observed early enough for a nowcast. Each changes the information set. Conditional scenario uncertainty omits predictor-path uncertainty by design. An unconditional forecast must account for predictor forecast errors and any relevant covariance with outcome shocks, as well as coefficient uncertainty.

Predict: does the method allowed to observe the target-period predictor answer the same forecasting question as a method that must forecast it?

set.seed(637610)
forecast_data <- ardl_simulator(n = 800)
forecast_data$y_lag1 <- c(NA, head(forecast_data$y, -1))
forecast_data$x_lag1 <- c(NA, head(forecast_data$x, -1))
train <- forecast_data[2:600, ]
holdout <- forecast_data[601:800, ]
forecast_fit <- lm(y ~ y_lag1 + x + x_lag1, data = train)
x_fit <- lm(x ~ x_lag1, data = train)
feasible_inputs <- holdout
feasible_inputs$x <- predict(x_fit, newdata = holdout)
feasible <- predict(forecast_fit, newdata = feasible_inputs)
oracle <- predict(forecast_fit, newdata = holdout)  # Knows target-period x.
knitr::kable(data.frame(
  information = c("Forecast x from its observed lag", "Oracle: observe target-period x"),
  RMSE = c(sqrt(mean((holdout$y - feasible)^2)),
           sqrt(mean((holdout$y - oracle)^2)))), digits = 3)
information RMSE
Forecast x from its observed lag 1.039
Oracle: observe target-period x 0.950

Both methods freeze training coefficients and roll observed lags. The oracle has future information. Its score could be a legitimate nowcast only if that predictor is actually released before the target at the stated issue date. The current contents of a CSV do not answer that question.

Modify the design on paper: move the issue date earlier than the previous month’s unemployment release. List the lag columns that are no longer observed. Explain what must be forecast or omitted before any new scoring.

7. Real data: do claims add usable information?

Let UtU_t be the monthly unemployment rate in percent. Let CtC_t be the mean of weekly initial claims whose week-ending Saturdays fall in month tt. Both are seasonally adjusted. Model

yt=ΔUt,xt=100ΔlogCt.y_t=\Delta U_t,\qquad x_t=100\Delta\log C_t.

yy is in percentage points; xx is approximately percentage growth in the mean weekly claims count. A monthly mean is not a monthly total. A permanent step in this xx would mean permanently higher claims growth, not one fixed increase in the level of claims. We use the empirical model for prediction and diagnostics, not a causal claims intervention.

The origin contract: wait until the previous month’s unemployment and its final claims Saturday have both been released, before target-month unemployment is released. Ignore partial target-month claims. This is one step ahead by monthly observation index and ordinarily a nowcast of the current calendar month. Write the information set as ℐtrelease\mathcal I_t^{\mathrm{release}}.

The caches are revised historical values collected on different dates. This is revised-data pseudo-out-of-sample evaluation, not a real-time-vintage experiment. Actual deployment requires checking the actual release calendar.

Prepare and inspect the training data

app_u <- read.csv(file.path(data_dir, "UNRATE.csv"))
app_c <- read.csv(file.path(data_dir, "ICSA.csv"))
app_u$observation_date <- as.Date(app_u$observation_date)
app_c$observation_date <- as.Date(app_c$observation_date)
app_u <- app_u[app_u$observation_date >= as.Date("1988-12-01") &
                 app_u$observation_date <= as.Date("2019-12-01"), ]
app_c <- app_c[order(app_c$observation_date), ]
stopifnot(
  !anyNA(app_u), !anyNA(app_c), all(app_c$ICSA > 0),
  !anyDuplicated(app_c$observation_date),
  all(as.POSIXlt(app_c$observation_date)$wday == 6),
  identical(app_c$observation_date,
            seq(as.Date("1989-01-07"), as.Date("2019-12-28"), by = "week")),
  identical(app_u$observation_date,
            seq(as.Date("1988-12-01"), as.Date("2019-12-01"), by = "month"))
)
# December 1988 unemployment supplies January 1989's change;
# this makes twelve outcome lags available by January 1990.
app_u$y <- c(NA_real_, diff(app_u$UNRATE))
app_u$U_previous <- c(NA_real_, head(app_u$UNRATE, -1L))
app_u <- app_u[app_u$observation_date >= as.Date("1989-01-01"), ]
app_c$month <- as.Date(format(app_c$observation_date, "%Y-%m-01"))
app_month_c <- aggregate(ICSA ~ month, data = app_c, FUN = mean)
app_nweeks <- table(app_c$month)
stopifnot(all(app_nweeks %in% c(4L, 5L)))
app_data <- merge(app_u, app_month_c,
                  by.x = "observation_date", by.y = "month", all = TRUE)
names(app_data)[match(c("observation_date", "UNRATE", "ICSA"),
                      names(app_data))] <- c("date", "U", "C")
stopifnot(nrow(app_data) == 372L, !anyNA(app_data))
app_data$x <- c(NA_real_, 100 * diff(log(app_data$C)))
for (j in 1:12) {
  app_data[[paste0("y_l", j)]] <-
    c(rep(NA_real_, j), head(app_data$y, -j))
  if (j <= 3) {
    app_data[[paste0("x_l", j)]] <-
      c(rep(NA_real_, j), head(app_data$x, -j))
  }
}
app_data <- app_data[app_data$date >= as.Date("1990-01-01"), ]
stopifnot(nrow(app_data) == 360L, !anyNA(app_data))
app_train <- app_data[app_data$date <= as.Date("2017-12-01"), ]
app_test <- app_data[app_data$date >= as.Date("2018-01-01"), ]
stopifnot(nrow(app_train) == 336L, nrow(app_test) == 24L)

Training has 336 monthly observations, January 1990–December 2017. Hold out the 24 months of 2018–2019. Padding from 1989 makes all lag menus use the same training outcomes. The code checks complete calendars before it merges or computes changes; missing weeks are not zero and must not be averaged away.

app_paths <- rbind(
  data.frame(date = app_train$date, value = app_train$U,
             series = "Unemployment rate (percent)"),
  data.frame(date = app_train$date, value = app_train$C / 1000,
             series = "Mean weekly claims (thousands)"),
  data.frame(date = app_train$date, value = app_train$y,
             series = "y: change in unemployment (percentage points)"),
  data.frame(date = app_train$date, value = app_train$x,
             series = "x: 100 times log change in mean weekly claims")
)
ggplot(app_paths, aes(date, value)) +
  geom_line(colour = "#16697A", linewidth = 0.4) +
  facet_wrap(~ series, scales = "free_y", ncol = 1) +
  labs(x = NULL, y = NULL) + theme_bw()

Training data: levels describe the economic setting; changes are the modeled variables. Each panel has its own units and scale.

Predict from the transformed plots: will the evidence support a stationary working baseline? The ADF null is a unit root; the KPSS null here is stationarity around a constant. The settings below are fixed in advance.

app_stationarity <- data.frame()
for (app_name in c("y", "x")) {
  app_adf <- urca::ur.df(app_train[[app_name]], type = "drift",
                       lags = 12, selectlags = "BIC")
  app_kpss <- urca::ur.kpss(app_train[[app_name]], type = "mu", use.lag = 12)
  app_stationarity <- rbind(app_stationarity, data.frame(
    Series = app_name,
    ADF_statistic = unname(app_adf@teststat[1, "tau2"]),
    ADF_5pct = unname(app_adf@cval["tau2", "5pct"]),
    KPSS_statistic = unname(app_kpss@teststat),
    KPSS_5pct = unname(app_kpss@cval[1, "5pct"])
  ))
}
knitr::kable(app_stationarity, digits = 3)
Series ADF_statistic ADF_5pct KPSS_statistic KPSS_5pct
y -3.692 -2.87 0.149 0.463
x -12.795 -2.87 0.105 0.463

Here both ADF tests reject at 5% and neither KPSS test rejects. This supports a working specification. It does not prove stationarity, constant variance, exogeneity or stability over future regimes.

Fit, diagnose, revise once

Fit a lag-only ARDL, setting β0=0\beta_0=0, beside an AR. The initial menu has 1–3 own lags and 1–3 claims lags. In this code q = 0 labels the AR benchmark with no claims term; it is not the theoretical ARDL(p,0)(p,0) with current xx. BIC chooses within each family on the same training rows. The holdout chooses neither the lag orders nor the coefficient estimates.

app_grid <- expand.grid(p = 1:3, q = 0:3)
app_grid$family <- ifelse(app_grid$q == 0, "AR", "ARDL")
app_grid$AIC <- app_grid$BIC <- NA_real_
app_fits <- vector("list", nrow(app_grid))
for (i in seq_len(nrow(app_grid))) {
  app_terms <- paste0("y_l", seq_len(app_grid$p[i]))
  if (app_grid$q[i] > 0) {
    app_terms <- c(app_terms, paste0("x_l", seq_len(app_grid$q[i])))
  }
  app_fits[[i]] <- lm(reformulate(app_terms, response = "y"),
                      data = app_train, na.action = na.fail)
  app_grid$AIC[i] <- AIC(app_fits[[i]])
  app_grid$BIC[i] <- BIC(app_fits[[i]])
  stopifnot(nobs(app_fits[[i]]) == nrow(app_train))
}
app_ar_rows <- which(app_grid$family == "AR")
app_ardl_rows <- which(app_grid$family == "ARDL")
app_ar_row <- app_ar_rows[which.min(app_grid$BIC[app_ar_rows])]
app_ardl_row <- app_ardl_rows[which.min(app_grid$BIC[app_ardl_rows])]
app_ar <- app_fits[[app_ar_row]]
app_ardl <- app_fits[[app_ardl_row]]
app_grid$selected <- seq_len(nrow(app_grid)) %in%
  c(app_ar_row, app_ardl_row)
knitr::kable(app_grid[order(app_grid$family, app_grid$BIC),
                     c("family", "p", "q", "AIC", "BIC", "selected")],
             digits = 2, row.names = FALSE)
family p q AIC BIC selected
AR 3 0 -326.58 -307.50 TRUE
AR 2 0 -311.87 -296.60 FALSE
AR 1 0 -295.81 -284.36 FALSE
ARDL 3 3 -362.42 -331.88 TRUE
ARDL 2 3 -351.75 -325.03 FALSE
ARDL 1 3 -347.93 -325.02 FALSE
ARDL 3 2 -343.53 -316.81 FALSE
ARDL 3 1 -332.16 -309.25 FALSE
ARDL 2 2 -329.17 -306.27 FALSE
ARDL 1 2 -320.36 -301.27 FALSE
ARDL 2 1 -315.68 -296.59 FALSE
ARDL 1 1 -302.14 -286.87 FALSE 
app_residuals <- data.frame(date = app_train$date, e = residuals(app_ardl))
ggplot(app_residuals, aes(date, e)) +
  geom_hline(yintercept = 0, colour = "grey65") +
  geom_line(colour = "#16697A", linewidth = 0.4) +
  labs(x = NULL, y = "Residual (percentage points)") + theme_bw()

Training residuals from the selected ARDL. Inspect both their time pattern and their dependence.
forecast::ggAcf(app_residuals$e, lag.max = 24) + theme_bw()

Training residuals from the selected ARDL. Inspect both their time pattern and their dependence.

The Breusch–Godfrey-style auxiliary regression includes the original regressors and 12 residual lags, using the same auxiliary sample after dropping the first 12 rows. Its nR2nR^2 chi-squared reference is asymptotic and relies on maintained regression conditions, including the appropriate homoskedastic variance restriction; it is not robust to arbitrary conditional heteroskedasticity. Selection also makes this an exploratory diagnostic. It is a specification warning, not an exogeneity test.

app_bg_data <- app_train
app_bg_data$e <- as.numeric(residuals(app_ardl))
for (j in 1:12) {
  app_bg_data[[paste0("e_l", j)]] <-
    c(rep(NA_real_, j), head(app_bg_data$e, -j))
}
app_bg_data <- app_bg_data[13:nrow(app_bg_data), ]
app_bg_terms <- c(attr(terms(app_ardl), "term.labels"), paste0("e_l", 1:12))
app_bg_fit <- lm(reformulate(app_bg_terms, response = "e"), data = app_bg_data)
app_bg_stat <- nobs(app_bg_fit) * summary(app_bg_fit)$r.squared
app_bg_p <- pchisq(app_bg_stat, df = 12, lower.tail = FALSE)
knitr::kable(data.frame(Test = "BG, residual lags 1 through 12",
                       Statistic = app_bg_stat, df = 12,
                       p_value = format.pval(app_bg_p, digits = 3)),
             digits = 4)
Test Statistic df p_value
BG, residual lags 1 through 12 45.4938 12 8.48e-06

The short specification selects the edge of the lag menu and leaves serial dependence. Make one training-only revision: allow up to 12 own lags in both families, keeping at most 3 claims lags and the same 336 training outcomes.

app_short_ardl <- app_ardl
app_short_bg_p <- app_bg_p
app_expanded <- expand.grid(p = 1:12, q = 0:3)
app_expanded$family <- ifelse(app_expanded$q == 0, "AR", "ARDL")
app_expanded$BIC <- NA_real_
app_expanded_fits <- vector("list", nrow(app_expanded))
for (i in seq_len(nrow(app_expanded))) {
  app_terms <- paste0("y_l", seq_len(app_expanded$p[i]))
  if (app_expanded$q[i] > 0) {
    app_terms <- c(app_terms, paste0("x_l", seq_len(app_expanded$q[i])))
  }
  app_expanded_fits[[i]] <- lm(reformulate(app_terms, response = "y"),
                               data = app_train, na.action = na.fail)
  app_expanded$BIC[i] <- BIC(app_expanded_fits[[i]])
  stopifnot(nobs(app_expanded_fits[[i]]) == nrow(app_train))
}
app_exp_ar_rows <- which(app_expanded$family == "AR")
app_exp_ardl_rows <- which(app_expanded$family == "ARDL")
app_exp_ar_row <- app_exp_ar_rows[which.min(app_expanded$BIC[app_exp_ar_rows])]
app_exp_ardl_row <- app_exp_ardl_rows[which.min(app_expanded$BIC[app_exp_ardl_rows])]
app_ar <- app_expanded_fits[[app_exp_ar_row]]
app_ardl <- app_expanded_fits[[app_exp_ardl_row]]
# Show each family's three smallest training BIC values.
app_show_rows <- c(
  head(app_exp_ar_rows[order(app_expanded$BIC[app_exp_ar_rows])], 3),
  head(app_exp_ardl_rows[order(app_expanded$BIC[app_exp_ardl_rows])], 3)
)
knitr::kable(app_expanded[app_show_rows, c("family", "p", "q", "BIC")],
             digits = 2, row.names = FALSE)
family p q BIC
AR 6 0 -314.45
AR 5 0 -312.45
AR 7 0 -308.75
ARDL 6 3 -356.51
ARDL 7 3 -355.72
ARDL 8 3 -350.93 
app_bg_data$e <- as.numeric(residuals(app_ardl))[13:nrow(app_train)]
for (j in 1:12) {
  app_bg_data[[paste0("e_l", j)]] <-
    c(rep(NA_real_, j), head(as.numeric(residuals(app_ardl)), -j))[
      13:nrow(app_train)]
}
app_bg_terms <- c(attr(terms(app_ardl), "term.labels"), paste0("e_l", 1:12))
app_bg_fit <- lm(reformulate(app_bg_terms, response = "e"), data = app_bg_data)
app_bg_stat <- nobs(app_bg_fit) * summary(app_bg_fit)$r.squared
app_bg_p <- pchisq(app_bg_stat, df = 12, lower.tail = FALSE)
app_phi <- coef(app_ardl)[grep("^y_l", names(coef(app_ardl)))]
app_min_root <- min(Mod(polyroot(c(1, -app_phi))))
knitr::kable(data.frame(
  Diagnostic = c("Initial short ARDL: BG p-value", "Revised ARDL: BG p-value",
                 "Revised ARDL: smallest AR root modulus"),
  Value = c(format.pval(app_short_bg_p, digits = 3),
            format.pval(app_bg_p, digits = 3),
            format(round(app_min_root, 4), nsmall = 4))
), digits = 4)
Diagnostic Value
Initial short ARDL: BG p-value 8.48e-06
Revised ARDL: BG p-value 0.0031
Revised ARDL: smallest AR root modulus 1.1412 
forecast::ggAcf(residuals(app_ardl), lag.max = 24) + theme_bw()

The choices are AR(6) and lag-only ARDL (p=p= 6, q=q= 3). The ARDL’s smallest lag-polynomial root modulus is 1.141. Its revised BG pp-value is 0.0031; serial dependence remains detectable. Report that limitation and stop this illustration’s revision here. Do not keep searching until a diagnostic passes. HAC standard errors would not repair invalid moments or a misspecified mean.

app_revised_coef <- data.frame(Term = names(coef(app_ardl)),
                               Estimate = unname(coef(app_ardl)))
knitr::kable(app_revised_coef, digits = 4, row.names = FALSE)
Term Estimate
(Intercept) -0.0001
y_l1 -0.1646
y_l2 -0.0029
y_l3 0.0933
y_l4 0.0976
y_l5 0.2197
y_l6 0.2092
x_l1 0.0083
x_l2 0.0108
x_l3 0.0101

The x_l1 coefficient is an association between approximately one percentage point faster prior-month claims growth and the current unemployment-rate change, in percentage points, holding the other included lags fixed. It is not the effect of forcing people to file claims. Different conditioning sets can give different coefficient estimates.

Freeze the coefficients and roll observed information

Predict: does the ARDL have to beat the AR or no-change benchmark on this short holdout just because it has extra information? Commit before scoring.

Predict only from released lag columns. Reconstruct each level using Ũt∣ℐtrelease=Ut−1+ỹt∣ℐtrelease\tilde U_{t\mid\mathcal I_t^{\mathrm{release}}} =U_{t-1}+\tilde y_{t\mid\mathcal I_t^{\mathrm{release}}}. Do not accumulate all 24 changes from the December 2017 level; that would be a different multi-step exercise. No change in the level means a zero prediction for yty_t.

app_pred <- data.frame(
  AR = as.numeric(predict(app_ar, newdata = app_test)),
  ARDL = as.numeric(predict(app_ardl, newdata = app_test)),
  No_change = rep(0, nrow(app_test))
)
app_scores <- data.frame()
app_forecast_paths <- data.frame()
for (app_model in names(app_pred)) {
  app_error <- app_test$y - app_pred[[app_model]]
  app_level_forecast <- app_test$U_previous + app_pred[[app_model]]
  stopifnot(isTRUE(all.equal(app_error,
                            app_test$U - app_level_forecast)))
  app_scores <- rbind(app_scores, data.frame(
    Model = app_model, N = length(app_error),
    RMSE = sqrt(mean(app_error^2)), MAE = mean(abs(app_error))
  ))
  app_forecast_paths <- rbind(app_forecast_paths, data.frame(
    date = app_test$date, U = app_level_forecast, Model = app_model
  ))
}
knitr::kable(app_scores, digits = 4, row.names = FALSE)
Model N RMSE MAE
AR 24 0.1129 0.0929
ARDL 24 0.1057 0.0831
No_change 24 0.1099 0.0875 
ggplot(app_forecast_paths, aes(date, U, colour = Model, linetype = Model)) +
  geom_line(linewidth = 0.6) +
  geom_line(data = app_test, aes(date, U), inherit.aes = FALSE,
            colour = "black", linewidth = 0.8) +
  geom_point(data = app_test, aes(date, U), inherit.aes = FALSE,
             colour = "black", size = 1.2) +
  scale_colour_manual(values = c(AR = "#16697A", ARDL = "#D1495B",
                                 No_change = "#8A7B37")) +
  scale_x_date(breaks = as.Date(c("2018-01-01", "2018-07-01",
                                "2019-01-01", "2019-07-01")),
               date_labels = "%Y-%m") +
  labs(x = NULL, y = "Unemployment rate (percent)",
       subtitle = "Black line and points: actual unemployment rate") +
  theme_bw() + theme(legend.position = "bottom")

January 2018–December 2019: one-month-index-ahead predictions at the stated post-release origin, using fixed training coefficients.

On these dates the ARDL’s RMSE is 0.1057 percentage points; the AR’s is 0.1129 and the no-change benchmark’s is 0.1099. Change and reconstructed-level errors coincide because they add the same observed prior level to both sides. MAE is also in percentage points; MAPE is unsuitable for changes that can be zero or negative.

This ranking describes the chosen specifications and these dates. It does not establish exogeneity, causal effects, or a universal advantage of adding predictors. Do not retune on this holdout and still call it untouched.

Modify: state a different lag menu before running it, keep common samples and the same issue dates, and report the result as an additional exploratory comparison. Explain why it does not overturn the existing specification warning merely because the score improves.

8. Looking Ahead — three questions, optional

These previews are outside Module 8’s assessable Core and do not add a new required simulation/application cycle.

  • Granger causality: do past xx values add linear predictive information conditional on the included histories? In the first lag-only model, the null is β1=⋯=βq=0\beta_1=\cdots=\beta_q=0, jointly. Conditioning on current xx changes that question; predictive precedence is not causal identification.
  • Cointegration: can two I(1)I(1) series have an I(0)I(0) linear combination? An ARDL-to-ECM algebraic rewrite does not establish that property. Formal testing and error correction arrive in Module 10.
  • VAR: give both variables equations to model joint dynamics and feedback. Structural shock identification still needs additional assumptions.

Ask the tutor to build any of these on request, labelled as a preview.

9. Explorations — optional, guided to open

These investigations are outside the assessable surface. They are starting points, not extra required questions or official assessment answers.

  1. Sample size versus target. Repeat section 4 with n = 100, 400, and 2000. Plot dispersion and mean coefficient separately. Which discrepancy shrinks, and which target is the mean approaching?
  2. Follow the feedback path. Keep independent current shocks and vary feedback modestly around zero. Derive the joint transition matrix first and check its eigenvalues before claiming stationarity. Measure how ϵt\epsilon_t relates to xt+1x_{t+1} and reconcile with sequential exogeneity.
  3. Choose your input experiment. Compare an imposed pulse, a step, and an innovation-induced persistent xx path under no feedback and independent current shocks, holding outcome innovations fixed. Then label the extra assumptions needed to call any of them a causal intervention.
  4. Redesign the calendar. Propose an earlier issue date for the claims application. List available releases, replace unavailable regressors with justified forecasts, and evaluate without target-period leakage. Do not claim historical real-time performance from these revised caches.
  5. Invent a common driver. Add a persistent unobserved state to both equations. Choose parameters and state the intended target before running. Can lagged xx still correlate with the outcome disturbance? Design a comparison that shows why a dated lag is not an exogeneity certificate.

10. Run the referee and reconcile

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

source("R/checks.R", local = new.env())
Module 8 checks
---------------
PASS 1. Complete calendars, mean weekly aggregation, dated lags, 336 training and 24 holdout rows  [weeks=1617; training=336; holdout=24]
PASS 2. Simulator equations, seed reproducibility, burn_in=0 length and joint feedback stability  [max recursion gap=2.2e-16; largest eigenvalue=0.8868]
PASS 3. Structural targets in independent-shock cases; distinct conditional targets with shared shocks  [max target gap=0.0052; shared beta0 estimate=1.1037; remainder ACF1=-0.0029]
PASS 4. Pulse recursion, cumulative step, stable limit and unstable nonconvergence; persistent-x distinction  [pulse h1=0.44; persistent-x h1=0.64; stable limit=1.50]
PASS 5. Unavailable-predictor error contains beta0 epsilon_x; oracle changes the information set  [feasible-minus-oracle MSE=0.1535; population gap=0.16]
PASS 6. Common-sample AR(6)/ARDL(6,3), frozen lag-only predictions, matched errors and remaining BG warning  [RMSE=0.112942/0.105747/0.109924; BG p=0.003095; min root=1.1412]
---------------
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 assumption, information set, or input experiment changed your mind.