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

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

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

}

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

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