#				
## 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)
source("gestCov_txRand_VESS.R")
source("gestCov_txRand_ZIPI.R")
source("gestCov_txRand_PUR.R")
source("../gestCov_txRand.R")
#source("../gestCov_txRandtemp.R")
#source("temp.R")
  
psi10 <- 1.6	
psi11 <- 0	
psi20 <- -3	
psi21 <- 2
psi22 <- 1	
psi23 <- 0	

bias.mat <- NULL

psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)
get.vess <- 0
get.zipi <- 1
get.pur <- 0

#for(psi20 in seq(-4,-2,.1)) {
#bias.vec <- NULL
for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
#psi<-c(psi10,psi11,psi20,psi21)
results <- psi
   psi.est90 <- var.est90 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.est80 <- var.est80 <- NULL
   psi.estPUR <- var.estPUR <- NULL
   psi.estJR <- var.estJR <- NULL
   wrong.rule <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(200,n.sim,200)) print(j)
      L1 <- rnorm(n,0,3)
      A1 <- rbinom(n,1,alpha[1])  
#      L2 <- rnorm(n,.75,.5)
      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      

      if(get.vess) {
       temp <- VESS_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.95)) # 90% CI exclusion
       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 <- VESS_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.975)) # 95% CI exclusion
       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 <- VESS_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.90)) # 80% CI exclusion
       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,])) 
      }
      if(get.zipi) {
       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.90)) # 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,])) 
      }
      if(get.pur) {
       temp <- PUR_rand.fcn(Y,L,A,return.U=F) 
       psi.estPUR <- cbind(psi.estPUR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
       var.estPUR <- cbind(var.estPUR,c(temp$psi.var[2,1:2],temp$psi.var[1,])) 
      }
      temp <- g_estCov.fcn(Y,L,A)
      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,])) 
   }

   print(paste("n = ",n,sep=""))
   if(get.vess) {
    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
    t190.bias <- abs(apply(psi.est90-psi,1,mean))
    t190.t <- abs(apply((psi.est90-psi)/sqrt(var.est90),1,mean))
    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
    t195.bias <- abs(apply(psi.est95-psi,1,mean))
    t195.t <- abs(apply((psi.est95-psi)/sqrt(var.est95),1,mean))
    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
    t180.bias <- abs(apply(psi.est80-psi,1,mean))
    t180.t <- abs(apply((psi.est80-psi)/sqrt(var.est80),1,mean))
    res.vess <- rbind(cbind(t1.90,t2.90,t190.t,t3.90,t4.90),cbind(t1.95,t2.95,t195.t,t3.95,t4.95))
   }
   if(get.zipi) {
    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
    t190.bias <- abs(apply(psi.est90-psi,1,mean))
    t190.t <- abs(apply((psi.est90-psi)/sqrt(var.est90),1,mean))
    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
    t195.bias <- abs(apply(psi.est95-psi,1,mean))
    t195.t <- abs(apply((psi.est95-psi)/sqrt(var.est95),1,mean))
    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
    t180.bias <- abs(apply(psi.est80-psi,1,mean))
    t180.t <- abs(apply((psi.est80-psi)/sqrt(var.est80),1,mean))
    res.zipi <- rbind(cbind(t1.90,t2.90,t190.t,t3.90,t4.90),cbind(t1.95,t2.95,t195.t,t3.95,t4.95))
   }
   if(get.pur) {
    t1p <- apply(psi.estPUR,1,mean)
    t0p <- apply(psi.estPUR,1,sd)
    t2p <- apply(sqrt(var.estPUR),1,mean) # calculated SE
    t3p <- apply(sqrt(var.estPUR + (psi.estPUR - psi)^2),1,mean)  # rMSE
    t4p <- apply((psi.estPUR - 1.96*sqrt(var.estPUR) < psi) & (psi.estPUR + 1.96*sqrt(var.estPUR) > psi),1,sum)/n.sim*100 # Wald CI coverage
    t1P.bias <- abs(apply(psi.estPUR-psi,1,mean))
    t1P.t <- abs(apply((psi.estPUR-psi)/sqrt(var.estPUR),1,mean))
    res.pur <- cbind(t1p,t2p,t1P.t,t3p,t4p)
   }

   t1g <- apply(psi.estJR,1,mean)
   t1g.bias <- abs(apply(psi.estJR-psi,1,mean))
   t1g.t <- abs(apply((psi.estJR-psi)/sqrt(var.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 - qt(.975,100)*sqrt(var.estJR) < psi) & (psi.estJR + qt(.975,100)*sqrt(var.estJR) > psi),1,sum)/n.sim*100 # Wald CI coverage

   if(get.vess) {
    results <- cbind(t1.80,t2.80,t180.t,t3.80,t4.80,t1g,t2g,t1g.t,t3g,t4g)
    print(round(results,3))
    print(round(res.vess,3))
   }
   if(get.zipi) {
    results <- cbind(t1.80,t2.80,t180.t,t3.80,t4.80,t1g,t2g,t1g.t,t3g,t4g)
    print(round(results,3))
    print(round(res.zipi,3))
   }

   if(get.pur) {
    results <- cbind(t1p,t2p,t1P.t,t3p,t4p,t1g,t2g,t1g.t,t3g,t4g)
    print(round(results,3))
   }

print("Abs. Bias")
   if(get.vess) {
    print("VESS")
    print(round(t180.bias,3))
    print(round(t190.bias,3))
    print(round(t195.bias,3))
   }
   if(get.zipi) {
    print("ZIPI")
    print(round(t180.bias,3))
    print(round(t190.bias,3))
    print(round(t195.bias,3))
   }
   if(get.pur) {
    print("URW")
    print(round(t1P.bias,3))
   }
    print("g-est")
    print(round(t1g.bias,3))


  if(!get.vess & !get.zipi & !get.pur) {
#    results <- cbind(results,cbind(t1g,t2g,t3g,t4g))
    results <- cbind(t1g,t2g,t1g.bias,t1g.t,t3g,t4g)
    print(round(results,3))
#    bias.vec <- c(bias.vec,t1g[1] - psi10)
  }
}
#bias.mat <- cbind(bias.mat,bias.vec)
#}
#bias.mat <- abs(bias.mat)

#postscript("bias_by_psi20a.eps")
#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),zlim=c(-.07,.07))
#dev.off()
#postscript("bias_by_psi21VESSb.eps")
#filled.contour(seq(250,1000,50),seq(-4,-2,.1),bias.mat,xlab="n",ylab=expression(psi[20]),col=rainbow(35),levels=seq(0,.07,.0025))
#dev.off()
#filled.contour(seq(250,1000,50),seq(1,3,.1),bias.mat,xlab="n",ylab=expression(psi[21]),col=rainbow(35),levels=seq(0,.07,.0025))
#round(cbind(t1p,t2p,t3p,t4p,t1g,t2g,t3g,t4g),3)