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

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

for(n in n.sample) {
psi<-c(psi10,psi11,psi20,psi21,psi22,psi23)
results <- psi
   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,])) 

      start.search <- c(temp$psi.hat[2,1:2],temp$psi.hat[1,])
      temp2 <- nlm(AAOg_estSeqNS.fcn,p=start.search,Y=Y,action=A,status=L)
      psi.estAAO <- cbind(psi.estAAO,temp2$estimate)
   }
   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

   t1a <- apply(psi.estAAO,1,mean)
   t2a <- apply(sqrt(var.estJR),1,mean) # calculated se FROM RECURSIVE METHOD!
   t3a <- apply(sqrt(var.estJR + (psi.estAAO - psi)^2),1,mean)  # rMSE
   t4a <- apply((psi.estAAO - 1.96*sqrt(var.estJR) < psi) & (psi.estAAO + 1.96*sqrt(var.estJR) > psi),1,sum)/n.sim*100 # Wald CI coverage

    results <- cbind(results,cbind(t1g,t2g,t3g,t4g,t1a,t2a,t3a,t4a))
    print(round(results,3))
}
n
apply((psi.estAAO-psi.estJR)[1:2,],1,summary)

par(mfrow=c(2,3))
apply(psi.estAAO,1,hist,xlab="psi",main="")
apply(psi.estAAO,1,median)
print(round(results,3))

######################################################################################################

## start.search <- c(temp$psi.hat[2,1:2],temp$psi.hat[1,])
start.search <- c(1.6,0,-3,2,1,0)
temp <- nlm(AAOg_estSeqNS.fcn,p=start.search,Y=Y,action=A,status=L)
psi.hat <- temp$estimate
#psi.min <- temp$estimate
psi.min <- psi.estAAO[,min.val==min(min.val)] 

par(mfrow=c(2,3))
  temp <- NULL
  for(i in seq(psi.hat[1]-1,psi.hat[1]+1,.05)) 
     temp <- cbind(temp,AAOg_estSeqNS.fcn(PSI=c(i,psi.hat[2:6]),Y,L,A,return.all=T))
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[1,],type="l",xlab="psi10",ylab="U(1)")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[2,],type="l",xlab="psi10",ylab="U(2)")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[3,],type="l",xlab="psi10",ylab="U(3)")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[4,],type="l",xlab="psi10",ylab="U(4)")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[5,],type="l",xlab="psi10",ylab="U(5)")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[6,],type="l",xlab="psi10",ylab="U(6)")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)
#dev.off()

  temp <- NULL
  for(i in seq(psi.hat[2]-1,psi.hat[2]+1,.05)) 
     temp <- cbind(temp,AAOg_estSeqNS.fcn(PSI=c(psi.hat[1],i,psi.hat[3:6]),Y,L,A,return.all=T))
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[1,],type="l",xlab="psi11",ylab="U(1)")
  abline(v=psi.min[2],lty=2)
  abline(v=psi11,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[2,],type="l",xlab="psi11",ylab="U(2)")
  abline(v=psi.min[2],lty=2)
  abline(v=psi11,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[3,],type="l",xlab="psi11",ylab="U(3)")
  abline(v=psi.min[2],lty=2)
  abline(v=psi11,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[4,],type="l",xlab="psi11",ylab="U(4)")
  abline(v=psi.min[2],lty=2)
  abline(v=psi11,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[5,],type="l",xlab="psi11",ylab="U(5)")
  abline(v=psi.min[2],lty=2)
  abline(v=psi11,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[6,],type="l",xlab="psi11",ylab="U(6)")
  abline(v=psi.min[2],lty=2)
  abline(v=psi11,lwd=2)
  abline(h=0)

  temp <- NULL
  for(i in seq(psi.hat[3]-1,psi.hat[3]+1,.05))
     temp <- cbind(temp,AAOg_estSeqNS.fcn(PSI=c(psi.hat[1:2],i,psi.hat[4:6]),Y,L,A,return.all=T))
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[1,],type="l",xlab="psi20",ylab="U(1)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi20,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[2,],type="l",xlab="psi20",ylab="U(2)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi20,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[3,],type="l",xlab="psi20",ylab="U(3)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi20,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[4,],type="l",xlab="psi20",ylab="U(4)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi20,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[5,],type="l",xlab="psi20",ylab="U(5)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi20,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[6,],type="l",xlab="psi20",ylab="U(6)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi20,lwd=2)
  abline(h=0)

### ETC