JCochrane_dog_did_not_bark.Rmd

# The Dog That Did Not Bark: A Defense of Return Predictability (Replication)

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Check packages and install if needed
  check_packages = c("readxl", "xtable", "data.table", "dyn", "kableExtra", "mvtnorm", "ggplot2", "rmarkdown")
  package.check = lapply(
    check_packages, FUN = function(x)
      {
      if (!require(x, character.only = TRUE)) {
        install.packages(x, dependencies = TRUE)
        library(x, character.only = TRUE)
        }
      }
    )

#Load library
  library(readxl)
  library(xtable)
  library(data.table)
  library(dyn)
  library(kableExtra)
  library(mvtnorm)
  library(ggplot2)
  library(rmarkdown)

set.seed(234)
```

## Table 1

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Get relevant data

  #Edit data source here
    data_link = "D:/University of Kent/01. Lectures/SPR - EC843 - Financial Econometrics_done/R/GWAnnual.csv"
    
    GWAnnual = read.csv(data_link)
    data = as.data.table(GWAnnual)
    
    GWAnnual = read.csv(data_link)
    data = as.data.table(GWAnnual)
  
    data = data[, DP := ((1+CRSP_SPvw)/(1+CRSP_SPvwx)-1)] # This means dividends are reinvested
    data = data[, Ret := (1+CRSP_SPvw)]
    data = data[, RF_Ret := (1+Rfree)]
    data = data[,list(Date, Ret, DP, RF_Ret, D12)]
    full_data = ts(data, start = 1926, end = 2017) #data is the whole data set

    
#Regressions
    
  #Estimate Return on Dividend Price Ratio
    OLS1 = dyn$lm(Ret ~ lag(DP, -1), data = full_data)
    OLS1_beta = summary(OLS1)$coef[2, 1]
    OLS1_sd = sd(OLS1[["fitted.values"]])*100
    OLS1_tstat = summary(OLS1)$coef[2, 3]
    OLS1_r2 = summary(OLS1)$r.squared*100
    
  #Estimate Equity Premium on Dividend Price Ratio
    OLS2 = dyn$lm((Ret-RF_Ret) ~ lag(DP, -1), data = full_data)
    OLS2_beta = summary(OLS2)$coef[2, 1]
    OLS2_sd = sd(OLS2[["fitted.values"]])*100
    OLS2_tstat = summary(OLS2)$coef[2, 3]
    OLS2_r2 = summary(OLS2)$r.squared*100
    
  #Estimate Dividend Growth on Dividend Price Ratio
    OLS3 = dyn$lm(D12/lag(D12, -1) ~ lag(DP, -1), data = full_data)
    OLS3_beta = summary(OLS3)$coef[2, 1]
    OLS3_sd = sd(OLS3[["fitted.values"]])*100
    OLS3_tstat = summary(OLS3)$coef[2, 3]
    OLS3_r2 = summary(OLS3)$r.squared*100
    
  #Estimate log return on log DP ratio
    OLS4 = dyn$lm(log(Ret) ~ log(lag(DP, -1)), data = full_data)
    OLS4_beta = summary(OLS4)$coef[2, 1]
    OLS4_sd = sd(OLS4[["fitted.values"]])*100
    OLS4_tstat = summary(OLS4)$coef[2, 3]
    OLS4_r2 = summary(OLS4)$r.squared*100
    OLS4_se = summary(OLS4)$coef[2, 2]
    
  #Estimate log dividend growth on log DP ratio
    OLS5 = dyn$lm(log(D12/lag(D12, -1)) ~ log(lag(DP, -1)), data = full_data)
    OLS5_beta = summary(OLS5)$coef[2, 1]
    OLS5_sd = sd(OLS5[["fitted.values"]])*100
    OLS5_tstat = summary(OLS5)$coef[2, 3]
    OLS5_r2 = summary(OLS5)$r.squared*100
    OLS5_se = summary(OLS5)$coef[2, 2]
    
  #Estimate log DP on its lag
    OLS6 = dyn$lm(log(DP) ~ log(lag(DP, -1)), data = full_data)
    OLS6_beta = summary(OLS6)$coef[2, 1]
    OLS6_sd = sd(OLS6$resid)
```

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Table contents
  Table1 = matrix(
      c(
        OLS1_beta, OLS2_beta, OLS3_beta, OLS4_beta, OLS5_beta,
        OLS1_tstat, OLS2_tstat, OLS3_tstat, OLS4_tstat, OLS5_tstat,
        OLS1_r2, OLS2_r2, OLS3_r2, OLS4_r2, OLS5_r2,
        OLS1_sd, OLS2_sd, OLS3_sd, OLS4_sd, OLS5_sd
      ),
      nrow = 5,
      ncol = 4
    )
  colnames(Table1) = c(
                       "$b$",
                       "$t$",
                       "$R^2(\\%)$",
                       "$\\sigma(bx)(\\%)$"
                     )
  rownames(Table1) = c(
                       "$R_{t+1} = a + b(D_t / P_t ) + \\epsilon_{t+1}$",
                       "$R_{t+1} - R^f_{t+1} = a + b(D_t / P_t ) + \\epsilon_{t+1}$",
                       "$D_{t+1} / D_t = a + b(D_t / P_t ) + \\epsilon_{t+1}$",
                       "$r_{t+1} = a + b(d_t - p_t ) + \\epsilon^r_{t+1}$",
                       "$\\Delta d_{t+1} = a + b(d_t - p_t ) + \\epsilon^{dp}_{t+1}$"
                     )
  
  Table1 %>%    
  kbl(row.names = TRUE,
      caption = "Table 1 - Forecasting regressions",
      escape = FALSE,
      digits = c(2, 2, 1, 1),
      ) %>%
  pack_rows("Regression", 1, 5) %>%
  pack_rows("", 1, 3) %>%
  pack_rows("", 4, 5) %>%
  kable_styling("striped", full_width = FALSE)
```

The result is quite similar to paper's result. The coefficients of dividend growth and log dividend growth are negative as expected, as high dividend/price ratio implies low price, which indicates that the dividend in the future would be low.

## Table 2

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#rho
  rho = exp(mean(log(1/full_data[,"DP"])))/(1+exp(mean(log(1/full_data[,"DP"]))))
  
#Table contents
  Table2 = matrix(
    c(
      OLS4_beta, OLS5_beta, OLS6_beta,
      summary(OLS4)$coef[2, 2], summary(OLS5)$coef[2, 2], summary(OLS6)$coef[2, 2],
      1-rho*OLS6_beta + OLS5_beta, OLS4_beta-(1-rho*OLS6_beta), (1+OLS5_beta-OLS4_beta)/rho,
      sd(residuals(OLS4))*100, cor(residuals(OLS4),residuals(OLS5))*100, cor(residuals(OLS4),residuals(OLS6))*100,
      cor(residuals(OLS4),residuals(OLS5))*100, sd(residuals(OLS5))*100, cor(residuals(OLS5),residuals(OLS6))*100,
      cor(residuals(OLS4),residuals(OLS6))*100, cor(residuals(OLS5),residuals(OLS6))*100, sd(residuals(OLS6))*100,
      0, rho*OLS6_beta-1, OLS6_beta,
      0, rho*0.99-1, 0.99
    ),
    nrow = 3,
    ncol = 8
  )
  colnames(Table2) = c(
                       "$\\hat{b}, \\hat{\\phi}$", "$\\sigma(\\hat{b})$", "$implied$",
                       "$r$", "$\\Delta d$", "$dp$",
                       "$b, \\phi$","$b, \\phi$"
                     )
  rownames(Table2) = c("$r$", "$\\Delta d$", "$dp$")
  
  Table2 %>%    
  kbl(row.names = TRUE,
      caption = "Table 2 - Forecasting regressions and null hypothesis",
      escape = FALSE,
      digits = c(3, 3, 3, 1, 1, 1, 4, 4)) %>%
  add_header_above(c(" " = 1, "Estimates"=3, "$\\epsilon$ s.d. (diagonal) and correlation"=3, "Null 1"=1, "Null 2"=1)) %>%
  kable_styling("striped", full_width = FALSE)
```

## Table 3

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Monte-Carlo Simulation for return & dividend growth

  #Simulation function
    MCSim = function(iteration, phi, rho, data_length, cov_matrix)
    {
      sim_result = matrix(0, nrow = iteration, ncol = 5)
      
        for (i in 1:iteration)
        {
          #generate errors from bi-variate normal
          sim_err = rmvnorm(data_length, c(0, 0), cov_matrix)
          
          #generate dividend price
          sim_dp = rep(0, data_length)
          
            #generate initial observation
            if(phi >= 1) {sim_dp[1] = 0} else {sim_dp[1] = rnorm(1, 0, sqrt(cov_matrix[1, 1]/(1-phi^2)))}
          
            #generate the rest of the observation
            for (j in 2:data_length)
            {
              sim_dp[j] = phi * sim_dp[j-1] + sim_err[j,1]
            }
          sim_dp = ts(sim_dp)
          
          #generate dividend growth
          sim_dg = rep(0, data_length)
          
            #generate initial observation
              sim_dg[1] = rnorm(1, 0, sqrt(cov_matrix[2, 2]/(1-(phi*rho-1)^2)))
              
            #generate the rest of the observation
              for (j in 2:data_length)
              {
                sim_dg[j] = (phi*rho-1) * sim_dp[j-1] + sim_err[j,2]
              }
              
          sim_dg = ts(sim_dg)
          
          #generate return
            sim_ret = sim_err[, 2]-rho*sim_err[, 1]
            sim_ret = ts(sim_ret)
          
          #estimate the regressions using simulated data
            sim_reg_ret = dyn$lm(sim_ret ~ lag(sim_dp, -1))
            sim_reg_dg = dyn$lm(sim_dg ~ lag(sim_dp, -1))
            sim_reg_dp = dyn$lm(sim_dp ~ lag(sim_dp, -1))

          #get results from the regressions:
            #return
              sim_result[i, 1] = summary(sim_reg_ret)$coef[2, 1]
              sim_result[i, 2] = summary(sim_reg_ret)$coef[2, 3]
              
            #dividend growth
              sim_result[i, 3] = summary(sim_reg_dg)$coef[2, 1]
              sim_result[i, 4] = summary(sim_reg_dg)$coef[2, 3]
              
            #dividend price ratio
              sim_result[i, 5] = summary(sim_reg_dp)$coef[2, 1]
        }
      
      sim_result = as.data.frame(sim_result)
      colnames(sim_result) = c("br", "t_br", "bd", "t_bd", "phi")
      return(sim_result)
    }
    
  #Simulation parameters (general)
    iteration = 1000
    data_length = nrow(data)
    
    #Real
      cov_matrix = cov(cbind(OLS6$resid, OLS5$resid))
      
      #Run Simulation 1
        phi = OLS6_beta
        sim_real1 = MCSim(iteration, phi, rho, data_length, cov_matrix)
        
      #Run Simulation 2
        phi = 0.99
        sim_real2 = MCSim(iteration, phi, rho, data_length, cov_matrix)
    
    #Excess
        
      #Log return minus log risk-free return
        OLS7 = dyn$lm(log(Ret)-log(lag(RF_Ret, -1)) ~ log(lag(DP, -1)), data = full_data)
        OLS7_beta = summary(OLS7)$coef[2, 1]
        OLS7_tstat = summary(OLS7)$coef[2, 3]
        OLS7_se = summary(OLS7)$coef[2, 2]
        
      #log dividend growth minus log risk-free return
        OLS8 = dyn$lm(log(D12/lag(D12, -1))-log(lag(RF_Ret, -1)) ~ log(lag(DP, -1)), data = full_data)
        OLS8_beta = summary(OLS8)$coef[2, 1]
        OLS8_tstat = summary(OLS8)$coef[2, 3]
        
      cov_matrix = cov(cbind(OLS6$resid, OLS8$resid))
      phi = OLS6_beta
      
      #Run Simulation
        sim_excess1 = MCSim(iteration, phi, rho, data_length, cov_matrix)
```

```{r, echo = FALSE, warning = FALSE, message = FALSE}
Table3 = matrix(
  c(
    sum(sim_real1$br > OLS4_beta)/iteration*100, sum(sim_excess1$br > OLS7_beta)/iteration*100,
    sum(sim_real1$t_br > OLS4_tstat)/iteration*100, sum(sim_excess1$t_br > OLS7_tstat)/iteration*100,
    sum(sim_real1$bd > OLS5_beta)/iteration*100, sum(sim_excess1$bd > OLS8_beta)/iteration*100,
    sum(sim_real1$t_bd > OLS5_tstat)/iteration*100, sum(sim_excess1$t_bd > OLS8_tstat)/iteration*100
  ),
  nrow = 2,
  ncol = 4
)
colnames(Table3) = c("$b_r$", "$t_r$", "$b_d$","$t_d$")
rownames(Table3) = c("Real", "Excess")

Table3 %>%
kbl(row.names = TRUE,
    caption = "Table 3 - Percent probability values under the $\\phi = 0.939$ null",
    escape = FALSE,
    digits = 1) %>%
kable_styling("striped", full_width = FALSE)
```

The result is the probability that the simulated values are larger than the sample values. Number of iteration: 1000

## Figure 1 - Joint distribution of the coefficients ($b_r$ and $b_d$) and their t-statistics

Number of iteration: 1000. The round blue dot is the sample values, the diamond red dot is the null used to simulate data.

```{r, echo = FALSE, message = FALSE, warning = FALSE}
#Figure 1a
  ggplot(sim_real1, aes(x = br, y = bd)) +
    ggtitle(expression(paste("Figure 1a - Coefficients, ", phi, " = 0.94"))) +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab(expression(b[r])) +
    ylab(expression(b[d])) +
    geom_jitter(size = 0.25) +
    geom_vline(xintercept = OLS4_beta, color = "blue") +
    geom_hline(yintercept = OLS5_beta, color = "blue") +
    geom_point(aes(x = 0, y = rho*OLS6_beta-1),shape = 18, size = 5, color = "red") +
    geom_point(aes(x = OLS4_beta, y = OLS5_beta),shape = 16, size = 5, color = "blue") +
    annotate("text", x = 0.3, y = 0.05, label = paste0(round(sum(sim_real1$br > OLS4_beta & sim_real1$bd > OLS5_beta)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = 0.3, y = -0.25, label = paste0(round(sum(sim_real1$br > OLS4_beta & sim_real1$bd <= OLS5_beta)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -0.1, y = 0.05, label = paste0(round(sum(sim_real1$br <= OLS4_beta & sim_real1$bd > OLS5_beta)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -0.1, y = -0.25, label = paste0(round(sum(sim_real1$br <= OLS4_beta & sim_real1$bd <= OLS5_beta)/iteration*100,1),"%"), color = "red")
  
#Figure 1b
  ggplot(sim_real1, aes(x = t_br, y = t_bd)) +
    ggtitle(expression(paste("Figure 1b - t-stats, ", phi, " = 0.94"))) +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab(expression(paste("t, ", b[r]))) +
    ylab(expression(paste("t, ", b[d]))) +
    geom_jitter(size = 0.25) +
    geom_vline(xintercept = OLS4_tstat, color = "blue") +
    geom_hline(yintercept = OLS5_tstat, color = "blue") +
    geom_point(aes(x = OLS4_tstat, y = OLS5_tstat),shape = 16, size = 5, color = "blue") +
    annotate("text", x = 4, y = 1, label = paste0(round(sum(sim_real1$t_br > OLS4_tstat & sim_real1$t_bd > OLS5_tstat)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = 4, y = -7, label = paste0(round(sum(sim_real1$t_br > OLS4_tstat & sim_real1$t_bd <= OLS5_tstat)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -3, y = 1, label = paste0(round(sum(sim_real1$t_br <= OLS4_tstat & sim_real1$t_bd > OLS5_tstat)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -3, y = -7, label = paste0(round(sum(sim_real1$t_br <= OLS4_tstat & sim_real1$t_bd <= OLS5_tstat)/iteration*100,1),"%"), color = "red")
```

```{r, echo = FALSE, message = FALSE, warning = FALSE}
#Figure 1c
  ggplot(sim_real2, aes(x = br, y = bd)) +
    ggtitle(expression(paste("Figure 1c - Coefficients, ", phi, " = 0.99"))) +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab(expression(b[r])) +
    ylab(expression(b[d])) +
    geom_jitter(size = 0.25) +
    geom_vline(xintercept = OLS4_beta, color = "blue") +
    geom_hline(yintercept = OLS5_beta, color = "blue") +
    geom_point(aes(x = 0, y = rho*0.99-1),shape = 18, size = 5, color = "red") +
    geom_point(aes(x = OLS4_beta, y = OLS5_beta),shape = 16, size = 5, color = "blue") +
    annotate("text", x = 0.3, y = 0.025, label = paste0(round(sum(sim_real2$br > OLS4_beta & sim_real2$bd > OLS5_beta)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = 0.3, y = -0.15, label = paste0(round(sum(sim_real2$br > OLS4_beta & sim_real2$bd <= OLS5_beta)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -0.05, y = 0.025, label = paste0(round(sum(sim_real2$br <= OLS4_beta & sim_real2$bd > OLS5_beta)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -0.05, y = -0.15, label = paste0(round(sum(sim_real2$br <= OLS4_beta & sim_real2$bd <= OLS5_beta)/iteration*100,1),"%"), color = "red")
    
#Figure 1d
  ggplot(sim_real2, aes(x = t_br, y = t_bd)) +
    ggtitle(expression(paste("Figure 1d - t-stats, ", phi, " = 0.99"))) +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab(expression(paste("t, ", b[r]))) +
    ylab(expression(paste("t, ", b[d]))) +
    geom_jitter(size = 0.25) +
    geom_vline(xintercept = OLS4_tstat, color = "blue") +
    geom_hline(yintercept = OLS5_tstat, color = "blue") +
    geom_point(aes(x = OLS4_tstat, y = OLS5_tstat),shape = 16, size = 5, color = "blue") +
    annotate("text", x = 4, y = 0, label = paste0(round(sum(sim_real2$t_br > OLS4_tstat & sim_real2$t_bd > OLS5_tstat)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = 4, y = -6, label = paste0(round(sum(sim_real2$t_br > OLS4_tstat & sim_real2$t_bd <= OLS5_tstat)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -2, y = 0, label = paste0(round(sum(sim_real2$t_br <= OLS4_tstat & sim_real2$t_bd > OLS5_tstat)/iteration*100,1),"%"), color = "red") +
    annotate("text", x = -2, y = -6, label = paste0(round(sum(sim_real2$t_br <= OLS4_tstat & sim_real2$t_bd <= OLS5_tstat)/iteration*100,1),"%"), color = "red")
```

## Table 4

```{r, echo = FALSE}
Table4 = matrix(
  c(
    OLS4_beta/(1-rho*OLS6_beta), OLS5_beta/(1-rho*OLS6_beta), OLS7_beta/(1-rho*OLS6_beta),
    OLS4_se/(1-OLS6_beta^2), OLS5_se/(1-OLS6_beta^2), OLS7_se/(1-OLS6_beta^2),
    OLS4_beta/(1-rho*OLS6_beta)/(OLS4_se/(1-OLS6_beta^2)),
    (OLS5_beta/(1-rho*OLS6_beta)+1)/(OLS5_se/(1-OLS6_beta^2)),
    OLS7_beta/(1-rho*OLS6_beta)/(OLS7_se/(1-OLS6_beta^2)),
    2*pt(-abs(OLS4_beta/(1-rho*OLS6_beta)/(OLS4_se/(1-OLS6_beta^2))), df = data_length-1)*100,
    2*pt(-abs((OLS5_beta/(1-rho*OLS6_beta)+1)/(OLS5_se/(1-OLS6_beta^2))), df = data_length-1)*100,
    2*pt(-abs(OLS7_beta/(1-rho*OLS6_beta)/(OLS7_se/(1-OLS6_beta^2))), df = data_length-1)*100
  ),
  nrow = 3,
  ncol = 4
)
colnames(Table4) = c("$\\hat{b^{lr}}$", "$s. e.$", "$t$", paste("$\\%$", "$p$", "$value$"))
rownames(Table4) = c("$r$", "$\\Delta d$", "Excess $r$")

Table4 %>%
kbl(row.names = TRUE,
    caption = "Table 4 - Long-run regression coefficients",
    escape = FALSE,
    digits = 2) %>%
kable_styling("striped", full_width = FALSE)
```

The p-values are retrieved from the simulation with known $b_r = 0$ and $\phi = 0.940$. I did not change the selection of the coefficients to reduce calculation time, so the p-value is not a range.

## Figure 2 - Distribution of the $b_r^{lr}$

Number of iteration: 1000. The vertical line is the sample value.

```{r, echo = FALSE, message = FALSE}
#Figure 2a
  #phi = 0.94
    ggplot(sim_real1, aes(x = br/(1-rho*OLS6_beta))) +
      geom_histogram(binwidth = 0.1, color = "black", fill = "white") +
      ggtitle(expression(paste(phi, " = 0.94"))) +
      theme(plot.title = element_text(hjust = 0.5)) +
      xlab(expression(paste(b[r], "/(1-", rho, phi, ")"))) +
      ylab("") +
      geom_vline(xintercept = OLS4_beta/(1-rho*OLS6_beta), color = "red") +
      annotate("text", x = 1.25, y = 500, label = "Data", color = "red")
```

```{r, echo = FALSE, message = FALSE}
#Figure 2b
  #phi = 0.99
    ggplot(sim_real2, aes(x = br/(1-rho*0.99))) +
      geom_histogram(binwidth = 0.1, color = "black", fill = "white") +
      ggtitle(expression(paste(phi, " = 0.99"))) +
      theme(plot.title = element_text(hjust = 0.5)) +
      xlab(expression(paste(b[r], "/(1-", rho, phi, ")"))) +
      ylab("") +
      geom_vline(xintercept = OLS4_beta/(1-rho*0.99), color = "red") +
      annotate("text", x = 2.5, y = 300, label = "Data", color = "red")
```

## Figure 3 - Joint distribution of $b_r$ and $\phi$

Number of iteration: 1000. The dashed line ($b_d$) is $b_r = 1- \rho\phi + \hat{b_d}$, the solid line $b_r^{lr}$ is $b_r/(1-\rho\phi) = \hat{b_r}/(1-\rho\hat{\phi})$.

```{r, echo = FALSE, warning = FALSE}
brlr = OLS4_beta/(1-rho*OLS6_beta)
#Figure 3a
  ggplot(sim_real1, aes(x = br, y = phi)) +
    ggtitle(expression(paste(b[r], " and ", phi, ", ", phi, " = 0.94"))) +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab(expression(b[r])) +
    ylab(expression(phi)) +
    geom_jitter(size = 0.25) +
    geom_vline(xintercept = OLS4_beta, color = "blue") +
    geom_hline(yintercept = OLS6_beta, color = "blue") +
    geom_line(aes(x=1-rho*sim_real1$phi+OLS5_beta, y=sim_real1$phi), color = "red", linetype = "dashed") +
    geom_line(aes(x=brlr*(1-rho*sim_real1$phi), y=sim_real1$phi), color = "blue") +
    geom_point(aes(x = 0, y = OLS6_beta),shape = 18, size = 5, color = "red") +
    geom_point(aes(x = OLS4_beta, y = OLS6_beta),shape = 16, size = 5, color = "blue") +
    annotate("text", x = 0.4, y = 1, label = paste0(round(sum(sim_real1$br > OLS4_beta & sim_real1$phi > OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = 0.4, y = 0.5, label = paste0(round(sum(sim_real1$br > OLS4_beta & sim_real1$phi <= OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = -0.1, y = 1, label = paste0(round(sum(sim_real1$br <= OLS4_beta & sim_real1$phi > OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = -0.1, y = 0.5, label = paste0(round(sum(sim_real1$br <= OLS4_beta & sim_real1$phi <= OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = 0.35, y = 0.7, label = expression(b[r]^{lr}), color = "blue") +
    annotate("text", x = 0.3, y = 0.65, label = expression(b[d]), color = "red")
```

```{r, echo = FALSE, warning = FALSE}
brlr = OLS4_beta/(1-rho*OLS6_beta)
#Figure 3b
  ggplot(sim_real2, aes(x = br, y = phi)) +
    ggtitle(expression(paste(b[r], " and ", phi, ", ", phi, " = 0.99"))) +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab(expression(b[r])) +
    ylab(expression(phi)) +
    geom_jitter(size = 0.25) +
    geom_vline(xintercept = OLS4_beta, color = "blue") +
    geom_hline(yintercept = OLS6_beta, color = "blue") +
    geom_point(aes(x = 0, y = 0.99),shape = 18, size = 5, color = "red") +
    geom_point(aes(x = OLS4_beta, y = OLS6_beta),shape = 16, size = 5, color = "blue") +
    geom_line(aes(x=1-rho*sim_real2$phi+OLS5_beta, y=sim_real2$phi), color = "red", linetype = "dashed") +
    geom_line(aes(x=brlr*(1-rho*sim_real2$phi), y=sim_real2$phi), color = "blue") +
    annotate("text", x = 0.4, y = 1, label = paste0(round(sum(sim_real2$br > OLS4_beta & sim_real2$phi > OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = 0.4, y = 0.5, label = paste0(round(sum(sim_real2$br > OLS4_beta & sim_real2$phi <= OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = -0.15, y = 1, label = paste0(round(sum(sim_real2$br <= OLS4_beta & sim_real2$phi > OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = -0.15, y = 0.5, label = paste0(round(sum(sim_real2$br <= OLS4_beta & sim_real2$phi <= OLS6_beta)/iteration*100,1),"%")) +
    annotate("text", x = 0.4, y = 0.65, label = expression(b[r]^{lr}), color = "blue") +
    annotate("text", x = 0.3, y = 0.65, label = expression(b[d]), color = "red")
```

## Table 5

```{r, echo = FALSE, message = FALSE, warning = FALSE}
#Simulation for table 5
  iteration = 1000
  data_length = nrow(data)
  
    #Real
      cov_matrix = cov(cbind(OLS6$resid, OLS5$resid))
      #Run Simulation 3
        phi = 0.9
        sim_real3 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 4
        phi = 0.96
        sim_real4 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 5
        phi = 0.98
        sim_real5 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 6
        phi = 1
        sim_real6 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 7
        phi = 1.01
        sim_real7 = MCSim(iteration, phi, rho, data_length, cov_matrix)
        
    #Excess
      cov_matrix = cov(cbind(OLS6$resid, OLS8$resid))
      phi = OLS6_beta
      #Run Simulation 2
        phi = 0.9
        sim_excess2 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 3
        phi = 0.96
        sim_excess3 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 4
        phi = 0.98
        sim_excess4 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 5
        phi = 0.99
        sim_excess5 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 6
        phi = 1
        sim_excess6 = MCSim(iteration, phi, rho, data_length, cov_matrix)
      #Run Simulation 7
        phi = 1.01
        sim_excess7 = MCSim(iteration, phi, rho, data_length, cov_matrix)
```

```{r, echo = FALSE, message = FALSE, warning = FALSE}
Table5 = matrix(
  c(
    round(sum(sim_real3$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real1$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real4$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real5$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real2$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real6$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real7$br > OLS4_beta)/iteration*100,1),
    round(sum(sim_real3$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_real1$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_real4$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_real5$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_real2$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_real6$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_real7$bd > OLS5_beta)/iteration*100,1),
    round(sum(sim_excess2$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess1$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess3$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess4$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess5$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess6$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess7$br > OLS7_beta)/iteration*100,1),
    round(sum(sim_excess2$bd > OLS8_beta)/iteration*100,1),
    round(sum(sim_excess1$bd > OLS8_beta)/iteration*100,1),
    round(sum(sim_excess3$bd > OLS8_beta)/iteration*100,1),
    round(sum(sim_excess4$bd > OLS8_beta)/iteration*100,1),
    round(sum(sim_excess5$bd > OLS8_beta)/iteration*100,1),
    round(sum(sim_excess6$bd > OLS8_beta)/iteration*100,1),
    round(sum(sim_excess7$bd > OLS8_beta)/iteration*100,1),
    round(OLS6_sd/sqrt(1-0.9^2),2),
    round(OLS6_sd/sqrt(1-OLS6_beta^2),2),
    round(OLS6_sd/sqrt(1-0.96^2),2),
    round(OLS6_sd/sqrt(1-0.98^2),2),
    round(OLS6_sd/sqrt(1-0.99^2),2),
    "$\\infty$",
    "$\\infty$",
    round(log(0.5, base = 0.9),1),
    round(log(0.5, base = OLS6_beta),1),
    round(log(0.5, base = 0.96),1),
    round(log(0.5, base = 0.98),1),
    round(log(0.5, base = 0.99),1),
    "$\\infty$",
    "$\\infty$"
  ),
  nrow = 7,
  ncol = 6
)

colnames(Table5) = c("$b_r$", "$b_d$", "$b_r$", "$b_d$", "$\\sigma (dp)$", "1/2 life")
rownames(Table5) = c("0.90", round(OLS6_beta, 3), "0.96", "0.98", "0.99", "1", "1.01")

Table5 %>%
  kbl(row.names = TRUE,
      caption = paste("Table 5 - The effects of dividend-yield autocorrelation", "$\\phi$"),
      escape = FALSE) %>%
  add_header_above(c("Null $\\phi$"=1, "Real returns"=2, "Excess returns"=2, "Statistics"=2)) %>%
  add_header_above(c(" "=1, "Percentage probability values"=4, "Other"=2)) %>%
  kable_styling("striped", full_width = FALSE)
```

The values in the table are the probability that the coefficients are larger than the sample values, with the number of iteration is 1000.

## Figure 5 - Joint distribution of $b_r$ and $\phi$ estimates, together with regions implied by long-run regressions

```{r, echo = FALSE, warning = FALSE}
#Figure 5
br5 = (1 + (rho*OLS6_beta)^1 + (rho*OLS6_beta)^2 + (rho*OLS6_beta)^3 + (rho*OLS6_beta)^4)*OLS4_beta
br10 = (1 + (rho*OLS6_beta)^1 + (rho*OLS6_beta)^2 + (rho*OLS6_beta)^3 + (rho*OLS6_beta)^4 + (rho*OLS6_beta)^5 + (rho*OLS6_beta)^6 + (rho*OLS6_beta)^7 + (rho*OLS6_beta)^8 + (rho*OLS6_beta)^9)*OLS4_beta
br20 = (1 + (rho*OLS6_beta)^1 + (rho*OLS6_beta)^2 + (rho*OLS6_beta)^3 + (rho*OLS6_beta)^4 + (rho*OLS6_beta)^5 + (rho*OLS6_beta)^6 + (rho*OLS6_beta)^7 + (rho*OLS6_beta)^8 + (rho*OLS6_beta)^9 + (rho*OLS6_beta)^10 + (rho*OLS6_beta)^11 + (rho*OLS6_beta)^12 + (rho*OLS6_beta)^13 + (rho*OLS6_beta)^14 + (rho*OLS6_beta)^15 + (rho*OLS6_beta)^16 + (rho*OLS6_beta)^17 + (rho*OLS6_beta)^18 + (rho*OLS6_beta)^19)*OLS4_beta
brlr = OLS4_beta/(1-rho*OLS6_beta)
brlr_unweighted = OLS4_beta/(1-OLS6_beta)

ggplot(sim_real1, aes(x = br, y = phi)) +
  ggtitle(expression(paste(b[r], " and ", phi, ", with long run regressions"))) +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab(expression(b[r])) +
  ylab(expression(phi)) +
  geom_jitter(size = 0.25) +
  geom_line(aes(x=OLS4_beta, y=sim_real1$phi), color = "blue") +
  geom_line(aes(x=br5/(1 + (rho*sim_real1$phi)^1 + (rho*sim_real1$phi)^2 + (rho*sim_real1$phi)^3 + (rho*sim_real1$phi)^4), y=sim_real1$phi), color = "blue") +
  geom_line(aes(x=br10/(1 + (rho*sim_real1$phi)^1 + (rho*sim_real1$phi)^2 + (rho*sim_real1$phi)^3 + (rho*sim_real1$phi)^4 + (rho*sim_real1$phi)^5 + (rho*sim_real1$phi)^6 + (rho*sim_real1$phi)^7 + (rho*sim_real1$phi)^8 + (rho*sim_real1$phi)^9), y=sim_real1$phi), color = "blue") +
  geom_line(aes(x=br20/(1 + (rho*sim_real1$phi)^1 + (rho*sim_real1$phi)^2 + (rho*sim_real1$phi)^3 + (rho*sim_real1$phi)^4 + (rho*sim_real1$phi)^5 + (rho*sim_real1$phi)^6 + (rho*sim_real1$phi)^7 + (rho*sim_real1$phi)^8 + (rho*sim_real1$phi)^9 + (rho*sim_real1$phi)^10 + (rho*sim_real1$phi)^11 + (rho*sim_real1$phi)^12 + (rho*sim_real1$phi)^13 + (rho*sim_real1$phi)^14 + (rho*sim_real1$phi)^15 + (rho*sim_real1$phi)^16 + (rho*sim_real1$phi)^17 + (rho*sim_real1$phi)^18 + (rho*sim_real1$phi)^19), y=sim_real1$phi), color = "blue") +
  geom_line(aes(x=brlr*(1-rho*sim_real1$phi), y=sim_real1$phi), color = "blue") +
  geom_line(aes(x=brlr_unweighted*(1-sim_real1$phi), y=sim_real1$phi), color = "blue", linetype = "dashed") +  
  geom_point(aes(x = 0, y = OLS6_beta),shape = 18, size = 5, color = "red") +
  geom_point(aes(x = OLS4_beta, y = OLS6_beta),shape = 16, size = 5, color = "blue") +
  annotate(geom="label", x = OLS4_beta, y = 0.6, label = "1", color = "blue", fill = "white") +
  annotate(geom="label", x = 0.165, y = 0.6, label = "5", color = "blue", fill = "white") +
  annotate(geom="label", x = 0.25, y = 0.6, label = "10", color = "blue", fill = "white") +
  annotate(geom="label", x = 0.345, y = 0.6, label = "20", color = "blue", fill = "white") +
  annotate(geom="label", x = 0.405, y = 0.6, label = expression(infinity), color = "blue", fill = "white") +
  annotate(geom="label", x = 0.45, y = 0.7, label = expression(paste(infinity, ", unweighted")), color = "blue", fill = "white")
```

## Figure 6 - Distribution of the Goyal - Welch statistic

Number of iteration: 5000. The vertical line is the sample value.

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Data's ΔRMSE
  FR = matrix(NA, data_length, 2)

  for (i in 20:data_length)
  {
    ret = ts(log(data[["Ret"]]), start = 1926, end = 1926 + i - 1)
    dp = ts(log(data[["DP"]]), start = 1926, end = 1926 + i - 1)
    
    #Forecast residual in the next period using mean:
    FR[i-19, 1] = log(data[["Ret"]][i+1]) - mean(ret)
    
    #Forecast residual in the next period using Prediction Model
    reg = dyn$lm(ret ~ lag(dp, -1))
    FR[i-19, 2] = log(data[["Ret"]][i+1]) - (reg$coef[1]+reg$coef[2]*log(data[["DP"]][i]))
  }
  
  delta_RMSE_sample = sqrt(mean(FR[,1]^2, na.rm = TRUE)) - sqrt(mean(FR[,2]^2, na.rm = TRUE))
  

#Monte-Carlo Simulation
  
  #Simulation function
    MCSim_GW = function(iteration, phi, rho, data_length, cov_matrix)
    {
      delta_RMSE_sim = matrix(NA, nrow = iteration, ncol = 1)
      
        for (i in 1:iteration)
        {
          #generate errors from bi-variate normal
          sim_err = rmvnorm(data_length, c(0, 0), cov_matrix)
          
          #generate dividend price
          sim_dp = rep(0, data_length)
          
            #generate initial observation
            sim_dp[1] = rnorm(1, 0, sqrt(cov_matrix[1, 1]/(1-phi^2)))
            
            #generate the rest of the observation
            for (j in 2:data_length)
            {
              sim_dp[j] = phi * sim_dp[j-1] + sim_err[j,1]
            }
          
          #generate dividend growth
          sim_dg = sim_err[, 2]
          
          #generate return
          sim_ret = rep(0, data_length)
          
            #generate initial observation
            sim_ret[1] = rnorm(1, 0, sd(OLS4$resid))
            
            #generate the rest of the observation
              for (j in 2:data_length)
              {
                sim_ret[j] = (phi*rho-1) * sim_dp[j-1] + (sim_err[j, 1] - rho*sim_err[j, 2])
              }

          #Out of sample RMSE
          FR_sim = matrix(NA, data_length, 2)

          for (k in 20:data_length)
          {
            ret = ts(sim_ret, start = 1, end = k)
            dp = ts(sim_dp, start = 1, end = k)
            
            #Forecast residual in the next period using mean:
            FR_sim[k-19, 1] = sim_ret[k+1] - mean(ret)
            
            #Forecast residual in the next period using Prediction Model
            reg = dyn$lm(ret ~ lag(dp, -1))
            FR_sim[k-19, 2] = sim_ret[k+1] - (reg$coef[1]+reg$coef[2]*sim_dp[k])
          }
  
          delta_RMSE_sim[i] = sqrt(mean(FR_sim[,1]^2, na.rm = TRUE)) - sqrt(mean(FR_sim[,2]^2, na.rm = TRUE))
        }
      
      delta_RMSE_sim = as.data.frame(delta_RMSE_sim)
      colnames(delta_RMSE_sim) = "delta_RMSE_sim"
      return(delta_RMSE_sim)
    }
    
  #Simulation parameters (general)
    iteration = 5000
    data_length = nrow(data)
    cov_matrix = cov(cbind(OLS6$resid, OLS5$resid))
    
      #Run Simulation 1
        phi = OLS6_beta
        simGW_real1 = MCSim_GW(iteration, phi, rho, data_length, cov_matrix)
        
      #Run Simulation 2
        phi = 0.99
        simGW_real2 = MCSim_GW(iteration, phi, rho, data_length, cov_matrix)
```

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Figure 6a
  #phi = 0.94
    ggplot(simGW_real1, aes(x = delta_RMSE_sim)) +
      geom_histogram(color = "black", fill = "white") +
      ggtitle(expression(paste(phi, " = 0.94"))) +
      theme(plot.title = element_text(hjust = 0.5)) +
      xlab(expression(paste(Delta, " rmse"))) +
      ylab("") +
      geom_vline(xintercept = delta_RMSE_sample, color = "red") +
      annotate("text", x = -0.0075, y = 700, label = "Data", color = "red") +
      annotate("text", x = -0.015, y = 200, label = paste(round(sum(simGW_real1 < delta_RMSE_sample)/iteration*100, 0), "%"), color = "red")
```

```{r, echo = FALSE, warning = FALSE, message = FALSE}
#Figure 6b
  #phi = 0.99
    ggplot(simGW_real2, aes(x = delta_RMSE_sim)) +
      geom_histogram(color = "black", fill = "white") +
      ggtitle(expression(paste(phi, " = 0.99"))) +
      theme(plot.title = element_text(hjust = 0.5)) +
      xlab(expression(paste(Delta, " rmse"))) +
      ylab("") +
      geom_vline(xintercept = delta_RMSE_sample, color = "red") +
      annotate("text", x = -0.0075, y = 700, label = "Data", color = "red") +
      annotate("text", x = -0.015, y = 200, label = paste(round(sum(simGW_real2 < delta_RMSE_sample)/iteration*100, 0), "%"), color = "red")
```