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("pi",run,".dat",sep="") 
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  
  pi <- 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("pi",run,".dat",sep="")
      tmp <- scan(filename)
      rows <- tmp[1]
      cols <- tmp[2]  
      pii <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
      pi <- cbind(pi,pii)
    }
  }

  idx <- seq(1,ncol(pi),skip)
  pi <-pi[idx]
}

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);
    }
  }

  idx <- seq(1,ncol(post),skip);
  post <- post[,idx];

  post <- rbind(post,pi)

  rows <- rows + 1

  rm(tmp)

  tpost <- t(post);

  theta <- c("g","R11","R12","R22","phi","delta","gamma","pi")
  prior.lo <- c(-0.005,  0.00, -1.00,  0.00, -1.00, 0.90,  0.25, -500)
  prior.hi <- c( 0.005,  1.00,  1.00,  1.00,  1.00, 1.05, 10.00,  500)
  
  run <- paste("0",j,sep="") 

  dimnames(tpost) <- list(NULL,theta);
  tpost <- as.data.frame(tpost)
  print("  ")
  print(mean(tpost))
  print("  ")
  print(cor(tpost))
  
  print("  ")
  print("Minima and Maxima")
  for (j in (1:7)) {
    print(paste(theta[j],"   ",prior.lo[j],"  ",min(tpost[,j]),"  ",
            max(tpost[,j]),"  ",prior.hi[j],sep=""))
  }
  print("  ")

  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[,j],lag.max=alag,type="correlation",plot=F,demean=T)
    y <- aut$acf
    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(1,1),mar=c(4.5,4.0,1.5,1)+0.1,pty="s",xaxt="n",yaxt="n")
  pairs(tpost)
  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 variance","stock mean","stock variance") 
  prior.lo <- c( 0.00, 0,     0.00, 0)
  prior.hi <- c( 6.15, 0.025, 6.15, 0.25)
  
  spost <- t(post)
  dimnames(spost) <- list(NULL,stats);
  spost <- as.data.frame(spost)
  print("  ")
  print(mean(spost))
  print("  ")
  print(cor(spost))
  print("  ")
  
  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="")
      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(1,1),mar=c(4.5,4.0,1.5,1)+0.1,pty="s",xaxt="n",yaxt="n")      
# pairs(spost      
# 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() 

}