########################################################################
# 
# BAYESIAN ESTIMATION OF THE AR(P) MODEL
#
########################################################################
rm(list=ls())
set.seed(98765)

ns   = c(120,600,1200)
beta = c(0.987,0.974,-0.098,-0.051,-0.155,0.045,
        -0.0702,-0.036,0.179,0.025,0.138,-0.288)
sig2 = 0.025
sig  = sqrt(sig2)
nsim = 16
pmax = 20

pdf(file="ARp-bayesianselection.pdf",width=15,height=10)

for (n in ns){
  bestmodel    = rep(0,pmax)
  bestmodelAIC = rep(0,pmax)
  bestmodelBIC = rep(0,pmax)
  index        = 1:pmax
  par(mfrow=c(4,4))
  for (sim in 1:nsim){
    y = rep(0,100*n)
    for (t in 12:(100*n))
      y[t] = rnorm(1,beta[1]+sum(beta[2:12]*y[(t-1):(t-11)]),sig)
    y = y[(99*n+1):(100*n)]  
    lpred = rep(0,pmax)
    AIC   = rep(0,pmax)
    BIC   = rep(0,pmax)
    for (p in 1:pmax){
      Y = y[(p+1):n]
      X = rep(1,n-p)
      for (j in 1:p){
        X = cbind(X,y[(p+1-j):(n-j)]) 
      }
      n1       = nrow(X)
      In       = diag(1,n1)
      d        = ncol(X)
      b.ols    = solve(t(X)%*%X)%*%t(X)%*%Y
      error    = Y - X%*%b.ols
      s.ols    = sqrt(mean(error^2))
      AIC[p]   = -2*sum(dnorm(error,0,s.ols,log=TRUE))+2*(d+1) 
      BIC[p]   = -2*sum(dnorm(error,0,s.ols,log=TRUE))+log(n1)*(d+1)
      b0       = rep(0,d)
      B0       = diag(10,d)
      s20      = 0.025
      nu0      = 6
      par1     = (nu0+n1)/2
      B1       = solve(t(X)%*%X+solve(B0))
      b1       = B1%*%(t(X)%*%Y+solve(B0)%*%b0)
      par2     = (nu0*s20+t(y)%*%y+t(b0)%*%solve(B0)%*%b0-t(b1)%*%solve(B1)%*%b1)/2
      P        = (In-X%*%solve(solve(B0)+t(X)%*%X)%*%t(X))/s20
      logdet   = n1*log(s20)+log(det(solve(B0)+t(X)%*%X))+log(det(B0))
      lpred[p] = -0.5*logdet-0.5*(nu0+n1)*log(1+t(Y-X%*%b0)%*%P%*%(Y-X%*%b0)/nu0)
    }
    pred = exp(lpred-max(lpred))
    pred = pred/sum(pred)
    bestmodel[index[pred==max(pred)]]  = bestmodel[index[pred==max(pred)]]  + 1
    bestmodelAIC[index[AIC==min(AIC)]] = bestmodelAIC[index[AIC==min(AIC)]] + 1
    bestmodelBIC[index[BIC==min(BIC)]] = bestmodelBIC[index[BIC==min(BIC)]] + 1
    plot(pred,xlab="Lag",ylab="PMP",type="h",lwd=5,main="",axes=FALSE)
    axis(2);axis(1,at=1:pmax);box()
    title(paste("dataset ",sim,"\n n=",n,sep=""))
  }
  par(mfrow=c(1,3))
  plot(bestmodel/nsim,xlab="LAG",type="h",lwd=5,ylab="Frequency",axes=FALSE)
  axis(2);axis(1,at=1:pmax);box()
  title(paste("BAYES \n n=",n,sep=""))
  plot(bestmodelAIC/nsim,xlab="LAG",type="h",lwd=5,ylab="Frequency",axes=FALSE)
  axis(2);axis(1,at=1:pmax);box()
  title(paste("AIC \n n=",n,sep=""))
  plot(bestmodelBIC/nsim,xlab="LAG",type="h",lwd=5,ylab="Frequency",axes=FALSE)
  axis(2);axis(1,at=1:pmax);box()
  title(paste("BIC \n n=",n,sep=""))
}
dev.off()