## bootstrap CPo results, but use one-month returns
## - bootstrap both PCs and infl-trend as VARs

rm(list=ls())
source("R/robust_fns.r")
source("R/var_fns.r")
require(sandwich)
set.seed(616)
M <- 5000
df <- getOwnCPoData()

X1.names <- c("PC1", "PC2", "PC3")
X2.names <- "trend"
fmla1 <- get(depvar) ~ PC1 + PC2 + PC3
fmla2 <- get(depvar) ~ PC1 + PC2 + PC3 + trend
depvar <- "xr1m.avg"
vcovfn <- function(lm.) vcovHC(lm., type="HC") # White standard errors

dgp <- getBootDGP(X1.names, X2.names, df, control=list(BC=TRUE, weightedReturns=TRUE, monthlyReturns=TRUE))
rval <- bootstrapTest(fmla1, fmla2, df, dgp, M=M, h=1, vcovfn=vcovfn)
tbl2 <- rbind(rval$tblCoef, rval$tblSize)[,1:4]

## include rows for restricted model
tbl1 <- matrix(NA, 3, 4)
rownames(tbl1) <- c("Coefficient", "t-statistic", "p-value")
colnames(tbl1) <- c("PC1", "PC2", "PC3", "trend")
mod1 <- lm(fmla1, df)
tbl1[1, 1:3] <- mod1$coef[-1]
tstats <- (mod1$coef/sqrt(diag(vcovfn(mod1))))[-1]
tbl1[2, 1:3] <- tstats
tbl1[3, 1:3] <- pt(abs(tstats), mod1$df, lower.tail=FALSE)*2

print(round(rbind(tbl1, tbl2), 3))

tbl <- rval$tblSize[1,,drop=FALSE]
print(tbl)

## no trends
cat("## No drift in X1 - initializing at unconditional mean\n")
dgp$X1.init <- dgp$EX1
rval <- bootstrapTest(fmla1, fmla2, df, dgp, M=M, h=1, vcovfn=vcovfn)
tbl <- rbind(tbl, rval$tblSize[1,])
print(tbl)
cat("## And no drift in X2\n")
dgp$X2.init <- dgp$EX2
rval <- bootstrapTest(fmla1, fmla2, df, dgp, M=M, h=1, vcovfn=vcovfn)
tbl <- rbind(tbl, rval$tblSize[1,])
print(tbl)











