#				
## Simulated data to test variance  - random treatment, 2 intervals (comparing at my covariance with (3.14))
#				

rm(list=ls())
library(MASS)

n.sim <- 1000
n.sample <- c(250,500,1000)
n.boot <- 500
source("../gestCov_txRand.R")
source("gestCov_txRand_ZIPI.R")
 
psi10 <- 1.6	
psi11 <- 0	
psi20 <- -3	
psi21 <- 2
psi22 <- 1	
psi23 <- 0	

psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)

for(n in n.sample) {
results <- psi
   psi.est80 <- var.est80 <- NULL
   psi.est90 <- var.est90 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.est99 <- var.est99 <- NULL
   psi.estJR <- var.estJR <- NULL
   psi.estAAO <- var.estAAO <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(100,n.sim,100)) print(j)
      L1 <- rnorm(n,0,3)
      A1 <- rbinom(n,1,alpha[1])  
      L2 <- round(rnorm(n,.75,.5))
      A2 <- rbinom(n,1,alpha[2])  
      L <- cbind(L1,L2)
      A <- cbind(A1,A2)
      Y.opt <- rnorm(n,1,1) 
      Y <- Y.opt - abs(psi10 + psi11*L1)*(A1 - as.numeric(psi10 + psi11*L1>0))^2 - abs(psi20 + psi21*L2 + psi22*A1 + psi23*L2*A1)*(A2 - as.numeric(psi20 + psi21*L2 + psi22*A1 + psi23*L2*A1>0))^2      

      temp <- g_estCov.fcn(Y,L,A,return.U=F)
      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
#      var.estJR <- cbind(var.estJR,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 

      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80% CI zeroing
      psi.est80 <- cbind(psi.est80,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
#       var.est80 <- cbind(var.est80,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 
      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.95)) # 90% CI zeroing
      psi.est90 <- cbind(psi.est90,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
#       var.est90 <- cbind(var.est90,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 
      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.975)) # 95% CI zeroing
      psi.est95 <- cbind(psi.est95,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
#       var.est95 <- cbind(var.est95,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 
      temp <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.995)) # 99% CI zeroing
      psi.est99 <- cbind(psi.est99,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
#       var.est99 <- cbind(var.est99,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 

## Bootstrap the variance!
      boot.psi <- boot.psi80 <- boot.psi90 <- boot.psi95 <- boot.psi99 <- NULL
      for(i in 1:n.boot) {
       boot.samp <- sample(1:n,replace=T)
       Y.boot <- Y[boot.samp]
       L.boot <- L[boot.samp,]
       A.boot <- A[boot.samp,]

       temp <- g_estCov.fcn(Y.boot,L.boot,A.boot,return.U=F)
       boot.psi <- cbind(boot.psi,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))

       temp <- ZIPI_rand.fcn(Y.boot,L.boot,A.boot,return.U=F,z.val=qnorm(0.9)) # 80% CI zeroing
       boot.psi80 <- cbind(boot.psi80,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   

       temp <- ZIPI_rand.fcn(Y.boot,L.boot,A.boot,return.U=F,z.val=qnorm(0.95)) # 90% CI zeroing
       boot.psi90 <- cbind(boot.psi90,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   

       temp <- ZIPI_rand.fcn(Y.boot,L.boot,A.boot,return.U=F,z.val=qnorm(0.975)) # 95% CI zeroing
       boot.psi95 <- cbind(boot.psi95,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   

       temp <- ZIPI_rand.fcn(Y.boot,L.boot,A.boot,return.U=F,z.val=qnorm(0.995)) # 99% CI zeroing
       boot.psi99 <- cbind(boot.psi99,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   

      }
      var.estJR <- cbind(var.estJR,apply(boot.psi,1,var)) 
      var.est80 <- cbind(var.est80,apply(boot.psi80,1,var)) 
      var.est90 <- cbind(var.est90,apply(boot.psi90,1,var)) 
      var.est95 <- cbind(var.est95,apply(boot.psi95,1,var)) 
      var.est99 <- cbind(var.est99,apply(boot.psi99,1,var)) 
   }

   print(paste("n = ",n,sep=""))
   t1.80 <- apply(psi.est80,1,mean)
   t2.80 <- apply(sqrt(var.est80),1,mean) # calculated SE
   t3.80 <- apply(sqrt(var.est80 + (psi.est80 - psi)^2),1,mean)  # rMSE
   t4.80 <- apply((psi.est80 - 1.96*sqrt(var.est80) < psi) & (psi.est80 + 1.96*sqrt(var.est80) > psi),1,sum)/n.sim*100 # Wald CI coverage

   t1.90 <- apply(psi.est90,1,mean)
   t2.90 <- apply(sqrt(var.est90),1,mean) # calculated SE
   t3.90 <- apply(sqrt(var.est90 + (psi.est90 - psi)^2),1,mean)  # rMSE
   t4.90 <- apply((psi.est90 - 1.96*sqrt(var.est90) < psi) & (psi.est90 + 1.96*sqrt(var.est90) > psi),1,sum)/n.sim*100 # Wald CI coverage

   t1.95 <- apply(psi.est95,1,mean)
   t2.95 <- apply(sqrt(var.est95),1,mean) # calculated SE
   t3.95 <- apply(sqrt(var.est95 + (psi.est95 - psi)^2),1,mean)  # rMSE
   t4.95 <- apply((psi.est95 - 1.96*sqrt(var.est95) < psi) & (psi.est95 + 1.96*sqrt(var.est95) > psi),1,sum)/n.sim*100 # Wald CI coverage

   t1.99 <- apply(psi.est99,1,mean)
   t2.99 <- apply(sqrt(var.est99),1,mean) # calculated SE
   t3.99 <- apply(sqrt(var.est99 + (psi.est99 - psi)^2),1,mean)  # rMSE
   t4.99 <- apply((psi.est99 - 1.96*sqrt(var.est99) < psi) & (psi.est99 + 1.96*sqrt(var.est99) > psi),1,sum)/n.sim*100 # Wald CI coverage

   t1g <- apply(psi.estJR,1,mean)
   t2g <- apply(sqrt(var.estJR),1,mean) # calculated SE
   t3g <- apply(sqrt(var.estJR + (psi.estJR - psi)^2),1,mean)  # rMSE
   t4g <- apply((psi.estJR - 1.96*sqrt(var.estJR) < psi) & (psi.estJR + 1.96*sqrt(var.estJR) > psi),1,sum)/n.sim*100 # Wald CI coverage

   results <- cbind(results,cbind(t1.80,t2.80,t3.80,t4.80,t1g,t2g,t3g,t4g))
   print(round(results,3))
   print(round(cbind(t1.90,t2.90,t3.90,t4.90)[1:2,],3))
   print(round(cbind(t1.95,t2.95,t3.95,t4.95)[1:2,],3))
   print(round(cbind(t1.99,t2.99,t3.99,t4.99)[1:2,],3))

   results <- cbind(results,cbind(t1.90,t2.90,t3.90,t4.90,t1g,t2g,t3g,t4g))
   print(round(results,3))
}