#				
## 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 <- seq(250,600,50)
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
bias.mat80 <- NULL
bias.mat95 <- NULL

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

for(psi20 in seq(-4,-2,.1)) {
bias.vec <- NULL
bias.vec80 <- NULL
bias.vec95 <- NULL
for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- var.est80 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.estPUR <- var.estPUR <- NULL
   psi.estJR <- var.estJR <- NULL
   wrong.rule <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(500,n.sim,500)) 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 <- 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,]))   
       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,]))   

      temp <- g_estCov.fcn(Y,L,A)
      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
   }

   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi20,sep=""))
   t1.80 <- apply(psi.est80,1,mean)
   bias.vec80 <- c(bias.vec80,t1.80[1] - psi10)
   t1.95 <- apply(psi.est95,1,mean)
   bias.vec95 <- c(bias.vec95,t1.95[1] - psi10)

   t1g <- apply(psi.estJR,1,mean)
   bias.vec <- c(bias.vec,t1g[1] - psi10)
}
bias.mat <- cbind(bias.mat,bias.vec)
bias.mat80 <- cbind(bias.mat80,bias.vec80)
bias.mat95 <- cbind(bias.mat95,bias.vec95)
}
bias.mat <- abs(bias.mat)
bias.mat80 <- abs(bias.mat80)
bias.mat95 <- abs(bias.mat95)

print("Max bias")
print(paste("g-est ",max(as.vector(bias.mat)),sep=""))
print(paste("ZIPI 80 ",max(as.vector(bias.mat80)),sep=""))
print(paste("ZIPI 95 ",max(as.vector(bias.mat95)),sep=""))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for(psi21 in seq(1,3,.1)) {
bias.vec <- NULL
bias.vec80 <- NULL
bias.vec95 <- NULL
for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- var.est80 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.estPUR <- var.estPUR <- NULL
   psi.estJR <- var.estJR <- NULL
   wrong.rule <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(500,n.sim,500)) 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 <- 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,]))   
       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,]))   

      temp <- g_estCov.fcn(Y,L,A)
      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
   }

   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi21,sep=""))
   t1.80 <- apply(psi.est80,1,mean)
   bias.vec80 <- c(bias.vec80,t1.80[1] - psi10)
   t1.95 <- apply(psi.est95,1,mean)
   bias.vec95 <- c(bias.vec95,t1.95[1] - psi10)

   t1g <- apply(psi.estJR,1,mean)
   bias.vec <- c(bias.vec,t1g[1] - psi10)
}
bias.mat <- cbind(bias.mat,bias.vec)
bias.mat80 <- cbind(bias.mat80,bias.vec80)
bias.mat95 <- cbind(bias.mat95,bias.vec95)
}
bias.mat <- abs(bias.mat)
bias.mat80 <- abs(bias.mat80)
bias.mat95 <- abs(bias.mat95)

print("Max bias")
print(paste("g-est ",max(as.vector(bias.mat)),sep=""))
print(paste("ZIPI 80 ",max(as.vector(bias.mat80)),sep=""))
print(paste("ZIPI 95 ",max(as.vector(bias.mat95)),sep=""))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for(psi22 in seq(0,2,.1)) {
bias.vec <- NULL
bias.vec80 <- NULL
bias.vec95 <- NULL
for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- var.est80 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.estPUR <- var.estPUR <- NULL
   psi.estJR <- var.estJR <- NULL
   wrong.rule <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(500,n.sim,500)) 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 <- 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,]))   
       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,]))   

      temp <- g_estCov.fcn(Y,L,A)
      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
   }

   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi22,sep=""))
   t1.80 <- apply(psi.est80,1,mean)
   bias.vec80 <- c(bias.vec80,t1.80[1] - psi10)
   t1.95 <- apply(psi.est95,1,mean)
   bias.vec95 <- c(bias.vec95,t1.95[1] - psi10)

   t1g <- apply(psi.estJR,1,mean)
   bias.vec <- c(bias.vec,t1g[1] - psi10)
}
bias.mat <- cbind(bias.mat,bias.vec)
bias.mat80 <- cbind(bias.mat80,bias.vec80)
bias.mat95 <- cbind(bias.mat95,bias.vec95)
}
bias.mat <- abs(bias.mat)
bias.mat80 <- abs(bias.mat80)
bias.mat95 <- abs(bias.mat95)

print("Max bias")
print(paste("g-est ",max(as.vector(bias.mat)),sep=""))
print(paste("ZIPI 80 ",max(as.vector(bias.mat80)),sep=""))
print(paste("ZIPI 95 ",max(as.vector(bias.mat95)),sep=""))



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for(psi23 in seq(-1,1,.1)) {
bias.vec <- NULL
bias.vec80 <- NULL
bias.vec95 <- NULL
for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- var.est80 <- NULL
   psi.est95 <- var.est95 <- NULL
   psi.estPUR <- var.estPUR <- NULL
   psi.estJR <- var.estJR <- NULL
   wrong.rule <- NULL
   for(j in 1:n.sim) {
if(j %in% seq(500,n.sim,500)) 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 <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.90)) # 80% CI zeroing
       psi.est90 <- cbind(psi.est90,c(temp$psi.hat[2,1:2],temp$psi.hat[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,]))   

      temp <- g_estCov.fcn(Y,L,A)
      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
   }

   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi23,sep=""))
   t1.80 <- apply(psi.est80,1,mean)
   bias.vec80 <- c(bias.vec80,t1.80[1] - psi10)
   t1.95 <- apply(psi.est95,1,mean)
   bias.vec95 <- c(bias.vec95,t1.95[1] - psi10)

   t1g <- apply(psi.estJR,1,mean)
   bias.vec <- c(bias.vec,t1g[1] - psi10)
}
bias.mat <- cbind(bias.mat,bias.vec)
bias.mat80 <- cbind(bias.mat80,bias.vec80)
bias.mat95 <- cbind(bias.mat95,bias.vec95)
}
bias.mat <- abs(bias.mat)
bias.mat80 <- abs(bias.mat80)
bias.mat95 <- abs(bias.mat95)

print("Max bias")
print(paste("g-est ",max(as.vector(bias.mat)),sep=""))
print(paste("ZIPI 80 ",max(as.vector(bias.mat80)),sep=""))
print(paste("ZIPI 95 ",max(as.vector(bias.mat95)),sep=""))