## test alternative spanning hypothesis

rm(list=ls())
graphics.off()
source("R/robust_fns.r")
source("R/var_fns.r")
library(xtable)
set.seed(616)
M <- 5000  # number of bootstrap replications

tbl1 <- data.frame(matrix("", 6, 7), stringsAsFactors=FALSE) # R^2 - each study 2 rows
colnames(tbl1) <- c("Study", "a) N=3", "N=4", "N=5", "b) N=3", "N=4", "N=5")
tbl2 <- data.frame(matrix("", 9, 8), stringsAsFactors=FALSE) # p-values/size - each study 3 rows
colnames(tbl2) <- c("Study", "Statistic", "a) N=3", "N=4", "N=5", "b) N=3", "N=4", "N=5")
tbl2[[2]] <- as.character(tbl2[[2]])

for (latesample in c(FALSE, TRUE)) {
    coladd <- ifelse(latesample, 3, 0)
    ## JPS
    tbl1[1, 1] <- "JPS"
    tbl2[1, 1] <- "JPS"
    tbl2[1:3, 2] <- c("HAC $p$-value", "HAC size", "BC bootstrap $p$-value")
    jpsdata <- getJPSdata(end.year=ifelse(latesample, 2016, 2007))
    macro.vars <- c("GRO", "INF")
    for (N in 3:5) {
        PCs <- paste0("PC", 1:N)
        K <- 2+N
        dgp <- getBootDGP(PCs, macro.vars, jpsdata, control=list(BC=TRUE))
        fmla1 <- formula(paste("xr.avg ~ ", paste(PCs, collapse="+")))
        fmla2 <- formula(paste("xr.avg ~ ", paste(c(PCs, macro.vars), collapse="+")))
        rval <- bootstrapTest(fmla1, fmla2, jpsdata, dgp, depvar="xr.avg", adjR2=TRUE, M=M, estSize=FALSE)
        tbl1[1, N-1+coladd] <- sprintf("%4.2f", rval$tblR2[1, 5])
        tbl1[2, N-1+coladd] <- sprintf("(%4.2f, %4.2f)", rval$tblR2[3, 5], rval$tblR2[3, 6])
        tbl2[1:3, N+coladd] <- sprintf("%5.3f", c(rval$tblCoef[3, K+1], rval$tblSize[1, K+1], rval$tblCoef[5, K+1]))
        print(tbl1)
        print(tbl2)
    }
    ## LN
    tbl1[3, 1] <- "LN"
    tbl2[4, 1] <- "LN"
    tbl2[4:6, 2] <- c("HAC $p$-value", "HAC size", "Bootstrap $p$-value")
    if (latesample) {
        lndata <- getLNdata(beg=198501, end=201612)  ## later sample period, our data
    } else {
        lndata <- getLNdata()
    }
    macro.vars <-  paste0("F", 1:8)
    for (N in 3:5) {
        PCs <- paste0("PC", 1:N)
        K <- 8+N
        dgp <- getBootDGP(PCs, macro.vars, lndata)
        fmla1 <- formula(paste("xr.avg ~", paste(PCs, collapse="+")))
        fmla2 <- formula(paste("xr.avg ~", paste(c(PCs, macro.vars), collapse="+")))
        rval <- bootstrapTest(fmla1, fmla2, lndata, dgp, depvar="xr.avg", adjR2=TRUE, M=M, estSize=FALSE)
        tbl1[3, N-1+coladd] <- sprintf("%4.2f", rval$tblR2[1, 5])
        tbl1[4, N-1+coladd] <- sprintf("(%4.2f, %4.2f)", rval$tblR2[3, 5], rval$tblR2[3, 6])
        tbl2[4:6, N+coladd] <- sprintf("%5.3f", c(rval$tblCoef[3, K+1], rval$tblSize[1, K+1], rval$tblCoef[5, K+1]))
        print(tbl1)
        print(tbl2)
    }
    ## CPO
    tbl1[5, 1] <- "CPO"
    tbl2[7, 1] <- "CPO"
    tbl2[7:9, 2] <- c("RR $p$-value", "RR size", "BC bootstrap $p$-value")
    if (latesample) {
        df <- getOwnCPoData(198501, 201612)
    } else {
        df <- getOwnCPoData()
    }
    for (N in 3:5) {
        PCs <- paste0("PC", 1:N)
        K <- 1+N
        dgp <- getBootDGP(PCs, "trend", df, control=list(BC=TRUE, weightedReturns=TRUE, monthlyReturns=TRUE))
        fmla1 <- formula(paste("xr.avg ~", paste(PCs, collapse="+")))
        fmla2 <- formula(paste("xr.avg ~", paste(c(PCs, "trend"), collapse="+")))
        rval <- bootstrapTestRR(fmla1, fmla2, df, dgp, M=M, estBootstrapSize=FALSE)
        tbl1[5, N-1+coladd] <- sprintf("%4.2f", rval$tblR2[1, 5])
        tbl1[6, N-1+coladd] <- sprintf("(%4.2f, %4.2f)", rval$tblR2[3, 5], rval$tblR2[3, 6])
        tbl2[7:9, N+coladd] <- sprintf("%5.3f", c(rval$tblCoef[5, K], rval$tblSize[2, K], rval$tblCoef[9, K]))
        print(tbl1)
        print(tbl2)
    }
}

print(tbl1)

print(tbl2)

