################################ 
##    Regret-plotting code    ## 
################################ 

source("AAOgest_txS.R")

n <- 3000
psi20 <- -3	
psi21 <- 2
psi22 <- 1	
psi23 <- 0	
psi10 <- psi20	
psi11 <- psi21
alpha <- c(.5,.5)
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      

psi.hat <- c(-3,2,1,0)
AAOg_estSeq.fcn(PSI=psi.hat,Y,L,A)
AAOg_estSeq.fcn(PSI=psi.hat,Y,L,A,return.all=T)

## NB need to return a univariate value for this to work [max{mean(U^2)}]
#temp <- nlm(AAOg_estSeq.fcn,p=c(rnorm(1,-3,1),rnorm(1,2,1),rnorm(1,1,1),rnorm(1)),Y=Y,action=A,status=L,iterlim=1000) 
#temp <- nlm(AAOg_estSeq.fcn,p=c(-3,2,1,0),Y=Y,action=A,status=L,iterlim=1000) 
#psi.hat <- temp$estimate
#psi.min <- temp$estimate

# for these plots, we return all 4 values [mean(U^2) the vector]
#postscript("AAO_UplotsPSI0.eps")
par(mfrow=c(2,2))
  temp <- NULL
  for(i in seq(psi.hat[1]-1,psi.hat[1]+1,.05)) 
     temp <- cbind(temp,AAOg_estSeq.fcn(PSI=c(i,psi.hat[2:4]),Y,L,A,return.all=T))
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp[1,],type="l",xlab="psi0",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="psi0",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="psi0",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="psi0",ylab="U(4)")
  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_estSeq.fcn(PSI=c(psi.hat[1],i,psi.hat[3:4]),Y,L,A,return.all=T))
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp[1,],type="l",xlab="psi1",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="psi1",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="psi1",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="psi1",ylab="U(4)")
  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_estSeq.fcn(PSI=c(psi.hat[1:2],i,psi.hat[4]),Y,L,A,return.all=T))
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[1,],type="l",xlab="psi22",ylab="U(1)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi22,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[2,],type="l",xlab="psi22",ylab="U(2)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi22,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[3,],type="l",xlab="psi22",ylab="U(3)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi22,lwd=2)
  abline(h=0)
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[4,],type="l",xlab="psi22",ylab="U(4)")
  abline(v=psi.min[3],lty=2)
  abline(v=psi22,lwd=2)
  abline(h=0)

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



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

par(mfrow=c(2,2))
  temp <- NULL
  for(i in seq(psi.hat[1]-1,psi.hat[1]+1,.05)) 
     temp <- cbind(temp,AAOg_estSeq.fcn(PSI=c(i,psi.hat[2:4]),Y,L,A))
  plot(seq(psi.hat[1]-1,psi.hat[1]+1,.05),temp,type="l",xlab="psi0",ylab="max U^2")
  abline(v=psi.min[1],lty=2)
  abline(v=psi10,lwd=2)
  abline(h=0)

  temp <- NULL
  for(i in seq(psi.hat[2]-1,psi.hat[2]+1,.05)) 
     temp <- cbind(temp,AAOg_estSeq.fcn(PSI=c(psi.hat[1],i,psi.hat[3:4]),Y,L,A))
  plot(seq(psi.hat[2]-1,psi.hat[2]+1,.05),temp,type="l",xlab="psi1",ylab="max U^2")
  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_estSeq.fcn(PSI=c(psi.hat[1:2],i,psi.hat[4]),Y,L,A))
  plot(seq(psi.hat[3]-1,psi.hat[3]+1,.05),temp[1,],type="l",xlab="psi22",ylab="max U^2")
  abline(v=psi.min[3],lty=2)
  abline(v=psi22,lwd=2)
  abline(h=0)

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