## simulation results for paper
## Table 1

rm(list=ls())
set.seed(616)
library(xtable)
source("R/robust_fns.r")  ## for getMueller and asymTstat

doWork <- function(n=100, delta=0, c1=0, c2=0, rho=0.99, rho1=rho, rho2=rho, sig1=1, sig2=1, sigu=1, trace=FALSE) {
    M <- 1000000
    b <- matrix(NA, M, 2)
    se <- matrix(NA, M, 2)
    tstatboot <- numeric(M)
    pvalsIM <- numeric(M)
    for (i in 1:M) {
        v <- matrix(rnorm((n+2)*3), n+2, 3)
        eps <- cbind(sig1*(delta*v[,3]+sqrt(1-delta^2)*v[,1]), # eps_1t
                     sig2*v[,2], # eps_2t
                     sigu*v[,3])   # u_t
        x1 <- numeric(n+2)
        x2 <- numeric(n+2)
        ## instead draw from unconditional distribution
        ## x1[1] <- c1/(1-rho1) + rnorm(1, sd=sig1/sqrt(1-rho1^2))
        ## x2[1] <- c2/(1-rho2) + rnorm(1, sd=sig2/sqrt(1-rho2^2))
        for (t in 2:(n+2)) {
            x1[t] <- c1 + rho1*x1[t-1] + eps[t, 1]
            x2[t] <- c2 + rho2*x2[t-1] + eps[t, 2]
        }
        x1 <- x1[-1]; x2 <- x2[-1]; eps <- eps[-1,]  # get rid of first obs, where x's are zero
        y <- eps[2:(n+1), 3]  # beta0=beta1=beta2=0
        df <- data.frame(y = y, x1 = x1[1:n], x2 = x2[1:n], trend = 1:n)
        mod12 <- lm(y ~ x1 + x2, data=df)
        b[i, ] <- mod12$coef[-1]
        se[i, ] <- sqrt(diag(vcov(mod12)))[-1]
        ## bootstrap
        ## 1. estimate AR
        ar.x1 <- ar.ols(df$x1, aic=FALSE, order.max=1, demean=FALSE, intercept=TRUE)
        ar.x2 <- ar.ols(df$x2, aic=FALSE, order.max=1, demean=FALSE, intercept=TRUE)
        ## 2. generate single bootstrap sample
        ind <- 1 + sample(n-1, n+2, TRUE)
        e1boot <- ar.x1$resid[ind]
        e2boot <- ar.x2$resid[ind]
        mod.r <- lm(y ~ x1, data=df)
        uboot <- mod.r$resid[ind-1]
        x1boot <- numeric(n+2)
        x2boot <- numeric(n+2)
        for (t in 2:(n+2)) {
            x1boot[t] <- ar.x1$x.intercept + ar.x1$ar[,,]*x1boot[t-1] + e1boot[t]
            x2boot[t] <- ar.x2$x.intercept + ar.x2$ar[,,]*x2boot[t-1] + e2boot[t]
        }
        x1boot <- x1boot[-1]; x2boot <- x2boot[-1]; e1boot <- e1boot[-1]; e2boot <- e2boot[-1]; uboot <- uboot[-1]
        yboot <- c(NA, ar.x1$ar[,,]*x1boot[1:n]) + uboot
        bootdf <- data.frame(y = yboot[2:(n+1)], x1 = x1boot[1:n], x2 = x2boot[1:n])
        ## 3. estimate predictive regression on bootstrap sample
        modboot <- lm(y ~ x1 + x2, data=bootdf)
        tstatboot[i] <- (modboot$coef/sqrt(diag(vcov(modboot))))[3]
    }
    t <- b/se
    cv <- qt(0.975, n-3)
    rval <- c(colMeans(b), 100*(mean(se[,2])/sd(b[,2])-1), mean(abs(t[,2])>cv))
    ## asymptotic approximation using local-to-unity
    c <- n*(rho-1)
    t.mc <- replicate(M, asymTstat(delta, c))
    rval <- c(rval, mean(abs(t.mc)>cv))
    ## bootstrap
    cvboot <- quantile(abs(tstatboot), 0.95)
    rval <- c(rval, mean(abs(t[,2])>cvboot))
}

cat("# Simulating...\n")
## first panel: no trends
tbl <- matrix(NA, 5, 8)
colnames(tbl) <- c("rho", "delta", "b(b1)", "b(b2)", "SE bias (%)", "size sim", "size asym", "bootstrap")
rval <- doWork(rho=0.99, delta=0)
tbl[1, ] <- c(0.99, 0, rval)
tbl[2, ] <- c(0, 1, doWork(rho=0, delta=1))
tbl[3, ] <- c(0.9, 1, doWork(rho=0.9, delta=1))
tbl[4, ] <- c(0.99, 0.8, doWork(rho=0.99, delta=0.8))
tbl[5, ] <- c(0.99, 1, doWork(rho=0.99, delta=1))
print(round(tbl, digi=3))

## filename <- "tables/sim_1.tex"
## cat("*** writing table ", filename, "\n")
## sink(filename)
## latex.tbl <- xtable(tbl, digits=c(NA, c(2,1,3,3,1,3,3,3,3)))
## print(latex.tbl, include.rownames=FALSE, include.colnames=FALSE, only.contents=TRUE, hline.after=NULL)
## sink()

## second panel: c_2 = 1
tbl[] <- NA
tbl[1, ] <- c(0.99, 0, doWork(rho=0.99, delta=0, c2=1))
tbl[2, ] <- c(0, 1, doWork(rho=0, delta=1, c2=1))
tbl[3, ] <- c(0.9, 1, doWork(rho=0.9, delta=1, c2=1))
tbl[4, ] <- c(0.99, 0.8, doWork(rho=0.99, delta=0.8, c2=1))
tbl[5, ] <- c(0.99, 1, doWork(rho=0.99, delta=1, c2=1))
print(round(tbl, digi=3))
tbl[,7] <- NA

## filename <- "tables/sim_2.tex"
## cat("*** writing table ", filename, "\n")
## sink(filename)
## latex.tbl <- xtable(tbl, digits=c(NA, c(2,1,3,3,1,3,3,3,3)))
## print(latex.tbl, include.rownames=FALSE, include.colnames=FALSE, only.contents=TRUE, hline.after=NULL)
## sink()

## third panel: c_1 = c_2 = 1
tbl[] <- NA
tbl[1, ] <- c(0.99, 0, doWork(rho=0.99, delta=0, c1=1, c2=1))
tbl[2, ] <- c(0, 1, doWork(rho=0, delta=1, c1=1, c2=1))
tbl[3, ] <- c(0.9, 1, doWork(rho=0.9, delta=1, c1=1, c2=1))
tbl[4, ] <- c(0.99, 0.8, doWork(rho=0.99, delta=0.8, c1=1, c2=1))
tbl[5, ] <- c(0.99, 1, doWork(rho=0.99, delta=1, c1=1, c2=1))
print(round(tbl, digi=3))
tbl[,7] <- NA

## filename <- "tables/sim_3.tex"
## cat("*** writing table ", filename, "\n")
## sink(filename)
## latex.tbl <- xtable(tbl, digits=c(NA, c(2,1,3,3,1,3,3,3,3)))
## print(latex.tbl, include.rownames=FALSE, include.colnames=FALSE, only.contents=TRUE, hline.after=NULL)
## sink()


