## Appendix -- Cooper-Priestley
## Table B.5

rm(list=ls())
library(sandwich)
source("R/robust_fns.r")
source("R/var_fns.r")
library(xtable)
library(vars)

cprdata <- getCPRdata()

M <- 5000

## CP factor
cprdata$cp <- fitted(lm(xr.avg~f1+f2+f3+f4+f5, data=cprdata, na.action="na.exclude"))

mod <- lm(gap ~ cp, data=cprdata)
print(cor(cprdata$gap, cprdata$cp, use="complete.obs"))
print(summary(mod)$adj.r.squared)

##generating orthogonalized and lagged cp factor
cprdata$cpresid <- residuals(lm(cp ~ gap, data=cprdata, na.action="na.exclude"))

cprdata$gap <- c(NA, cprdata$gap[1:(nrow(cprdata)-1)])
cprdata <- cprdata[!is.na(cprdata$gap),]

## Replication - their Table 9
depvars <- c("xr2", "xr3", "xr4", "xr5")
results <- matrix(NA, 5, 8)
rownames(results) <- c("gap","gap.t","cp","cp.t","adj.r.squared")
colnames(results) <- paste(rep(depvars, each=2), rep(c("gap", "cp"), 4), sep="-")
col <- 1
vcovfn <- function(lm.) NeweyWest(lm., lag=22, prewhite=FALSE)
##vcovfn <- function(lm.) vcov(lm.) ## OLS standard errors - this closely matches their t-stats
for (depvar in depvars) {
    lm.gap <- lm(get(depvar) ~ gap, data=cprdata)
    b.gap <- lm.gap$coef
    se.gap <- sqrt(diag(vcovfn(lm.gap)))
    lm.cp <- lm(get(depvar) ~ gap + cpresid, data=cprdata)
    b.cp <- lm.cp$coef
    se.cp <- sqrt(diag(vcovfn(lm.cp)))
    results[,col] <- c(b.gap[2], abs(b.gap[2]/se.gap[2]), NA, NA, summary(lm.gap)$adj.r.squared)
    results[,col+1] <- c(b.cp[2], abs(b.cp[2]/se.cp[2]), b.cp[3],
                         abs(b.cp[3]/se.cp[3]), summary(lm.cp)$adj.r.squared)
    col <- col+2
}
print(round(results, digi=3))

## our results
dgp <- getBootDGP(c("PC1", "PC2", "PC3"), "gap", cprdata, control=list(lag1=1, lag2=1))

fmla1 <- get(depvar) ~ PC1 + PC2 + PC3
fmla2 <- get(depvar) ~ PC1 + PC2 + PC3 + gap

depvar <- "xr5"
fmlas <- list(get(depvar) ~ gap,
              get(depvar) ~ gap + cpresid,
              get(depvar) ~ gap + cp,
              get(depvar) ~ gap + PC1 + PC2 + PC3)

vcovfn <- function(lm.) NeweyWest(lm., lag=22, prewhite=FALSE)
for (i in seq_along(fmlas)) {
    fmla <- fmlas[[i]]
    tbl <- matrix(NA, 6, 6)
    colnames(tbl) <- c("gap", "cpresid", "cp", "PC1", "PC2", "PC3")
    tblMueller <- getMuellerTable(fmla, cprdata)
    colind <- match(colnames(tblMueller), colnames(tbl))
    mod <- lm(fmla, data=cprdata)
    b <- mod$coef
    SEs <- sqrt(diag(vcov(mod)))
    tstats <- abs(b/SEs)
    tbl[1, colind] <- b[-1]
    tbl[2, colind] <- tstats[-1]
    SEs <- sqrt(diag(vcovfn(mod)))
    tstats <- abs(b/SEs)
    tbl[3, colind] <- tstats[-1]
    tbl[4, colind] <- pt(tstats[-1], mod$df, lower.tail=FALSE)*2
    tbl[5:6, colind] <- tblMueller[2:3,]
    rownames(tbl) <- c("Coefficients", "OLS $t$-statistic", "HAC $t$-statistic", "HAC $p$-value", "IM $q=8$", "IM $q=16$")
    if (i<4) {
        ## conventional tests for their specifications
        tbl <- tbl[-(5:6),]
    } else {
        ## bootstrap for specification with PCs
        rval <- bootstrapTest(fmla1, fmla2, cprdata, dgp, depvar=depvar, M=M)
        tblBoot <- rval$tblCoef[4:5,1:4]
        tbl <- rbind(tbl[1:4,], NA, NA, tbl[5:6,])
        colind <- match(colnames(tblBoot), colnames(tbl))
        tbl[5:6, colind] <- tblBoot
        rownames(tbl)[5:6] <- rownames(tblBoot)
    }
    print(round(tbl, 3))
    ## printLatexTable(tbl)
    ## xtbl <- xtable(tbl, digi=3)
    ## print(xtbl, include.rownames=TRUE, include.colnames=FALSE, only.contents=TRUE, hline.after=nrow(tbl),
    ##       sanitize.text.function=function(x){x})
}

