number <- 5
skip <- 5
alag <- 100

library(MASS)
library(ts)

ps.options(horizontal=FALSE,onefile=FALSE)
ps.options(pagecentre=TRUE)
ps.options(paper="letter")
ps.options(width=7.0,height=10.0)
  
list <- c("01","02","03","04","05","06","07","08","09","10","11","12","13")

list <- c("01")
  
for (j in list) { 

  run <- paste("0",j,"000",sep="") 
  
  filename <- paste("theta",run,".dat",sep="") 
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  
  post <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  
  if (number > 0) {
    for (i in (1:number) ) {
      if (i < 10) {
        run <- paste("0",j,"00",i,sep="")
      } else {
        run <- paste("0",j,"0",i,sep="") 
      }
      filename <- paste("theta",run,".dat",sep="")
      tmp <- scan(filename)
      rows <- tmp[1]
      cols <- tmp[2]
      posti <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
      post <- cbind(post,posti);
    }
  }

  rm(tmp)

  theta <- c("g","sig","rho","sig_w","phi","delta","gamma")
  prior.lo <- c(-3.0, 0.0, 0.0,  0.8, 0.5, 0.90, 1.0)
  prior.hi <- c( 5.0,25.0, 1.0, 50.0, 1.0, 0.99, 2.0)

  run <- paste("0",j,sep="") 
  
  idx <- seq(1,ncol(post),skip);
  post <- post[,idx];

  post[1,] <- 12*100*post[1,]

  for ( i in seq(1,ncol(post)) ) {
    R <- array(c(post[2,i],0,post[3,i],post[4,i]),dim=c(2,2))
    S <- R%*%t(R)
    sdev11 <- sqrt(S[1,1])
    sdev22 <- sqrt(S[2,2])
    sdev12 <- S[1,2]/(sdev11*sdev22)
    sig  <- sqrt(12)*100*sdev11
    rho  <- sdev12
    sigw <- sqrt(12)*100*sdev22
    post[2,i] <- sig
    post[3,i] <- rho
    post[4,i] <- sigw
  }

  post[5,] <- post[5,]^12   
     
  post[6,] <- post[6,]^12   
 
  tpost <- t(post);
  dimnames(tpost) <- list(NULL,theta);
  tpost <- as.data.frame(tpost)
  print(mean(tpost))
  print(cor(tpost))
  
  filename <- paste("thser",run,".eps",sep="")  
  postscript(file=filename)
  par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
  for (j in 1:rows)
  {
    plot(idx,post[j,],type="n",ylab="")
      lines(idx,post[j,],lty="solid")
      title(theta[j])
  }
  dev.off()
  
  filename <- paste("thden",run,".eps",sep="") 
  postscript(file=filename)
  par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
  for (j in 1:rows)
  {
    #den <- density(post[j,], from=prior.lo[j], to=prior.hi[j],adjust=4)
    #plot(den,main=theta[j])
    den <- density(post[j,],adjust=4)
    xidx <- ( prior.lo[j] <= den$x ) & ( den$x <= prior.hi[j] )
    x <- den$x[xidx]
    y <- den$y[xidx]
    plot(x,y,type="n",ylab="",yaxt="n")
      lines(x,y,lty="solid")
      title(theta[j])
  }
  dev.off()

  rm(den,x,y,xidx)

  filename <- paste("thaut",run,".eps",sep="")
  postscript(file=filename)
  par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
  for (j in 1:rows)
  {
    aut <- acf(tpost,lag.max=alag,type="correlation",plot=F,demean=T)
    y <- aut$acf[,j,j]
    aidx <- 0:alag
    plot(c(0,alag),c(0,1),type="n",ylab="")
      lines(aidx,y,lty="solid")
      title(theta[j])
  }
  dev.off()

  rm(aut,y,aidx)
 
  filename <- paste("thcor",run,".eps",sep="")  
  postscript(file=filename) 
  par(mfrow=c(rows,rows),mar=c(4.5,4.0,1.5,1)+0.1,pty="s",xaxt="n",yaxt="n")
  for (j in 1:rows)
  {
    for (i in 1:rows)
    {
       plot(x=post[i,],y=post[j,],xlab=theta[i],ylab=theta[j],pch=".") 
    }
  }
  dev.off()

# filename <- paste("thbiv",run,".eps",sep="") 
# postscript(file=filename) 
# par(mfrow=c(rows,rows),mar=c(0.0,0.0,0.0,0.0)+0.1)
# for (j in 1:rows)
# {  
#   for (i in 1:rows) 
#   {  
#      f2 <- kde2d(x=post[i,],y=post[j,],n=70)
#      persp(f2)
#   }
# }
# dev.off() 

}


for (j in list) { 

  run <- paste("0",j,"000",sep="") 
  
  filename <- paste("stats",run,".dat",sep="") 
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  
  post <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  
  if (number > 0) {
    for (i in (1:number) ) {
      if (i < 10) {
        run <- paste("0",j,"00",i,sep="")
      } else {
        run <- paste("0",j,"0",i,sep="") 
      }
      filename <- paste("stats",run,".dat",sep="")
      tmp <- scan(filename)
      rows <- tmp[1]
      cols <- tmp[2]  
      posti <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
      post <- cbind(post,posti);
    }
  }

  rm(tmp)

  run <- paste("0",j,sep="") 
  
  idx <- seq(1,ncol(post),skip);
  post <- post[,idx];
  
  stats <- c("bond mean","bond std dev","stock mean","stock std dev") 
  prior.lo <- c(-1.0,  0.0, -25.0,  0.0)
  prior.hi <- c( 3.0,  5.0,  50.0, 90.0)
  
  post[1,] <- 100*12*(post[1,]-1)

  post[2,] <- 100*sqrt(12*post[2,])

  post[3,] <- 100*12*(post[3,]-1)

  post[4,] <- 100*sqrt(12*post[4,])

  spost <- t(post)
  dimnames(spost) <- list(NULL,stats);
  spost <- as.data.frame(spost)
  print(mean(spost))
  print(cor(spost)) 

  filename <- paste("stser",run,".eps",sep="")  
  postscript(file=filename)
  par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
  for (j in 1:rows)
  {
    plot(idx,post[j,],type="n",ylab="")
      lines(idx,post[j,],lty="solid")
      title(stats[j])
  }
  dev.off()
  
  filename <- paste("stden",run,".eps",sep="") 
  postscript(file=filename)
  par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
  for (j in 1:rows)
  {
    #den <- density(post[j,], from=prior.lo[j], to=prior.hi[j],adjust=4)
    #plot(den,main=stats[j])
    den <- density(post[j,],adjust=4)
    xidx <- ( prior.lo[j] <= den$x ) & ( den$x <= prior.hi[j] )
    x <- den$x[xidx]
    y <- den$y[xidx]
    plot(x,y,type="n",ylab="",yaxt="n")
      lines(x,y,lty="solid")
      title(stats[j])
  }
  dev.off()

  rm(den,x,y,xidx)
 
  filename <- paste("stcor",run,".eps",sep="")  
  postscript(file=filename) 
  par(mfrow=c(rows,rows),mar=c(4.5,4.0,1.5,1)+0.1,pty="s",xaxt="n",yaxt="n")
  for (j in 1:rows)
  {
    for (i in 1:rows)
    {
       plot(x=post[i,],y=post[j,],xlab=stats[i],ylab=stats[j],pch=".") 
    }
  }
  dev.off()

# filename <- paste("stbiv",run,".eps",sep="") 
# postscript(file=filename) 
# par(mfrow=c(rows,rows),mar=c(0.0,0.0,0.0,0.0)+0.1) 
# for (j in 1:rows)
# {  
#   for (i in 1:rows) 
#   {  
#      f2 <- kde2d(x=post[i,],y=post[j,],n=70)
#      persp(f2)
#   }
# }
# dev.off() 

}