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

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

n.sim <- 500#1000
n.sample <- 250#c(250,500,1000)#seq(250,1000,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 <- 2#1	
psi23 <- 0	

bias.mat <- NULL

alpha<-c(.5,.5)

g.2 <- NULL
psi.estJR <- 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      

      temp <- g_estCov.fcn(Y,L,A)
      psi.estJR <- cbind(psi.estJR,c(temp$psi.hat[2,1:2],temp$psi.hat[1,]))   
      g.2 <- c(g.2,unique((cbind(1,L2,A1,L2*A1)%*%temp$psi.hat[1,])[L2==1 & A1==1]))
#      g.2 <- c(g.2,unique((cbind(1,L2,A1,L2*A1)%*%temp$psi.hat[1,])))
}
print(paste("n = ",n,sep=""))

postscript("dens_g2_n250_noEL.eps")
hist(g.2*(g.2>0),main="",xlab="",freq=F)
abline(v=0,lwd=3)
abline(v=mean(g.2*(g.2>0)),lwd=3,lty=2)
lines(density(g.2*(g.2>0)))
dev.off()

postscript("dens_psi10_n250_noEL.eps")
hist(psi.estJR[1,],main="",xlab="",freq=F,ylim=c(0,2.5))
abline(v=1.6,lwd=3)
abline(v=mean(psi.estJR[1,]),lwd=3,lty=2)
lines(density(psi.estJR[1,]))
dev.off()

postscript("dens_psi11_n250_noEL.eps")
hist(psi.estJR[2,],main="",xlab="",freq=F,ylim=c(0,9))
abline(v=0,lwd=3)
abline(v=mean(psi.estJR[2,]),lwd=3,lty=2)
lines(density(psi.estJR[2,]))
dev.off()