---
title: "Replicating Clark & McCracken (2001)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replicating Clark & McCracken (2001)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 4
)
```

This article reproduces the Phillips-curve forecasting illustration of Clark and McCracken (2001, Section 5) on the bundled `cm2001` dataset. The benchmark forecasts US unemployment from its own lags (univariate AR). The alternative adds lagged inflation (ARX). `enc_new()` tests whether the AR forecasts *encompass* those of the larger ARX model, i.e. whether lagged inflation adds any predictive content beyond what lagged unemployment already contains.

```{r setup, message = FALSE}
library(forecastdom)
data(cm2001)

# CM (2001) sample: 1957:01-1997:08.
cm <- subset(cm2001,
             date >= as.Date("1957-01-01") &
             date <= as.Date("1997-08-01"))

nrow(cm)
```

## Recursive forecasts

Both models are estimated by OLS on an expanding window with `p`
lags of each variable. At each origin we form a one-step-ahead
forecast, then roll forward by one month.

```{r helper}
make_lags <- function(z, p) {

  n <- length(z)
  sapply(seq_len(p), function(k) c(rep(NA, k), z[seq_len(n - k)]))

}

recursive_arx <- function(y, x, p, R) {

  n <- length(y)
  YL <- make_lags(y, p)
  XL <- make_lags(x, p)
  target <- y[(p + 1):n]
  YL <- YL[(p + 1):n, , drop = FALSE]
  XL <- XL[(p + 1):n, , drop = FALSE]
  P <- length(target) - R
  
  e_ar <- e_arx <- numeric(P)
  
  for (j in seq_len(P)) {

    idx     <- seq_len(R + j - 1)
    fit_ar  <- lm.fit(cbind(1, YL[idx, ]),         target[idx])
    fit_arx <- lm.fit(cbind(1, YL[idx, ], XL[idx, ]), target[idx])
    pred_ar  <- sum(coef(fit_ar)  * c(1, YL[R + j, ]))
    pred_arx <- sum(coef(fit_arx) * c(1, YL[R + j, ], XL[R + j, ]))
    e_ar[j]  <- target[R + j] - pred_ar
    e_arx[j] <- target[R + j] - pred_arx

  }
  
  list(e_ar = e_ar, e_arx = e_arx)

}
```

## ENC-NEW across lag orders

CM (2001) consider lag orders `p = 1, 3, 6, 12`. The initial
estimation window is `R = 120` months (10 years); the remaining
observations form the out-of-sample period.

```{r encnew-table}
R <- 120L

rows <- lapply(c(1L, 3L, 6L, 12L), function(p) {

  fc  <- recursive_arx(cm$unrate, cm$infl, p = p, R = R)
  enc <- enc_new(fc$e_ar, fc$e_arx)
  msfe1 <- mean(fc$e_ar ^ 2)
  msfe2 <- mean(fc$e_arx ^ 2)
  data.frame(p        = p,
             n_oos    = length(fc$e_ar),
             pi_ratio = round(length(fc$e_ar) / R, 2),
             MSFE_AR  = msfe1,
             MSFE_ARX = msfe2,
             R2OS_pct = 100 * (1 - msfe2 / msfe1),
             ENC_NEW  = unname(enc$statistic))
})

tab <- do.call(rbind, rows)

knitr::kable(
  tab, digits = 3, row.names = FALSE, format = "html",
  table.attr = "style='width:auto;'", escape = FALSE,
  col.names = c("\\(p\\)", "\\(T_{OOS}\\)", "\\(\\pi = P/R\\)",
                "\\(\\mathrm{MSFE}_{AR}\\)", "\\(\\mathrm{MSFE}_{ARX}\\)",
                "\\(R^2_{OS}\\) (%)", "ENC-NEW"))
```

At every lag order the ARX model achieves a lower MSFE than the pure AR. That is, knowing past inflation reduces out-of-sample unemployment-forecast error. The ENC-NEW values are large relative to the Clark-McCracken (2001) Table 2 asymptotic 5% critical values, which are tabulated in $k_2$ (number of extra regressors, equal to $p$) and $\pi = P/R$. For $k_2 = 6$ and $\pi \approx 2$ the 5% critical value is about 8.6; for $k_2 = 12$ and $\pi \approx 2$ it is about 12.0. Our statistics exceed these benchmarks decisively at every specification, so the AR forecasts do **not** encompass the ARX: lagged inflation carries real predictive content for unemployment.

## Comparison with Clark-West MSFE-adjusted

`enc_new()` answers the encompassing question. `cw_test()` answers the closely related "equal MSFE" question using a $t$-statistic that is approximately $N(0, 1)$ under the null. The two tests typically point the same way, and reporting both is informative.

```{r cw-comparison}
rows2 <- lapply(c(1L, 3L, 6L, 12L), function(p) {

  fc <- recursive_arx(cm$unrate, cm$infl, p = p, R = R)
  f1 <- cm$unrate[(p + 1 + R):nrow(cm)] - fc$e_ar
  f2 <- cm$unrate[(p + 1 + R):nrow(cm)] - fc$e_arx
  cw <- cw_test(fc$e_ar, fc$e_arx, f1, f2)

  data.frame(p         = p,
             CW_stat   = unname(cw$statistic),
             CW_pvalue = unname(cw$pvalue))
             
})

knitr::kable(
  do.call(rbind, rows2), digits = 3, row.names = FALSE,
  col.names = c("$p$", "CW stat", "CW $p$-value"))
```

CW p-values are far below 5% for every lag order, agreeing with the
ENC-NEW decision: the ARX model has significantly lower out-of-sample
MSFE than the pure AR.

## References

- Clark, T. E. and McCracken, M. W. (2001). Tests of equal forecast
  accuracy and encompassing for nested models. *Journal of
  Econometrics*, 105(1), 85-110.
- Clark, T. E. and West, K. D. (2007). Approximately normal tests
  for equal predictive accuracy in nested models. *Journal of
  Econometrics*, 138(1), 291-311.