## CPo case study
## How is inference on coefficient on PC1 affected by including a trend?

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

df <- getOwnCPoData()   # original sample

X1.names <- c("PC1", "PC2", "PC3")
X2.names <- "trend"

dgp <- getBootDGP(X1.names, X2.names, df, control=list(BC=TRUE, weightedReturns=TRUE, monthlyReturns=TRUE))

### bootstrap simulations and results
rval <- bootstrapTestRRPC1(df, dgp, M=M)
attach(rval)

## tables with results
nms <- c("True beta1", "Mean(b1)", "Bias(b1)", "Size beta1=beta1_true", "Size beta1=0", "Mean SE(b1)", "RMSE(b1)", "SD(b1)")
tbl <- matrix(NA, length(nms), 2)  ## 3 rows for data, 2 rows for bootstrap
colnames(tbl) <- c("w/o x2", "w/ x2")
rownames(tbl) <- nms
tbl[1, ] <- beta1
tbl[2, ] <- colMeans(bootCoefs)
tbl[3, ] <- colMeans(bootCoefs - beta1)
## tests
tbl[4, ] <- colMeans(abs(bootRRCoefs - beta1)/bootSEs > qnorm(0.975)) # Size of test for b1=beta1
tbl[5, ] <- colMeans(abs(bootRRCoefs)/bootSEs > qnorm(0.975)) # Size of test for b1=beta1
## standard errors
tbl[6, ] <- colMeans(bootSEs)
tbl[7, ] <- sqrt(colMeans((bootCoefs - beta1)^2))
tbl[8, ] <- apply(bootCoefs, 2, sd)
print(round(tbl,4))

## tbl[2, ] <- colMeans(bootRRCoefs)
## tbl[3, ] <- colMeans(bootRRCoefs - beta1)
## tbl[7, ] <- sqrt(colMeans((bootRRCoefs - beta1)^2))
## tbl[8, ] <- apply(bootRRCoefs, 2, sd)
## print(round(tbl,4))

tbl2 <- matrix(NA, 1, 2)
colnames(tbl2) <- c("PCs and tau", "tau only")
rownames(tbl2) <- c("rejections")
tbl2[1,] <- colMeans(bootTauStats > qnorm(0.975))
print(round(tbl2, 4))
