#				
## 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,1000,50)
source("gestCov_txRand_ZIPI.R")
 
psi20 <- -3	
psi21 <-  2
psi22 <-  1	
psi23 <-  0	
psi10 <- psi20	
psi11 <- psi21
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)
bias0.mat80 <- NULL

#for(psi20 in seq(-4,-2,.1)) {
#psi10 <- psi20	
#psi11 <- psi21
#bias0.vec80 <- NULL
#for(n in n.sample) {
#psi<-c(psi20,psi21,psi22,psi23)
#results <- psi
#   psi.est80 <- 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 <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80%CI
#      icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
#      psi.est80 <- cbind(psi.est80,c(icw.est,temp$psi.hat[1,3:4]))
#   }

#   t1.80 <- apply(psi.est80,1,mean)
#   bias0.vec80 <- c(bias0.vec80,t1.80[1] - psi20)
#   print(paste("n = ",n,sep=""))
#   print(paste("psi = ",psi20,sep=""))
#}
#bias0.mat80 <- cbind(bias0.mat80,abs(bias0.vec80))
#}

#postscript("bias0_psi20ZIPI80.eps")
#filled.contour(seq(250,1000,50),seq(-4,-2,.1),bias0.mat80,xlab="n",
#    ylab=expression(psi[20]),col=rainbow(50),levels=seq(0,.18,.005))
#dev.off()


psi20 <- -3	
psi21 <-  2
psi22 <-  1	
psi23 <-  0	
psi10 <- psi20	
psi11 <- psi21
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)
bias0.mat80 <- NULL

for(psi21 in seq(1,3,.1)) {
psi10 <- psi20	
psi11 <- psi21
bias0.vec80 <- NULL
for(n in n.sample) {
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- 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 <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80%CI
      icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
      psi.est80 <- cbind(psi.est80,c(icw.est,temp$psi.hat[1,3:4]))

   }

   t1.80 <- apply(psi.est80,1,mean)

   bias0.vec80 <- c(bias0.vec80,t1.80[1] - psi20)
   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi21,sep=""))

}
bias0.mat80 <- cbind(bias0.mat80,abs(bias0.vec80))
}

postscript("bias0_psi21ZIPI80.eps")
filled.contour(seq(250,1000,50),seq(1,3,.1),bias0.mat80,xlab="n",
    ylab=expression(psi[21]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()

psi20 <- -3	
psi21 <-  2
psi22 <-  1	
psi23 <-  0	
psi10 <- psi20	
psi11 <- psi21
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)
bias0.mat80 <- NULL

for(psi22 in seq(0,2,.1)) {
psi10 <- psi20	
psi11 <- psi21
bias0.vec80 <- NULL
for(n in n.sample) {
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- 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 <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80%CI
      icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
      psi.est80 <- cbind(psi.est80,c(icw.est,temp$psi.hat[1,3:4]))

   }

   t1.80 <- apply(psi.est80,1,mean)

   bias0.vec80 <- c(bias0.vec80,t1.80[1] - psi20)
   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi22,sep=""))

}
bias0.mat80 <- cbind(bias0.mat80,abs(bias0.vec80))
}

postscript("bias0_psi22ZIPI80.eps")
filled.contour(seq(250,1000,50),seq(0,2,.1),bias0.mat80,xlab="n",
    ylab=expression(psi[22]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()


psi20 <- -3	
psi21 <-  2
psi22 <-  1	
psi23 <-  0	
psi10 <- psi20	
psi11 <- psi21
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
alpha<-c(.5,.5)
bias0.mat80 <- NULL

for(psi23 in seq(-1,1,.1)) {
psi10 <- psi20	
psi11 <- psi21
bias0.vec80 <- NULL
for(n in n.sample) {
psi<-c(psi20,psi21,psi22,psi23)
results <- psi
   psi.est80 <- 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 <- ZIPI_rand.fcn(Y,L,A,return.U=F,z.val=qnorm(0.9)) # 80%CI
      icw.est <- (temp$psi.hat[2,1:2]/temp$psi.var[2,1:2] + temp$psi.hat[1,1:2]/temp$psi.var[1,1:2])/(1/temp$psi.var[2,1:2] + 1/temp$psi.var[1,1:2])
      psi.est80 <- cbind(psi.est80,c(icw.est,temp$psi.hat[1,3:4]))

   }

   t1.80 <- apply(psi.est80,1,mean)

   bias0.vec80 <- c(bias0.vec80,t1.80[1] - psi20)
   print(paste("n = ",n,sep=""))
   print(paste("psi = ",psi23,sep=""))

}
bias0.mat80 <- cbind(bias0.mat80,abs(bias0.vec80))
}

postscript("bias0_psi23ZIPI80.eps")
filled.contour(seq(250,1000,50),seq(-1,1,.1),bias0.mat80,xlab="n",
    ylab=expression(psi[23]),col=rainbow(50),levels=seq(0,.18,.005))
dev.off()