### Section 4 - Ludvigson & Ng
## Appendix results: predictors are return-forecasting factors
## Table B.2

rm(list=ls())
source("R/robust_fns.r")
source("R/var_fns.r")
library(xtable)
set.seed(616)

## lndata <- getLNdata()
## estSize <- TRUE

lndata <- getLNdata(beg=198501, end=201612)  ## later sample period, our data
estSize <- FALSE

macro.vars <-  paste0("F", 1:8)
PCs <- paste0("PC", 1:3)
M <- 5000 # number of bootstrap replications

## construct return-forecasting factors
modLN <- lm(xr.avg ~ F1 + F2 + F3 + F4 + F5 + F6 + F7 + F8 + I(F1^3), data=lndata)
lndata$H8 <- c(modLN$fitted.values, rep(NA,12))
modCP <- lm(xr.avg ~ f1 + f2 + f3 + f4 + f5, data=lndata)
lndata$CP <- c(modCP$fitted.values, rep(NA,12))

dgp <- getBootDGP(PCs, macro.vars, lndata)

fmla1 <- get(depvar) ~ CP
fmla2 <- get(depvar) ~ CP + H8
regvars <- c("CP", "H8")
depvars <- c("xr2", "xr3", "xr4", "xr5")
tblCoef <- array(NA, c(10, 2, length(depvars)))
tblR2 <- array(NA, c(4, 6, length(depvars)))
for (i in seq_along(depvars)) {
    depvar <- depvars[i]
    cat("# Dependent variable:", depvar, "\n")
    rval2 <- bootstrapTest(fmla1, fmla2, lndata, dgp, depvar=depvar, adjR2=TRUE, M=M, estSize=estSize, estIMSize=FALSE)
    tblCoef[1:7, , i] <- rval2$tblCoef[,1:2]
    tblCoef[8:9, , i] <- rval2$tblSize[1:2, 1:2]
    tblR2[, , i] <- rval2$tblR2
    ## Power
    rval3 <- bootstrapTest(fmla1, fmla2, lndata, dgp, depvar=depvar, adjR2=TRUE, M=M, H0=FALSE, estSize=estSize, estIMSize=FALSE)
    rownames(rval3$tblSize) <- gsub("Size", "Power", rownames(rval3$tblSize))
    tblCoef[10, , i] <- rval3$tblSize[2, 1:2]
}
dimnames(tblCoef) <- list(c(rownames(rval2$tblCoef), rownames(rval2$tblSize)[1:2], rownames(rval3$tblSize)[2]), regvars, depvars)
dimnames(tblR2) <- list(rownames(rval2$tblR2), c("R^2_1", "", "R^2_2", "", "R^2_2 - R^2_1", ""), depvars)

## tables
## R^2
print(round(tblR2, digi=2))
for (i in seq_along(depvars)) {
    depvar <- depvars[i]
    ## cat("\\multicolumn{4}{l}{\\emph{", switch(depvar, xr2="Two", xr3="Three", xr4="Four", xr5="Five"), "-year bond}} \\\\ \n", sep="")
    cat("\\emph{", switch(depvar, xr2="Two", xr3="Three", xr4="Four", xr5="Five"), "-year bond} \\\\ \n", sep="")
    cat("                          ", sprintf("& %3.2f ", tblR2[1, c(1,3,5), i]), "\n")
    cat("                          ", sprintf("& %3.2f", tblR2[2, c(1,3,5), i]), "\n")
    cat("                          ", sprintf("& (%3.2f, %3.2f)", tblR2[3, c(1,3,5), i], tblR2[3, c(2,4,6), i]), "\n")
}

## Coefficients
## two columns per dependent variable
tbl <- cbind(tblCoef[,,1], tblCoef[,,2], tblCoef[,,3], tblCoef[,,4])
print(round(tbl, 3))
for (depvar in depvars)
    cat("& \\multicolumn{2}{l}{\\emph{", switch(depvar, xr2="Two", xr3="Three", xr4="Four", xr5="Five"), "-year bond}}", sep="")
cat("\\\\ \n")
xtbl <- xtable(tbl, digi=3)
print(xtbl, include.rownames=TRUE, include.colnames=TRUE, only.contents=TRUE, sanitize.text.function=function(x){x})

