---
title: "Replicating Welch & Goyal (2008)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replicating Welch & Goyal (2008)}
  %\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 out-of-sample analysis of Welch and Goyal (2008, *RFS*) for the **log dividend-price ratio** as a predictor of the annual log equity premium. The bundled `wg2008` dataset is built from WG's original `PredictorData.xls` (annual sheet), the data vintage shipped with the published paper. The effective sample is **1872-2005 (134 annual observations)**, matching WG: the file itself begins in 1871, but the one-year lag in the predictor consumes that first row.

The benchmark is the prevailing historical mean (NULL). The alternative is a predictive regression on the lagged predictor (ALTERNATIVE). The predictor is

$$x_{t-1} = \log(D12_{t-1}) - \log(\text{Index}_{t-1}),$$

per WG Section 1, where $D12$ is the 12-month moving sum of dividends and $\text{Index}$ is the S&P 500 price level. (The public plotting script `goyal-welch-plots.R` uses the level D/P, but WG's paper text and Table 1 use the log form.)

WG report five OOS statistics per predictor in Table 1, all of which we compute on the same data:

- $R^2_{OS}$: out-of-sample $R^2$ (Campbell-Thompson 2008).
- $\Delta\text{RMSE}$: $\sqrt{\mathrm{MSE}_N} - \sqrt{\mathrm{MSE}_A}$.
- MSE-F: McCracken (2007) $F$-statistic for equal MSE (`mse_f_test()`).
- ENC-NEW: Clark-McCracken (2001) encompassing test (`enc_new()`).
- CW MSFE-adj: Clark-West (2007) MSFE-adjusted $t$-statistic (`cw_test()`), which WG report in footnote 2.

```{r 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 (WG procedure)

The recursive setup matches WG's `goyal-welch-plots.R` exactly:

- At each year `t`, refit `lm(logeqp ~ log_dp_lag)` on years 1...t-1.
- ALTERNATIVE forecast for year t = fitted value at the
  contemporaneous `log_dp_lag[t]`.
- NULL forecast = mean of `logeqp` over years 1...t-1.

```{r 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])

}
```

## Table 1: five tests across three OOS specifications

WG explore three OOS-start dates: 20 years after the data begins (around 1892), 1965, and the most recent 30 years (1976-2005). The reported $\bar R^2_{OS}$ is the *adjusted* out-of-sample $R^2$ that WG quote in Table 1, applied to the OOS sample of size $T$:

$$\bar R^2_{OS} = 1 - (1 - R^2_{OS})\,\frac{T-1}{T-k},$$

with $k = 2$ parameters (intercept plus predictor).

```{r 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"))
```

Compared to WG's reported d/p numbers:

| Spec  | $T$ | Ours $\bar R^2_{OS}$ | WG $\bar R^2_{OS}$ | Gap  |
|:------|---:|----:|----:|---:|
| 1892+ | 114 | **-2.06** | **-2.06** | 0.00 |
| 1965+ |  41 | -3.73 | -3.69 | 0.04 |
| 1976+ |  30 | -15.22 | -15.14 | 0.09 |

The WG figures come from p. 1474 of the published paper (in-text table for 1892+ and 1976+) and from Table 1 column "Forecasts begin 1965" (1965+).

The longest window matches WG **exactly to two decimals**. The two shorter windows are within 0.1 percentage points. The pattern of the residual gap (zero on the long sample, small on the short samples that are entirely post-1965) is consistent with minor revisions to Goyal's annual data file between the 2007 vintage that fed the published paper and the version currently distributed on Goyal's website. The longest window draws most of its weight from pre-1965 data that has not been revised. The shorter windows are entirely post-1965 and show small drift in proportion to how concentrated they are. Other plausible reconstructions (WG's `log(1 + R - Rfree)` plotting-script formula, hybrid Shiller/CRSP returns, alternative OOS-start boundaries) do not close the gap.

The deterioration of DP through time is unmistakable: a small negative $R^2_{OS}$ in the long sample, more negative in 1965+, and substantially negative for 1976-2005. The McCracken (2004) and Clark-McCracken (2001) asymptotic 5% critical values for $k_2 = 1$ extra regressor are:

| $\pi$ = P/R | MSE-F (5%) | ENC-NEW (5%) |
|---:|---:|---:|
| 0.6  | 1.62 | 2.37 |
| 1.0  | 1.71 | 2.52 |
| 2.0  | 1.82 | 2.70 |

For the 1892+ window ($\pi \approx 6.4$) MSE-F is small and ENC-NEW is below the 5% threshold. For the recent 30 years both statistics are firmly negative or near zero. No window supports a "DP beats the mean" conclusion under WG's tests.

## WG Figure 1: cumulative SSE difference for DP

WG's signature visual is the cumulative squared-error difference

$$\sum_{s \le t}(e_{N,s}^2 - e_{A,s}^2).$$

A rising line means the ALTERNATIVE beats the NULL up to that date, a falling line the opposite. The plot below mirrors the *d/p* panel of WG Figure 1 (IS = dotted, OOS = solid; Oil Shock 1973-1975 shaded in red).

```{r 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")
```

The replicated picture matches the *d/p* panel of WG Figure 1: a quiet first half-century, a climb from WW II to the early 1970s where DP modestly beats the historical mean, a peak around the Oil Shock, and a steep decline through the 1990s as the dividend yield collapsed during the dot-com bull market. The IS line (dashed) sits steadily above zero, so DP looks like a useful in-sample predictor. The OOS line (solid) eventually crashes through zero in the late 1990s. This gap between IS and OOS is what motivated WG's "comprehensive look".

## Takeaway

For the dividend-price ratio at annual frequency, applied with WG's own data and procedure:

- **Long sample (1892+)**: every OOS statistic agrees that DP does not significantly beat the historical mean; the encompassing evidence is weak.
- **1965+ window**: both $R^2_{OS}$ and MSE-F turn clearly negative.
- **Recent 30 years (1976+)**: the DP-augmented forecast is decisively worse than the historical mean.

This is the central message of Welch and Goyal about the dividend-price ratio: in-sample significance does not survive an honest out-of-sample evaluation, and the cumulative-SSE plot makes the structural break around the Oil Shock and the late-1990s decline immediate.

## 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.
- McCracken, M. W. (2007). Asymptotics for out of sample tests of
  Granger causality. *Journal of Econometrics*, 140(2), 719-752.
- Welch, I. and Goyal, A. (2008). A comprehensive look at the
  empirical performance of equity premium prediction. *Review of
  Financial Studies*, 21(4), 1455-1508.