Module 8: Autoregressive Distributed Lag Models

Econ 6376 — Applied Time Series Econometrics

How to Use These Notes

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

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

If you missed lecture, read the main text, run the Core code, and complete the Core Practice problems. Then return to the optional sections that address your questions or interests. The previews of Granger causality and cointegration are orientation for Modules 9–10, not additional assessable material in Module 8.

Core Learning Objectives

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

1. Translate a time-unrolled diagram into an ARDL equation, distinguishing observed-variable lags, innovation lags, and contemporaneous relationships.
2. Write and interpret an ARDL$(p,q)$, including its restricted cases, and distinguish it from a regression with ARMA disturbances.
3. Explain the condition $E[\epsilon_t\mid\mathcal F_{t-1},x_t]=0$, why it involves the whole regressor information set, and how common current shocks differ from lagged feedback.
4. Use simulation to distinguish sampling variability from failure to recover structural coefficients, and explain why clean residuals do not establish exogeneity or causal identification.
5. Calculate impact, pulse, cumulative, and stable long-run multipliers. Distinguish an imposed path for $x$ from a shock to a persistent $x$ process and from a causal intervention.
6. Fit a stationary ARDL by OLS under the maintained assumptions, compare modest lag specifications on a common sample, and inspect regression residuals.
7. Construct a forecast using information available at its origin, identify when future predictor values must themselves be forecast, and explain the uncertainty a conditional forecast leaves out.
8. Evaluate an ARDL against univariate benchmarks on the same dates and explain the limitations of a forecast exercise using revised historical data.

---

8.1 Letting a Second Variable into the Room — Core

The Module 7 Handoff

Module 7 compared and combined forecasts made from a series’ own history. Combining forecasts can rearrange the information in that pool; it cannot supply information absent from every member. The next question is therefore concrete: could initial unemployment-insurance claims tell us something about unemployment that unemployment’s own history does not?

The autoregressive distributed lag model, abbreviated ARDL or ADL, gives us a first way to put that question into an equation. It also gives us a first encounter with the complications of multivariate time series: feedback, common shocks, reporting delays, and long-run relationships.

We begin with pictures. A directed acyclic graph (DAG) contains directed arrows and no path that follows those arrows back to its starting node. Here we unroll time: $y_{t-1}$ and $y_t$ are different nodes. That lets us depict feedback across dates without drawing a cycle within a date. The diagrams are small fragments of a longer time sequence.

What an ARMA Diagram Contains

For an ARMA$(1,1)$,

$$y_t=\alpha+\phi_1y_{t-1}+\theta_1\epsilon_{t-1}+\epsilon_t,$$

the outcome remembers past outcomes and past innovations. The following diagram shows those two routes into the current outcome. Deterministic terms are omitted from the picture.

Two kinds of memory in an ARMA model. Innovations are unobserved DGP quantities, not estimated residuals.

The arrow $\epsilon_{t-1}\to y_t$ corresponds to $\theta_1\epsilon_{t-1}$. The route $\epsilon_{t-1}\to y_{t-1}\to y_t$ also carries part of last period’s innovation forward. These are not competing definitions of memory; they are two routes through the equation.

What Changes in an ARDL Diagram

Now replace the explicit MA term with current and lagged observations of another series:

$$y_t=\alpha+\phi_1y_{t-1}+\beta_0x_t+\beta_1x_{t-1}+\epsilon_t.$$

An ARDL outcome equation. The dashed arrow supplies context about x's persistence; this one equation does not model x's own evolution.

The new information is observed $x$, not a past innovation. Its coefficient can differ across dates: $\beta_0$ need not equal $\beta_1$. In a monthly labor-market application, employers’ layoffs and unemployed workers’ job finding need not affect the unemployment rate on the same schedule.

Read the arrows carefully. These pictures represent posited dependencies in a model. Drawing $x_t\to y_t$ is not evidence that changing $x_t$ would cause a change in $y_t$. The within-period ordering is an assumption to investigate. We will make that issue central in §8.3.

Pause and predict. Which diagram term could still be present if the current value of $x$ were unavailable? Which arrow would disappear if the outcome had no own persistence? Write the corresponding restrictions before moving on.

The Mixing Board

The familiar master equation was

$$y_t=\alpha+\delta t+\sum_{j=1}^{p}\phi_jy_{t-j} +\sum_{\ell=1}^{q}\theta_\ell\epsilon_{t-\ell}+\epsilon_t.$$

We expand the board with a new observed-variable block. The baseline ARDL sets the explicit MA block and deterministic trend to zero:

Slider	Earlier univariate models	Module 8 baseline
Intercept $\alpha$	Available	On
Deterministic trend $\delta t$	Available	Off in our stationary examples
Outcome lags $\phi_j y_{t-j}$	On in AR/ARMA	On
Innovation lags $\theta_\ell\epsilon_{t-\ell}$	On in MA/ARMA	Off
Predictor block $\sum_{\ell=0}^{q}\beta_\ell x_{t-\ell}$	Absent	New
Current disturbance $\epsilon_t$	Present	Present; its relation to $x$ now needs attention

This module’s sequence is Concept → Math → Simulate → Real Data: the diagrams and equations establish the model; §8.3 states the assumptions; §8.4 deliberately breaks them; §§8.5–8.7 trace responses and apply the model. The short optional previews at the end do not repeat a full application workflow.

8.2 One Conditional Equation — Core

Start with ARDL(1,1)

In

$$y_t=\alpha+\phi_1y_{t-1}+\beta_0x_t+\beta_1x_{t-1}+\epsilon_t,$$

$\phi_1$ describes outcome persistence, $\beta_0$ is the contemporaneous $x$ coefficient, and $\beta_1$ is the direct lagged-$x$ contribution holding the other regressors fixed. The intercept is $\alpha$, not automatically the unconditional mean of $y$.

For example, if $\alpha=0.2$, $\phi_1=0.6$, $\beta_0=0.4$, and $\beta_1=0.2$, then $y_{t-1}=2$, $x_t=1$, and $x_{t-1}=0.5$ imply a conditional mean of $0.2+0.6(2)+0.4(1)+0.2(0.5)=1.9$, provided the remaining disturbance has conditional mean zero. The realization is $1.9+\epsilon_t$, not necessarily 1.9.

After estimation, we instead write

$$y_t=\hat\alpha+\hat\phi_1y_{t-1} +\hat\beta_0x_t+\hat\beta_1x_{t-1}+e_t.$$

The fitted value $\hat y_t$ is the estimated right-hand side without $e_t$. The residual is $e_t=y_t-\hat y_t$. Hats denote estimates; the true DGP parameters and $\epsilon_t$ remain unobserved outside a simulation. From Module 7 onward, an actual forecast carries a tilde, $\tilde y_{t+h\mid t}$.

Generalize Only After the Pieces Are Familiar

An ARDL$(p,q)$ with one predictor is

$$y_t=\alpha+\sum_{j=1}^{p}\phi_jy_{t-j} +\sum_{\ell=0}^{q}\beta_\ell x_{t-\ell}+\epsilon_t.$$

Here $p$ is the maximum outcome lag and $q$ is the maximum predictor lag. In an ARMA$(p,q)$, $q$ meant the MA order. The model name tells you which meaning applies. Our ARDL convention includes current $x$; an application may restrict $\beta_0=0$ and use only lagged predictors.

Restriction	Result
All $\beta_\ell=0$	AR$(p)$
All $\phi_j=0$	Distributed-lag regression
$q=0$	Current $x$ plus outcome lags
$\beta_0=0$	ARDL with only lagged $x$ terms

There can be several predictors, each with its own lag length. We use one so that timing and interpretation remain visible. Adding five series and twenty lags does not remove the need to justify each information set.

ARDL and Regression with ARMA Disturbances

An ARDL is not simply another name for every regression involving dependent observations. Compare

$$\begin{aligned} \text{ARDL:}\qquad y_t&=\alpha+\phi_1y_{t-1}+\beta_0x_t+\beta_1x_{t-1}+\epsilon_t,\\ \text{Regression with AR(1) disturbances:}\qquad y_t&=\alpha+\beta x_t+u_t,\qquad u_t=\phi_1u_{t-1}+\epsilon_t. \end{aligned}$$

In the second model, serial dependence is specified for $u_t$, the disturbance around the static regression relationship. Substituting its equation gives

$$y_t=(1-\phi_1)\alpha+\phi_1y_{t-1}+\beta x_t-\phi_1\beta x_{t-1}+\epsilon_t.$$

That is a restricted ARDL: the lagged predictor coefficient must equal $-\phi_1\beta$. An unrestricted ARDL estimates the two predictor coefficients separately. More general ARMA disturbances introduce additional structure. Check the software model before translating an xreg coefficient into an ARDL multiplier.

The Stationary Baseline

For the main workflow, we work with stationary variables and a stable dynamic equation. The AR lag polynomial remains

$$\Phi(L)=1-\phi_1L-\cdots-\phi_pL^p.$$

All its roots must lie outside the unit circle for stable conditional outcome dynamics. For ARDL$(1,1)$ this reduces to $|\phi_1|<1$. With stationary $x$ and suitable stationary disturbances, this permits a stationary solution for $y$. If $x$ also responds to past $y$, stability of the joint system must be checked; the outcome coefficient alone is not enough.

Do not feed trending levels into a conventional stationary-regression workflow and interpret the usual output automatically. Modules 1–2 showed why unrelated random walks can look related. Cointegration provides an important alternative for some nonstationary series; §8.9 previews it. It is not a blanket permission to regress one wandering series on another.

8.3 Exogeneity: What Is Left in the Disturbance? — Core

Start with the Economic Timing

Suppose monthly inflation is $y_t$ and the federal funds rate is $x_t$. A positive association could reflect policy responding to inflation, a delayed effect of policy, or both responding to other news. A regression table does not choose among these stories for us.

Calling a column an “independent variable” only identifies where we put it in the regression. It does not establish independence or exogeneity. We need an assumption about what is left out of the equation and how that omitted part relates to the included information.

Let $\mathcal F_{t-1}$ contain the relevant histories of $y$ and $x$ through period $t-1$. A sufficient conditional-mean assumption for our baseline equation is

$$\boxed{E[\epsilon_t\mid\mathcal F_{t-1},x_t]=0.}$$

Read it as: once we account for the history and current $x$, the remaining disturbance has mean zero. This is a statement about the whole conditioning information set, not a collection of separate pairwise independence claims. Under it, the ARDL right-hand side without $\epsilon_t$ really is the specified conditional mean.

For OLS we need valid regressor–disturbance moments. If $z_t$ collects the intercept, included $y$ lags, and included current/lagged $x$ values, this conditional assumption implies $E[z_t\epsilon_t]=0$. Together with appropriate stability, moments, and a nonsingular regressor second-moment matrix, it supports large-sample coefficient consistency. It does not promise exact recovery in a finite sample.

Three Different Relationships

Relationship	What it says	Does it settle the OLS issue?
$x_t$ and $y_t$ are correlated	The variables move together	No; we often include $x$ precisely because they are related.
$\epsilon_t$ is serially uncorrelated	The disturbance does not linearly predict itself across dates	No; it can still be correlated with current $x_t$.
$E[\epsilon_t\mid\mathcal F_{t-1},x_t]=0$	Current regressors and included history do not reveal a remaining conditional-mean component	Yes, it supplies the relevant mean restrictions, with other estimation conditions still required.

A disturbance can be white noise and contemporaneously related to another white-noise shock. “No memory” is not “unrelated to every variable in the model.”

Common Current News and Lagged Feedback

We can extend the opening diagram without creating a same-period cycle:

Common news can connect current x with the outcome disturbance. Lagged feedback is a different route. The diagram depicts possible mechanisms, not identified effects.

Common current shock. Suppose unexpected layoffs move both initial claims and unemployment within a month. If the outcome disturbance contains part of that news, current claims may be related to the disturbance. A coefficient intended to describe a particular structural effect need not be recovered by OLS. Including the observed predictor can still improve a conditional forecast; prediction and structural interpretation are different questions.

Lagged feedback. Suppose

$$x_t=0.5x_{t-1}+0.2y_{t-1}+\epsilon_{x,t}.$$

If $\epsilon_{x,t}$ is independent of the current outcome innovation and both are independent of the past, this dependence on $y_{t-1}$ does not itself invalidate the current outcome equation. The conditioning set already contains that history. Feedback tomorrow is different from correlation with the disturbance today.

Strict exogeneity is stronger. Requiring today’s disturbance to be unrelated to the entire path of $x$, including future values, rules out some feedback that our conditional-mean assumption allows. In the feedback example, today’s outcome innovation helps determine $y_t$, which enters $x_{t+1}$. Future $x$ therefore contains information about today’s disturbance. We do not need that stronger whole-path assumption merely to state the sequential conditional model.

What Would Not Repair the Problem?

• Lagging $x$ automatically. A lag is known earlier, but it may still be related to an omitted persistent driver or a serially dependent disturbance. A lag is a date, not an exogeneity certificate.
• Adding many lags. This may approximate missing dynamics, but it can also exhaust the sample and create collinearity. It does not establish the missing moment restrictions.
• A clean residual ACF. This checks a particular implication for residual dependence; it cannot verify that current regressors are uncorrelated with the structural disturbance.
• Heteroskedasticity/autocorrelation-robust standard errors. These address variance estimation when the coefficient estimator has valid moments and suitable asymptotics. They do not remove inconsistency caused by invalid moments.

A reporting habit. State the goal before interpreting the coefficient: “conditional prediction given this information,” “response along this specified input path,” or “causal response to an intervention.” The last claim requires additional causal assumptions, including a defensible source of variation and a relationship that remains valid under the intervention. The conditional-mean equation alone does not supply them.

Technical Note — Consistency Is Not Exact Finite-Sample Inference

Write $y_t=z_t'\vartheta+\epsilon_t$. The OLS estimation error satisfies

$$\hat\vartheta-\vartheta= \left(T^{-1}\sum_tz_tz_t'\right)^{-1} \left(T^{-1}\sum_tz_t\epsilon_t\right).$$

Under stationary/ergodic sampling, suitable finite moments, $E[z_t\epsilon_t]=0$, and positive definite $E[z_tz_t']$, the second factor converges to zero and the first to a finite matrix. Asymptotic normality and a valid standard error require additional central-limit and variance conditions. Conditional homoskedasticity yields the familiar covariance expression; heteroskedasticity requires an appropriate robust estimate.

Because lagged $y$ contains earlier disturbances, the full sample regressor matrix is not generally independent of the full disturbance vector. Do not transfer the exact finite-sample unbiasedness and textbook exact $t/F$ distributions from a strictly exogenous Gaussian regression without checking their conditions. The Core claim is large-sample justification under stated restrictions, not an exact small-sample guarantee.

8.4 Build the DGP, Then Break Its Assumptions — Core

A Persistent Predictor, Not Five Unrelated Columns

We will generate both series, so we know what the outcome coefficients are. Let

$$\begin{aligned} x_t&=0.5x_{t-1}+c\,y_{t-1}+\epsilon_{x,t},\\ y_t&=0.2+0.6y_{t-1}+0.4x_t+0.2x_{t-1}+\epsilon_t. \end{aligned}$$

The symbol $c$ is the feedback coefficient in this demonstration. The two current shocks have variance one and correlation $r$; shock pairs are independent across dates. We vary $c$ and $r$ while keeping the structural outcome coefficients fixed.

Case	Feedback $c$	Current-shock correlation $r$	Prediction before running
Baseline	0	0	OLS should approach the DGP coefficients as the sample grows.
Shared current shocks	0	0.7	Current $x$ reveals part of the structural disturbance; structural coefficients need not be recovered.
Lagged feedback only	0.2	0	Conditional mean restriction can still hold, despite dependence between future $x$ and today’s innovation.

All three parameter configurations have stable joint dynamics. In the feedback case the two transition eigenvalues are approximately 0.887 and 0.293, inside one in absolute value. This is the system analogue of the characteristic-root condition, equivalent to lag-polynomial roots outside the unit circle.

The shared-shock case deliberately breaks the innovation interpretation: there $\epsilon_t$ remains the structural DGP disturbance, but it is no longer an innovation relative to $\mathcal F_{t-1}$ and $x_t$. Keep that distinction in mind when comparing its known coefficients with regression estimates.

Build the Simulator Inline

Run the code in order. It needs the course packages ggplot2 and forecast; the empirical section also uses urca. Cached CSV files make the chapter runnable without a network connection or an API key. The two functions built here have identical canonical copies in helpers/ardl_utils.R for use in later modules.

library(ggplot2)
theme_set(theme_bw(base_size = 12))
Sys.setenv(OMP_NUM_THREADS = "1", OPENBLAS_NUM_THREADS = "1",
           MKL_NUM_THREADS = "1")
# Run from the project root or from notes/ (Quarto's default R directory).
data_dir <- if (dir.exists("data")) "data" else "../data"
stopifnot(dir.exists(data_dir))

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 pair of equations. We compute $x_t$ first because this DGP specifies it before current $y_t$. We then discard a burn-in period so the arbitrary zero starting values have little influence. Burn-in is useful for a stable system; it cannot turn an unstable system into a stationary one.

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

We construct lags with an explicit missing value at the beginning. This is more transparent than relying on a package’s lag() convention: for a plain numeric vector, stats::lag() changes time attributes rather than necessarily producing the shifted column you intended. lm() drops the first unavailable-lag row here.

One Draw Is Not the Distribution

An estimate can miss the truth through ordinary sampling variation. To see the pattern, repeat the same small experiment. This is a short, serial Monte Carlo: one seed at the beginning, successive random draws, no parallel workers.

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 simulation standard deviation describes variation across samples. The Monte Carlo standard error of the mean describes the numerical precision of this simulation summary. Neither is a correction for estimating the wrong target.

The baseline and feedback-only cases should center near the structural coefficients, with finite-sample error. The shared-shock case should systematically miss some of them. Increasing the sample size makes that discrepancy clearer; it does not make the structural interpretation valid.

Why the Shared-Shock Regression Can Still Look Clean

Because the simulator uses jointly Gaussian shocks with unit variances,

$$E[\epsilon_t\mid\epsilon_{x,t}]=0.7\epsilon_{x,t}.$$

In the shared-shock case, $\epsilon_{x,t}=x_t-0.5x_{t-1}$. Consequently the conditional mean is

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

The structural contemporaneous coefficient is 0.4, but the conditional projection coefficient is $0.4+0.7=1.1$. The structural lag coefficient is 0.2, but its conditional counterpart is $0.2-0.7(0.5)=-0.15$. OLS has not randomly malfunctioned: it is estimating a different conditional relationship than the structural one we wanted to interpret.

After accounting for that projection, the remaining disturbance is $\epsilon_t-0.7\epsilon_{x,t}$. In this constructed example it is independent of current $x$ innovations and serially independent. Residuals can therefore look clean despite failure to recover the structural coefficients.

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

Break it again. Increase the sample size to 2,000 while keeping the shared-shock design. Before running, predict whether $\hat\beta_0$ will move toward 0.4 or toward 1.1. Then set shock_corr = 0 and keep feedback = 0.2. Explain why that change is economically different.

Application callback. A current policy rate or current claims measure may both reflect news about the outcome. An empirical coefficient could be a useful conditional association without measuring the causal response to an isolated policy or labor-market intervention. In §8.7 we make a forecasting claim and use dated predictor lags; we do not relabel its coefficients as structural effects.

8.5 How a Change Travels Through the Equation — Core

Specify the Input Experiment First

A distributed lag and an autoregression both spread a change across time. Start with

$$y_t=0.6y_{t-1}+0.4x_t+0.2x_{t-1}+\epsilon_t.$$

Compare two paths with the same initial conditions and the same outcome shocks. Change only the imposed $x$ path. Subtracting the paths cancels the shocks and gives a deterministic recursion for the difference in $y$.

There are two different experiments:

1. Pulse: raise $x_t$ by one unit for one period; return later $x$ values to their baseline path.
2. Permanent step: raise $x_t,x_{t+1},\ldots$ by one unit in every subsequent period.

These are conditional input experiments. When $x$ itself has persistence or responds to $y$, the imposed path need not be the path generated by a one-unit shock to $x$’s own equation. We will distinguish those experiments below.

Impact and Horizon-Specific Pulse Multipliers

Let $m_h$ be the difference in $y_{t+h}$ caused by the imposed one-period unit pulse, holding all other inputs and innovations fixed. In an ARDL$(1,1)$,

$$\begin{aligned} m_0&=\beta_0,\\ m_1&=\phi_1\beta_0+\beta_1,\\ m_h&=\phi_1m_{h-1},\qquad h\ge2. \end{aligned}$$

The impact multiplier is 0.4. One period later the response is $0.6(0.4)+0.2=0.44$: part comes through inherited $y$, and part through lagged $x$. At horizon two it is $0.6(0.44)=0.264$.

Thus $\beta_1$ alone is not the total response after one period. It describes one route into the equation, not the sum of all paths in the time-unrolled graph.

Cumulative and Permanent-Step Multipliers

The cumulative pulse response through horizon $H$ is $\sum_{h=0}^{H}m_h$. For this linear, time-invariant equation, that sum also equals the response at horizon $H$ to a permanent unit step. Some texts call this an interim multiplier; we state “cumulative through $H$” to keep the experiment explicit.

Horizon	Pulse response $m_h$	Permanent-step response $\sum_{j=0}^{h}m_j$
0	0.4000	0.4000
1	0.4400	0.8400
2	0.2640	1.1040
3	0.1584	1.2624
Long run	0	1.5000

The intercept does not affect these differences because it is the same in both paths. If the outcome is measured in percentage points and the predictor in thousands of claims, the response is in percentage points per additional thousand claims. Units belong in the interpretation.

Long-Run Multiplier

If the conditional outcome dynamics are stable, a permanent unit increase in $x$ eventually changes $y$ by

$$\boxed{\frac{\sum_{\ell=0}^{q}\beta_\ell} {1-\sum_{j=1}^{p}\phi_j}.}$$

For the numerical example it is $(0.4+0.2)/(1-0.6)=1.5$. You can obtain this by summing the pulse responses, or by comparing two steady-state equations whose predictor levels differ by one.

This is a long-run change, not the mean of $y$. With constant predictor level $\bar x$, the steady-state level would instead solve $\bar y=[\alpha+(\sum\beta_\ell)\bar x]/[1-\sum\phi_j]$.

Even when the estimated dynamics are stable, a denominator close to zero makes the estimated long-run multiplier sensitive to small coefficient errors. Report that uncertainty rather than treating a large plug-in ratio as a precisely known response. The optional Technical Note below explains why the coefficient estimates must be considered jointly.

Build and Plot the Responses

The following function implements the recursion for any finite AR and predictor lag vectors. beta[1] is $\beta_0$; phi[1] is $\phi_1$.

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

A second worked example. Keep the direct predictor effects fixed but set $\phi_1=-0.5$. Then $m_0=0.4$, $m_1=-0.5(0.4)+0.2=0$, and all later pulse responses are zero. The long-run step response is $0.6/1.5=0.4$. A positive lag coefficient and negative autoregressive propagation exactly offset at horizon one. You cannot infer the response shape by reading the signs of the $\beta$ coefficients alone.

Break Stability

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

At $\phi_1=1$, the pulse leaves a lasting increment, so repeated increments from a permanent step accumulate. At $\phi_1=1.05$, the denominator in the algebraic ratio is nonzero, but the response does not converge. A nonzero denominator is not a substitute for stability. Reporting $0.6/(1-1.05)=-12$ as the eventual response would contradict the path we just computed.

An Innovation to Persistent $x$ Is a Different Experiment

In the baseline design with independent current shocks and no feedback, let $x_t=0.5x_{t-1}+\epsilon_{x,t}$. An isolated unit innovation to $x$, holding the outcome-shock path fixed, produces the input path $1,0.5,0.25,\ldots$. Its outcome response at horizon one is therefore $0.44+0.4(0.5)=0.64$, not the pulse response 0.44. Later responses combine the ARDL pulse multipliers with the entire innovation-induced $x$ path. With shared current shocks, conditioning on an $x$ innovation can also change the expected outcome disturbance: 0.64 then describes only the mechanical input-path contribution. With feedback, one must follow both equations jointly.

This is why a conditional multiplier is not automatically a structural impulse response. A causal intervention is a further interpretation requiring assumptions that justify changing one mechanism while holding others fixed. A fitted equation by itself does not establish that invariance.

Application callback. If §8.7 models monthly changes in the unemployment rate against growth in claims, a “permanent step in $x$” means permanently higher claims growth, not a once-and-for-all increase in the level of claims. Always translate the multiplier experiment back through the transformations before making an economic statement.

Deeper Dive — The General Recursion

Set $\beta_h=0$ for $h>q$ and $m_h=0$ for $h<0$. Equating the difference between two outcome paths gives

$$m_h=\beta_h+\sum_{j=1}^{p}\phi_jm_{h-j}.$$

The first term is the direct contribution of the original pulse through the distributed lag. The sum is its inherited contribution through past outcome responses. If the responses are absolutely summable, summing both sides over $h\ge0$ yields

$$\sum_{h=0}^{\infty}m_h =\sum_{\ell=0}^{q}\beta_\ell+ \left(\sum_{j=1}^{p}\phi_j\right)\sum_{h=0}^{\infty}m_h,$$

which gives the long-run multiplier. This also shows why simply summing the direct $\beta$ coefficients misses autoregressive propagation.

Technical Note — An Estimated Ratio Has Uncertainty

Replacing coefficients by estimates gives $\hat m_h$ and an estimated long-run multiplier. They are estimates, not known DGP responses. The long-run ratio can be particularly unstable when $1-\sum_j\hat\phi_j$ is close to zero. Uncertainty depends on the joint covariance of the estimated AR and predictor coefficients, not only their separate standard errors. Delta-method or suitable bootstrap calculations are possible under their own regularity conditions; we do not use a plug-in ratio as its own confidence statement.

8.6 Estimation, Diagnostics, and Forecast Timing — Core

Extend the Familiar Workflow

The practical sequence is still specify → estimate → diagnose → revise → evaluate. The new tasks are to justify predictor timing and to keep the goal of the equation explicit.

1. Inspect and transform. Establish units, frequency, seasonal adjustment, release timing, and whether the stationary baseline is appropriate.
2. Choose a modest candidate lag set. Use economic timing and the sampling frequency to bound the search. ACF/PACF patterns of $y$ alone cannot determine all cross-variable lags.
3. Build and inspect the lag columns. Verify that a regressor labelled $x_{t-1}$ actually refers to the previous observation on the intended calendar.
4. Use a common estimation sample. Candidate models with different maximum lags otherwise lose different initial observations. Information criteria must compare the same outcomes under the same likelihood convention.
5. Estimate, then diagnose. A lower BIC does not certify correct timing, exogeneity, or a stable relationship.
6. Evaluate the intended use. For forecasting, score predictions made without future information on the same dates as the benchmark.

A Common-Sample Lag Comparison

We first use the baseline simulated data, where the maintained assumptions are known. The real-data section repeats this workflow with the complications of an observed economic series.

d <- case_data[[1]]
for (j in 1:3) {
  d[[paste0("y_lag", j)]] <- c(rep(NA_real_, j), head(d$y, -j))
  d[[paste0("x_lag", j)]] <- c(rep(NA_real_, j), head(d$x, -j))
}
d_common <- d[complete.cases(d), ]  # All candidates lose the same three rows.
candidate_grid <- expand.grid(p = 1:3, q = 0:3)
candidate_grid$BIC <- NA_real_
candidate_fits <- vector("list", nrow(candidate_grid))
for (i in seq_len(nrow(candidate_grid))) {
  p <- candidate_grid$p[i]
  q <- candidate_grid$q[i]
  rhs <- c(paste0("y_lag", seq_len(p)), "x")
  if (q > 0) rhs <- c(rhs, paste0("x_lag", seq_len(q)))
  fit <- lm(reformulate(rhs, response = "y"), data = d_common)
  candidate_fits[[i]] <- fit
  candidate_grid$BIC[i] <- BIC(fit)
}
knitr::kable(candidate_grid[order(candidate_grid$BIC), ], digits = 2)

	p	q	BIC
4	1	1	1164.42
5	2	1	1169.62
7	1	2	1170.05
10	1	3	1173.59
6	3	1	1175.01
8	2	2	1175.39
11	2	3	1177.34
9	3	2	1180.70
1	1	0	1181.31
12	3	3	1182.00
2	2	0	1185.26
3	3	0	1190.45

selected_fit <- candidate_fits[[which.min(candidate_grid$BIC)]]
knitr::kable(coef(summary(selected_fit)), digits = 3)

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	0.223	0.055	4.070	0
y_lag1	0.556	0.036	15.518	0
x	0.354	0.052	6.823	0
x_lag1	0.274	0.057	4.828	0

In a finite sample, an information criterion may select a model other than the DGP specification. Its penalty trades fit against complexity; it is not an oracle. In BIC(lm(...)), the likelihood parameter count includes the regression coefficients and an estimated disturbance variance, as in Module 4.

Persistent $x$ values can make adjacent lag columns highly correlated. Their individual coefficients may then be imprecise even when the predictor block contributes useful information. Do not drop one lag solely because its individual $t$ statistic is small, without considering the joint lag structure and the intended forecast.

Residual Diagnostics for a Dynamic Regression

Plot the residuals, inspect their ACF, and ask whether the remaining structure points to missing dynamics or instability. A regression-appropriate serial-correlation check is the Breusch–Godfrey auxiliary regression: regress the fitted residual on the original regressors and several of its own lags. If those residual lags add predictive information, the fitted equation has left serial dependence behind.

residual_frame <- data.frame(t = d_common$t, e = residuals(selected_fit))
ggplot(residual_frame, aes(t, e)) + geom_line(color = "#315375") +
  geom_hline(yintercept = 0, linetype = 2) +
  labs(x = "Time", y = "Fitted residual", title = "Does the fitted equation leave structure?")

forecast::ggAcf(residual_frame$e, lag.max = 20) +
  labs(title = "Residual ACF of the selected simulated ARDL")

# Visible Breusch–Godfrey-style auxiliary regression, using four residual lags.
bg_lags <- 4L
bg_frame <- data.frame(e = residuals(selected_fit),
                       model.matrix(selected_fit)[, -1, drop = FALSE])
for (j in seq_len(bg_lags)) {
  bg_frame[[paste0("e_lag", j)]] <-
    c(rep(NA_real_, j), head(bg_frame$e, -j))
}
bg_frame <- bg_frame[complete.cases(bg_frame), ]
bg_fit <- lm(e ~ ., data = bg_frame)
bg_stat <- nobs(bg_fit) * summary(bg_fit)$r.squared
c(LM = bg_stat, df = bg_lags,
  approximate_p = pchisq(bg_stat, df = bg_lags, lower.tail = FALSE))

           LM            df approximate_p 
    4.0107116     4.0000000     0.4045581 

Under the null and the usual stationary regression regularity and variance conditions, the LM statistic has an asymptotic $\chi^2$ reference with degrees of freedom equal to the number of residual lags tested. This is an approximate diagnostic, not an exact small-sample certificate. The first four residual-lag rows were omitted consistently in the auxiliary regression.

Do not mechanically reuse an ARMA-specific Ljung–Box degrees-of-freedom correction for an ARDL. More fundamentally, neither this test nor a residual ACF establishes the relation between the regressors and the structural disturbance. The shared-shock simulation already showed the limit of that reasoning.

Break the specification. Refit the baseline simulated data with lm(y ~ x, data = d_common) and repeat the residual ACF. Predict the omitted paths in the diagram before plotting. If outcome lags belong in the conditional mean, leaving them out can move both coefficient estimates and residual dependence; a standard-error adjustment does not restore the missing dynamics.

The New Predictor Creates a New Forecasting Problem

For a correctly specified ARDL$(1,1)$, a one-step forecast made at time $t$ takes the form

$$\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 new object is $\tilde x_{t+1\mid t}$. If current $x$ enters the outcome equation, a future outcome generally requires a future predictor value. There are several ways to proceed, each answering a different question.

Available information	Interpretation
Only observed predictor lags enter the one-step equation	Forecast from already available history
Current predictor released before the current target outcome	Nowcast using additional within-period information
A stipulated future path for $x$	Conditional/scenario forecast
Future $x$ is unknown and must be predicted	Forecast requiring a predictor model, possibly a joint system

An ARDL with $\beta_0=0$ avoids needing future $x$ for a one-step forecast only if its lagged variables are actually released at the forecast origin. Longer horizons may still require forecasts of $x$.

A Timestamp Is Not a Release Date

Suppose a dataset labels an observation “October.” That may mean it measures October, not that it was available on October 1 or October 31. A forecast must use the information available on the date it is issued.

Reference period:       September              October
Calendar:       ----------------------|---------------------------->
                                      previous-month data released
                                      ^ possible forecast origin
                                                  target-month release later

A forecast of October unemployment made after September data are released can use September variables. It cannot use October information that arrives later simply because all rows are now sitting in one CSV file. Weekly claims also require an explicit aggregation rule before they can become a monthly regressor.

Conditional and Predictor Uncertainty

If the future $x$ path were known, forecast uncertainty would still include future outcome disturbances and estimated-coefficient uncertainty. If that path is itself estimated, an honest unconditional forecast interval must account for predictor uncertainty and any relevant dependence between predictor and outcome forecast errors. Treating a point forecast of $x$ as if it were known generally understates what is uncertain.

For example, even with known coefficients, the one-step error contains $\beta_0(x_{t+1}-\tilde x_{t+1\mid t})+\epsilon_{t+1}$. Its variance includes a predictor-error variance term and, when present, a covariance term. A scenario band conditional on a specified $x$ path omits those terms by design; label it accordingly.

Forecast availability experiment. In the baseline simulation, the following comparison deliberately gives one method future information. It illustrates how easily a favorable score can answer the wrong question.

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 keep estimated coefficients fixed and roll the observed lag information forward. The oracle has an informational advantage that a genuine earlier forecast does not have. Calling its score an ordinary out-of-sample forecasting gain would be leakage. It could instead be a legitimate nowcast score if the predictor really were available first; the release calendar decides.

8.7 Real Data — Do Initial Claims Add Usable Information? — Core

Module 7 ended with a limitation: combining forecasts can reshuffle the information in a pool, but every member of our pool used unemployment’s own history. Initial unemployment-insurance claims offer a plausible additional signal. More people entering the claims system may warn us that labor-market conditions are deteriorating. Whether that signal improves a particular forecast is an empirical question.

Define the variables before fitting anything

Let $U_t$ denote the monthly U.S. unemployment rate, measured in percent, and let $C_t$ denote the mean of the weekly initial-claims observations whose week-ending Saturdays fall in month $t$. Both source series are seasonally adjusted. $U_t$ comes from the household survey; claims are administrative records of initial benefit applications, so these are related labor-market measures with different populations and definitions. See FRED’s UNRATE metadata and FRED’s ICSA metadata.

We model

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

The outcome is a change in percentage points. The predictor is a log change multiplied by 100, approximately the percentage change in the monthly mean of weekly claims. Averaging four or five weekly counts does not create a monthly total. We assign a week crossing a month boundary to the month containing its Saturday, and use that convention consistently.

These transformations give us a stationary working specification to investigate. They do not declare that both original series are necessarily $I(1)$ in every sample, or that differencing automatically solves specification problems. We check the transformed series on training data before using ordinary stationary-regression reasoning.

The timing contract matters. A month’s label is its reference period, not its publication date. BLS usually releases the unemployment rate early in the following month; the BLS definitions distinguish reference, interview, and release dates. The DOL weekly claims release likewise reports a completed week after its Saturday ending date.

For this exercise, the forecast origin is after both the previous month’s unemployment release and the release covering its final claims Saturday have arrived, but before the target month’s unemployment release. Thus, when forecasting $y_t$, we use $y_{t-1}$ and $x_{t-1}$ and their earlier lags. We deliberately ignore partial target-month claims. This is one step ahead in the monthly observation index, and ordinarily a nowcast of the current calendar month; it is not a forecast made at the end of the preceding calendar month. We write its information set as $\mathcal I_t^{\mathrm{release}}$ and its forecast as $\tilde y_{t\mid\mathcal I_t^{\mathrm{release}}}$ to keep that distinction visible.

The cached observations are revised historical data, not releases archived at each forecast origin. The unemployment and claims caches were also collected on different dates. Our exercise therefore isolates temporal fitting and evaluation on a revised-data pseudo-out-of-sample basis. It does not reproduce a real-time forecaster’s historical data vintages or exact release calendar. A real deployment must check the releases actually available at its origin, especially during delays, and use a vintage archive to evaluate real-time performance.

Retrieve once, then work offline

The live retrieval is shown for students with a FRED key. It is not executed when the notes render. The stored claims cache covers 1989–2019; 1989 supplies the pre-sample observations needed for differences and lags. The existing unemployment cache is left unchanged.

library(fredr)
fredr_set_key(Sys.getenv("FRED_KEY"))
claims_live <- fredr(
  series_id = "ICSA",
  observation_start = as.Date("1989-01-01"),
  observation_end = as.Date("2019-12-31")
)
unemployment_live <- fredr(
  series_id = "UNRATE",
  observation_start = as.Date("1988-12-01"),
  observation_end = as.Date("2019-12-31")
)

Load the frozen files and check their calendars before merging. Missing weeks cannot be silently dropped from a monthly mean, and missing months cannot be silently treated as adjacent observations.

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 runs from January 1990 through December 2017. January 2018 through December 2019 is held out, before the pandemic’s extraordinary disruption. We fix the candidate models and selection rule before looking at that holdout. Restricting the illustration to this period limits its scope; it does not establish that the same model would work through 2020 or in every business-cycle phase.

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.

From Module 2, the ADF null is a unit root, while the KPSS null here is stationarity around a constant. Use a drift in the ADF regression, choose its augmentation lag by BIC from a maximum of 12, and use 12 lags in the KPSS long-run variance estimate. These choices are stated in advance rather than varied until a preferred conclusion appears.

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

Reject the ADF null when its statistic is more negative than its critical value; reject the KPSS null when its statistic is larger than its critical value. Here the training-sample evidence is consistent with treating both transformed series as stationary: both ADF tests reject at 5%, and neither KPSS test rejects. These finite-sample tests and plots support a working specification, not a proof of stationarity, exogeneity, constant variance, or stability across future regimes.

Compare a small AR with a lag-only ARDL

The benchmark uses $p$ own lags. The larger information set adds $q$ claims lags:

$$y_t=\alpha+\sum_{j=1}^{p}\phi_j y_{t-j} +\sum_{\ell=1}^{q}\beta_\ell x_{t-\ell}+\epsilon_t.$$

There is no current $x_t$ term. This is the earlier ARDL with the contemporaneous coefficient constrained to zero. It avoids needing a complete target-month claims average at our origin. It does not establish that lagged claims are structurally exogenous: omitted persistent labor-market shocks could still matter.

We consider $p=1,2,3$ and $q=1,2,3$ for the ARDL and $p=1,2,3$ for the AR benchmark. Every candidate uses the same 336 training rows. BIC selects one model within each family; the holdout selects neither lags nor coefficients. In the code, q = 0 is merely the table label for an AR with no claims term, not an ARDL containing current claims.

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

The selected benchmark has $p=$ 3; the selected ARDL has $p=$ 3 and $q=$ 3. AIC() and BIC() for lm count the innovation variance as a parameter, consistent with the course’s information-criterion convention. BIC measures penalized in-sample fit under its likelihood convention; it does not guarantee a win on future observations.

app_coef <- data.frame(
  Term = names(coef(app_ardl)),
  Estimate = unname(coef(app_ardl)),
  Conventional_SE = unname(sqrt(diag(vcov(app_ardl))))
)
knitr::kable(app_coef, digits = 4, row.names = FALSE)

Term	Estimate	Conventional_SE
(Intercept)	-0.0008	0.0076
y_l1	-0.0523	0.0531
y_l2	0.1120	0.0515
y_l3	0.1823	0.0513
x_l1	0.0069	0.0019
x_l2	0.0093	0.0020
x_l3	0.0091	0.0020

For example, x_l1 estimates $\hat\beta_1=$ 0.0069: a one-unit increase in $x_{t-1}$—approximately one percentage point faster claims growth—is associated with that many percentage points of predicted unemployment-rate change, holding the other included lags fixed. It is not an estimate of the effect of forcing people to file claims. The conventional standard errors presume the usual variance conditions and do not account for lag selection; we do not promote their ratios into post-selection causal tests.

Diagnose before celebrating fit

The multiplier formula also needs care after transformation: a permanent step in $x_t$ means permanently faster claims growth, not one extra claimant or a fixed increase in the claims level. The corresponding response is in unemployment-rate changes. We keep the simple level-path multiplier calculations in the simulation and use this observational application for prediction and diagnostics.

Plot the selected model’s training residuals and their ACF. A persistent pattern means the chosen short lag structure has left information behind.

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 following Breusch–Godfrey auxiliary regression tests for residual serial correlation through lag 12. Regress the OLS residuals on the original regressors and twelve residual lags, drop the same first twelve rows for that auxiliary regression, and use its effective sample size times $R^2$. The chi-squared reference is asymptotic and depends on the maintained regression conditions; this is a diagnostic rather than 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 diagnostic $p$-value is 8.48e-06. It rejects at 5%, so this short specification leaves evidence of serial dependence. Treat the forecast comparison below as an evaluation of an imperfect candidate, and revisit the lag range or dynamics using training data before deployment. Heteroskedasticity and structural change remain separate concerns. Improving standard errors would not remove an omitted dynamic term or an endogenous regressor.

One disciplined revision using training data

The small-grid winner is not the end of the workflow. It selects the largest allowed own-lag and claims-lag orders, and its residual diagnostics reject the short specification. Make one substantive expansion: allow up to twelve monthly outcome lags, keeping the claims lag maximum at three. Twelve gives the outcome a full year of memory; it does not assert that the already seasonally adjusted series needs seasonal differencing. Expand the AR benchmark’s own-lag menu too. All candidates still use the same 336 training observations, and BIC remains the selection rule. No forecast error from 2018–2019 enters this revision.

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 revised choices are AR(6) and a lag-only ARDL with $p=$ 6 and $q=$ 3. Its smallest AR lag-polynomial root modulus is 1.141; values above one satisfy the course’s AR stability condition. The revised BG $p$-value is 0.0031. Serial dependence remains detectable at 5%. Stop this illustration here and report that limitation; do not keep changing the model until the diagnostic becomes insignificant. Both the expanded menu and the diagnostic use training data only. Repeated testing and choosing the specification from the same data make these $p$-values exploratory checks, not untouched confirmatory evidence.

The forecast comparison now uses these two revised, training-selected specifications. The earlier coefficient table described the initial short model; coefficients in the revised model can differ because its conditioning set has changed.

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

Roll the information set, keep the coefficients fixed

The first evaluation origin uses observations through December 2017 once those data have been released. As January 2018 is subsequently released, it becomes an observed lag for the February forecast. The coefficients remain those estimated from the training window. Calling predict() on the held-out rows is safe here because the formulas refer only to lag columns: the target month’s y, x, and unemployment level are never regressors.

For either model,

$$\tilde U_{t\mid\mathcal I_t^{\mathrm{release}}} =U_{t-1}+\tilde y_{t\mid\mathcal I_t^{\mathrm{release}}}.$$

Add each predicted change to the latest observed unemployment level. Do not accumulate all 24 predicted changes from December 2017: that would be a different, recursive multi-step exercise. A no-change forecast for $U_t$ corresponds to forecasting $y_t=0$.

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

RMSE and MAE have units of unemployment percentage points. They are identical whether calculated for changes or the reconstructed one-step levels because the same observed $U_{t-1}$ is added to both sides. MAPE is inappropriate for changes that can be zero or negative, so it is omitted.

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 24 dates, the selected ARDL’s RMSE is 0.1057, compared with 0.1129 for the selected AR and 0.1099 for carrying the last unemployment rate forward. The ARDL has the lower RMSE than the AR in this particular window. A short holdout ranking is descriptive evidence about these specifications and dates, not a universal ranking of univariate and multivariate models. We do not use it to retune the lag orders and then report the same dates as untouched test data.

What have we learned about exogeneity? Even a forecast improvement would establish neither a policy effect of claims nor the absence of omitted labor-market shocks. Lagged claims can be a useful signal without being an externally controlled input. Our evidence consists of training diagnostics and forecast errors under a stated information convention. A structural interpretation would need an additional argument about the disturbance and the stability of the relationship under the proposed intervention.

Check Your Understanding

1. Why is averaging weekly claims different from constructing a monthly total? What happens if the final Saturday is missing?
2. Why is the forecast origin later than the last day of month $t-1$? Which data would be unavailable if we moved the origin to that day?
3. Interpret a coefficient of 0.01 on x_l1, stating both units and the conditioning variables. Which causal claim would go beyond it?
4. Why do the change and reconstructed-level RMSEs match? Why would a 24-month recursive forecast be different?
5. If claims improve the training BIC but not holdout RMSE, which result should determine a claim about performance on the holdout? Does either result establish exogeneity?

8.8 Looking Ahead — Granger Causality

The forecast exercise asks whether another variable supplies useful information. Granger causality formalizes a related question: does the past of $x$ improve prediction of $y$ beyond the past information already being used?

For a first linear stationary illustration, use only lagged predictors:

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

The null of no incremental linear predictive content from those $x$ lags is

$$H_0:\beta_1=\cdots=\beta_q=0.$$

It is a joint restriction: one insignificant lag coefficient is not a test of the entire predictor history. Module 9 develops the restricted/unrestricted comparison, its assumptions, and tests in both directions.

The conditioning set matters. Including current $x_t$ and then testing the lag coefficients asks whether past $x$ adds information conditional on current $x$. That is a different information-set question from the usual forecast based on the past. This is why the first Granger example omits the contemporaneous term.

Return to the graph. An arrow from past $x$ to current $y$ can describe predictive precedence. It does not tell us whether $x$ is a direct causal driver, an early indicator of a common driver, or a proxy for omitted information. Results may change when another variable is added, when the sampling frequency changes, or when different lags are used. There may also be predictive feedback in both directions across dates.

A useful forecasting relationship can therefore coexist with no identified intervention effect. Module 9 teaches what the predictive test says and what it leaves unanswered. This preview is optional and is not part of Module 8’s assessable surface.

8.9 Looking Ahead — Cointegration and Error Correction

Return to the Random Walks

Modules 1–2 showed that two unrelated wandering series can produce an impressive levels regression. Differencing helps us avoid treating unrelated trends as a stable relationship. But must every relationship between nonstationary levels be discarded?

Consider Harry and his dog, Mouse. Both wander; the leash keeps their separation from wandering without bound. Karen walks independently. Knowing Mouse’s location may tell us where to look for Harry even though neither location has a fixed mean.

We can generate that distinction with a shared random-walk component and a stationary separation:

set.seed(637611)
n_preview <- 300
harry <- cumsum(rnorm(n_preview))
karen <- cumsum(rnorm(n_preview))
separation <- as.numeric(arima.sim(model = list(ar = 0.6), n = n_preview))
mouse <- harry + separation
wandering <- data.frame(
  t = rep(seq_len(n_preview), 3),
  value = c(harry, mouse, karen),
  series = rep(c("Harry", "Mouse", "Karen"), each = n_preview))
ggplot(wandering, aes(t, value, color = series)) + geom_line() +
  labs(x = "Time", y = "Location", color = NULL,
       title = "Harry and Mouse share a wandering component") +
  theme(legend.position = "bottom")

separation_plot <- rbind(
  data.frame(t = seq_len(n_preview), gap = mouse - harry, pair = "Mouse minus Harry"),
  data.frame(t = seq_len(n_preview), gap = karen - harry, pair = "Karen minus Harry"))
ggplot(separation_plot, aes(t, gap)) + geom_line(color = "#315375") +
  facet_wrap(~ pair, ncol = 1, scales = "free_y") +
  labs(x = "Time", y = "Separation", title = "A shared trend can cancel in a difference")

Here Harry and Mouse are each $I(1)$, but mouse - harry is $I(0)$ by construction. Karen minus Harry is a random walk. Real data do not come with that construction label; a convincing picture alone is insufficient evidence.

For two $I(1)$ series, cointegration means that a nontrivial linear combination is stationary. In the simple normalization,

$$y_t\sim I(1),\quad x_t\sim I(1),\quad y_t-a-bx_t\sim I(0).$$

The stationary combination describes a long-run relation, with deviations around it. It does not mean the series are identical, never separate, or are necessarily linked by an identified causal effect.

The Familiar ARDL Contains an Error-Correction Form

Subtract $y_{t-1}$ from both sides of the ARDL$(1,1)$ and use $x_t=x_{t-1}+\Delta x_t$:

$$\Delta y_t =\alpha+(\phi_1-1)y_{t-1} +(\beta_0+\beta_1)x_{t-1}+\beta_0\Delta x_t+\epsilon_t.$$

If $\phi_1\ne1$, the same equation can be written

$$\Delta y_t =-(1-\phi_1) \left[y_{t-1}-\frac{\alpha}{1-\phi_1} -\frac{\beta_0+\beta_1}{1-\phi_1}x_{t-1}\right] +\beta_0\Delta x_t+\epsilon_t.$$

The last regression term describes short-run change. The bracket is the lagged deviation from a level relation. With $\phi_1<1$, a positive deviation contributes negatively to next period’s change, other inputs held fixed. The familiar long-run multiplier reappears inside that relation.

This is algebra, not a cointegration test. The rewriting is available for stationary variables too, and it does not establish that the bracket is stationary for an arbitrary pair of $I(1)$ series. Cointegration, dynamic adjustment, and admissible conditioning assumptions need their own evidence.

Module 10 develops Engle–Granger testing and error-correction modeling. In particular, a unit-root test applied to estimated cointegrating residuals does not use ordinary ADF critical values without adjustment. We leave the testing workflow, Johansen methods, and ARDL bounds testing outside this introduction. The lesson here is why differencing can remove spurious associations while also throwing away a genuine long-run restriction.

8.10 Looking Ahead — From a Conditional Equation to a System

The ARDL has treated $x$ as information supplied to the outcome equation. For a long-horizon forecast, or when feedback matters, we may need to model $x$ as well. A simple reduced-form VAR gives each variable an equation involving both histories:

$$\begin{aligned} y_t&=a_y+b_{yy}y_{t-1}+b_{yx}x_{t-1}+\epsilon_{y,t},\\ x_t&=a_x+b_{xy}y_{t-1}+b_{xx}x_{t-1}+\epsilon_{x,t}. \end{aligned}$$

The $a$ and $b$ subscripts distinguish the two equations in this preview. Each equation has own-lag and cross-lag arrows; these arrows point forward in time. Current reduced-form innovations may be correlated. That correlation is not drawn as two opposing same-period causal arrows and called a DAG.

A VAR models the variables’ joint evolution and can propagate a forecast through both equations. It does not automatically identify causal shocks or turn a correlation into a policy experiment. Modules 11–12 separate the mechanics of a system from the assumptions needed to interpret its shocks.

ARDL is useful precisely because one equation makes these issues visible. It can be a defensible conditional model for a particular question even when other questions call for a system.

---

Common Pitfalls and Misconceptions

1. “An ARDL is an ARMA with an extra column.” Our basic ARDL has outcome lags and observed predictor lags, with no explicit MA innovation terms. Regression with ARMA disturbances is related but may impose different restrictions.
2. “The diagram has an arrow, so the effect is causal.” The arrow records a posited model relationship. Identification requires evidence and assumptions beyond drawing or estimating it.
3. “An exogenous regressor must be uncorrelated with $y$.” The relevant restriction concerns the regressor’s relation to the disturbance after conditioning on the modeled information.
4. “White-noise errors guarantee exogeneity.” A serially independent structural disturbance can be correlated with a current predictor shock.
5. “Any feedback invalidates the outcome equation.” Lagged feedback can coexist with a valid sequential conditional mean. Common current shocks or simultaneous determination pose a different problem.
6. “Lagging $x$ solves endogeneity.” Persistent omitted information or serially dependent disturbances can still invalidate moments involving lagged regressors.
7. “A clean residual ACF or robust standard errors establish causality.” Neither does. Diagnostics and standard errors have narrower tasks.
8. “$\beta_1$ is the entire next-period response.” The response also inherits the current change through the AR terms.
9. “The cumulative multiplier is just the sum of the $\beta$ coefficients.” That omits autoregressive propagation. Sum the horizon-specific pulse responses.
10. “A nonzero long-run denominator is enough.” The response must converge. Check stability before interpreting the ratio as an eventual response.
11. “A future observed value is allowed because I did not use its outcome.” Future predictor information can leak too. Establish the issue date and release dates.
12. “Adding information must improve holdout performance.” Extra coefficients bring estimation error, and predictive relationships can change. A simple benchmark can win.

Connection to Enders and Other Texts

Enders, Applied Econometric Time Series, 4th edition, Chapter 5 §2, pp. 267–276, develops ADLs and transfer functions. It is the main companion for conditional dynamic responses. Chapter 5 §4, pp. 281–284, provides context for the difficulties of estimating and interpreting structural multivariate relationships. Chapter 5 §8 and Chapter 6 §§1–3 connect to the later Granger-causality and cointegration modules.

Hamilton, Time Series Analysis, §8.2, pp. 215–219, is a more formal reference for dynamic regression and its inferential qualifications. §11.2, pp. 302–309, develops Granger causality and its predictive interpretation. These are optional depth, not a request to memorize asymptotic arguments.

Hyndman, the 2014 Forecasting: Principles and Practice short-course handout, Chapter 10, especially pp. 109 and 113–114, connects regression to forecasting and lagged predictors. Hyndman and Athanasopoulos’s textbook discussion of regression forecasting also emphasizes the need for future predictor values and the distinction between scenario and ex-ante forecasts.

Textbook convention	Course convention and reading advice
ADL versus ARDL	Alternative names; check the maximum-lag convention and whether current $x$ is included.
$a_j$ for AR coefficients in Enders	We use $\phi_j$; $\beta_\ell$ labels predictor coefficients.
Characteristic roots inside the unit circle	We state lag-polynomial roots outside. They are reciprocal descriptions of the same stability condition.
Forecasts often carry hats	From Module 7 onward we use $\tilde y_{t+h\mid t}$; hats remain for estimates and fitted values.
Different uses of “multiplier”	Specify pulse, cumulative pulse, permanent step, or an innovation-induced input path.
Adjustment coefficient often called $\alpha$ in ECM/VECM texts	In this module $\alpha$ remains the intercept; we write the adjustment coefficient explicitly as $\phi_1-1$.

Practice Problems

Core Practice

1. Draw a time-unrolled ARDL$(1,2)$ diagram, including current and two lagged values of $x$, the first lag of $y$, and the current outcome disturbance. Identify which arrows disappear when $\beta_0=0$. Explain why none of these arrows is itself evidence of causal identification.
2. For $y_t=0.3+0.5y_{t-1}+0.6x_t-0.1x_{t-1}+\epsilon_t$, calculate the conditional mean when $y_{t-1}=2$, $x_t=1$, and $x_{t-1}=3$, under the stated zero-conditional-mean assumption. Distinguish this mean from a realization, an estimated fitted value, and a forecast.
3. Explain how $\epsilon_t$ could be white noise but correlated with $x_t$. Use an economic example and say what a residual ACF can and cannot diagnose.
4. Re-run the three-case simulation with 200 and 2,000 observations per sample. Compare average coefficient estimates and their dispersion. Explain which discrepancy more observations can reduce and which they cannot repair.
5. For $\phi_1=0.5$, $\beta_0=0.6$, and $\beta_1=-0.1$, compute pulse responses at horizons 0–3, cumulative responses through those horizons, and the stable long-run step multiplier. Explain each in words.
6. In the previous problem, suppose a unit shock to $x$ produces the path $1,0.8,0.8^2,\ldots$. Calculate the outcome response through horizon two and explain why it differs from an imposed one-period pulse.
7. Compare the lag-selection code using a common sample with a version that lets each model omit only its own missing lags. Explain why differences in their IC rankings would be difficult to interpret.
8. A model uses current monthly claims to predict current monthly unemployment. Describe the release information needed to call its fitted prediction a nowcast. Identify a design that would instead leak future information.
9. In the empirical exercise, choose an alternative lag order using training data only and report its 2018–2019 RMSE and MAE beside the supplied benchmarks. Retain the same forecast origins and dates. Describe the comparison without claiming a universal model ranking or causal effect.
10. A colleague reports a long-run multiplier of $-12$ from $\hat\phi_1=1.05$ and $\hat\beta_0+\hat\beta_1=0.6$. Use a response path to explain the problem. Separately explain why a near-one but stable estimate still warrants care about uncertainty.

Beyond Core (optional)

• Derive the conditional regression coefficients in the shared-shock simulation for a general shock correlation $r$ and predictor persistence $\rho_x$.
• Extend ardl_multipliers()’s examples to two outcome lags and explain the response shape using the general recursion.
• In the Harry–Mouse simulation, replace the stationary separation with another independent random walk. Explain which defining property of cointegration disappears.
• Write both directions of a lag-only predictive model and list the restrictions that would be examined in a bilateral Granger-causality exercise. Testing is developed in Module 9.

Key Takeaways

1. ARDL adds observed information to a dynamic outcome equation. Current and lagged predictors are different from the innovation lags of an MA model.
2. Exogeneity is about the disturbance and the conditioning information. A sequential conditional-mean restriction can allow lagged feedback while excluding correlation with current predictor innovations.
3. Prediction, conditional response, and causal response are different claims. A regression can forecast well and have clean residuals while failing to recover structural DGP coefficients.
4. Dynamic responses follow every route through the equation. The next-period pulse response includes both the direct lag coefficient and inherited outcome persistence. A stable long-run multiplier sums that propagation.
5. A forecast needs an honest calendar. Future predictors, reporting delays, aggregation, and revisions matter alongside model coefficients.
6. The familiar workflow still applies. Use sensible transformations, compare lag specifications on a common sample, diagnose residuals, and evaluate forecasts against simple benchmarks.

Looking Ahead — The Next Questions

Module 9 asks whether the history of one variable adds predictive information beyond the history already being used: Granger causality. It retains the stationary baseline and makes the information set explicit. Module 10 returns to nonstationary levels and asks when a stable relationship survives between them: cointegration and error correction.

The later VAR modules then model the joint evolution of the variables. The thread connecting these topics is already on the opening diagram: which information enters, when does it enter, and what interpretation do the assumptions justify?