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

## ----setup, message = FALSE---------------------------------------------------
library(forecastdom)
library(ggplot2)
data(wg2008)

# WG (2008) Table 1 covers 1872-2005, the entire bundled file.
wg <- wg2008
c(first_year = min(wg$year), last_year = max(wg$year), n = nrow(wg))

## ----helper-------------------------------------------------------------------
recursive_forecasts <- function(y, x, R) {

  n <- length(y)
  P <- n - R
  e_N <- e_A <- f_N <- f_A <- numeric(P)

  for (j in seq_len(P)) {

    idx <- seq_len(R + j - 1)
    f_N[j] <- mean(y[idx])
    fit    <- lm.fit(cbind(1, x[idx]), y[idx])
    f_A[j] <- sum(coef(fit) * c(1, x[R + j]))
    e_N[j] <- y[R + j] - f_N[j]
    e_A[j] <- y[R + j] - f_A[j]

  }

  list(e_N = e_N, e_A = e_A, f_N = f_N, f_A = f_A, year = wg$year[(R + 1):n])

}

## ----main-table---------------------------------------------------------------
specs <- list(
  list(label = "1892+",  R = 20L),
  list(label = "1965+",  R = which(wg$year == 1964)),
  list(label = "1976+",  R = which(wg$year == 1975))
)

run_spec <- function(spec) {

  fc    <- recursive_forecasts(wg$logeqp, wg$log_dp_lag, R = spec$R)
  MSE_N <- mean(fc$e_N ^ 2)
  MSE_A <- mean(fc$e_A ^ 2)
  T_oos <- length(fc$e_N)
  R2    <- 1 - MSE_A / MSE_N
  R2bar <- 1 - (1 - R2) * (T_oos - 1) / (T_oos - 2)
  dRMSE <- sqrt(MSE_N) - sqrt(MSE_A)
  msef  <- mse_f_test(fc$e_N, fc$e_A)
  enc   <- enc_new(fc$e_N, fc$e_A)
  cw    <- cw_test(fc$e_N, fc$e_A, fc$f_N, fc$f_A)

  data.frame(spec      = spec$label,
             T_oos     = T_oos,
             R2bar_pct = 100 * R2bar,
             dRMSE_pct = 100 * dRMSE,
             MSE_F     = unname(msef$statistic),
             ENC_NEW   = unname(enc$statistic),
             CW_stat   = unname(cw$statistic),
             CW_p      = unname(cw$pvalue))

}

tab <- do.call(rbind, lapply(specs, run_spec))
knitr::kable(
  tab, digits = 3, row.names = FALSE, format = "html",
  table.attr = "style='width:auto;'", escape = FALSE,
  caption = "\\(\\bar R^2_{OS}\\) and \\(\\Delta\\mathrm{RMSE}\\) in percent.",
  col.names = c("Spec", "\\(T\\)", "\\(\\bar R^2_{OS}\\)",
                "\\(\\Delta\\mathrm{RMSE}\\)",
                "MSE-F", "ENC-NEW", "CW stat", "\\(p\\)-value"))

## ----wg-fig1, fig.width = 7.5, fig.height = 4.5-------------------------------
# IS residuals from a single regression on the entire sample
fit_full <- lm(logeqp ~ log_dp_lag, data = wg)
is_xy    <- residuals(fit_full)
is_mean  <- wg$logeqp - mean(wg$logeqp)

# OOS residuals starting at year 21
R <- 20L
fc <- recursive_forecasts(wg$logeqp, wg$log_dp_lag, R = R)

is_imp  <- cumsum(is_mean^2) - cumsum(is_xy^2)
oos_imp <- c(rep(NA, R), cumsum(fc$e_N^2) - cumsum(fc$e_A^2))

df <- data.frame(year = wg$year, IS = is_imp, OOS = oos_imp)
df_long <- rbind(
  data.frame(year = df$year, kind = "IS",  value = df$IS),
  data.frame(year = df$year, kind = "OOS", value = df$OOS)
)

ggplot(df_long, aes(x = year, y = value, color = kind,
                    linetype = kind)) +
  annotate("rect",
           xmin = 1973, xmax = 1975,
           ymin = -Inf, ymax = Inf,
           fill = "red", alpha = 0.15) +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_line(linewidth = 0.9, na.rm = TRUE) +
  scale_color_manual(values = c(IS = "black", OOS = "steelblue4")) +
  scale_linetype_manual(values = c(IS = "dashed", OOS = "solid")) +
  labs(x = NULL,
       y = "Cumulative SSE difference (NULL − ALT)",
       title = sprintf("d/p, 1872 - %d", max(wg$year))) +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.position = "top")