###############################################################################################
#
#  Bayesian Analysis for the AR(1) model with Markov Switching 
#
###############################################################################################
#
#  Author: Hedibert Freitas Lopes
#  Affiliation: Graduate School of Business, University of Chicago
#
#  This version: May, 13th 2005
#
###############################################################################################
#
#  A few steps before you "source" this file into R:
#
#  [1] In R change the directory to the one you currently have the executable file
#
#                  "markovswitchingAR1.exe"
#
#
#  [2] Change the simulation specification accordinlgy.
#
#
#  [3] Once [1] and [2] are done, simply source this file:
#
#  > source("markovswitching.R")
#
#
#  [4] Several files will be created plus graphical summaries in a PDF file called
#
#                  "markovswitching.pdf"
#
###############################################################################################
set.seed(93726)

# Specification
n     = 401 
k     = 2
alpha = c(0.2,1)
phi   = 0.5
sig2  = 0.5
sig   = sqrt(sig2)
P     = matrix(c(0.99,0.01,0.01,0.99),2,2,byrow=T)
theta = alpha 
tau   = 100 
kappa = phi 
eta   = 1 
par1  = 5 
par2  = 0.15
u     = matrix(1,k,k)
M0    = 2000 
M     = 2000 
LAG   = 1  


# Data simulation 
# ---------------
S     =  rep(1,n)
y     = rep(alpha[S[1]]/(1-phi),n) 
for(i in 2:n){
  S[i] = sample(1:k,1,prob=P[S[i-1],])
  y[i] = rnorm(1,alpha[S[i]]+phi*y[i-1],sig)}

pdf(file="markovswitching.pdf")
par(mfrow=c(1,2)) 
ts.plot(y,xlab="time",ylab="",main="Observed time series")
ts.plot(S,xlab="time",ylab="",main="Regime latent indicator")


# Creating parameters.dat, which will be called by markovswitchingAR1.exe
# -----------------------------------------------------------------------
write(c(n,k,M0,M,LAG),"parameters.dat",ncolumns=5) 
write(c(alpha,phi,sig2),"parameters.dat",ncolumns=k+2,append=T) 
write(c(theta,kappa,par1),"parameters.dat",ncolumns=k+3,append=T) 
write(c(tau,eta,par2),"parameters.dat",ncolumns=4,append=T)
for (i in 1:k)   
  write(P[i,],"parameters.dat",ncolumns=k,append=T) 
for (i in 1:k) 
  write(u[i,],"parameters.dat",ncolumns=k,append=T) 
write(y,"Y.dat",ncolumns=1)
write(S,"S.dat",ncolumns=1)

# Calling Fortran 90 executable
# -----------------------------
system.time(system("markovswitchingAR1.exe"))

# Reading the MCMC output
# -----------------------
pars = matrix(scan("pars.out"),M,k+2,byrow=T)
mS   = matrix(scan("mS.out"),n,k,byrow=T)
Ps   = matrix(scan("Ps.out"),M,(k-1)*k,byrow=T)
par(mfrow=c(k,k))
ts.plot(y,xlab="time",ylab="",main="y(t)")
plot(S,xlab="time",ylab="S",main="True S x 1+Pr(S(i)=2|data)") 
lines(1+mS[,2],col=2)
l = 0
for (i in 1:k)
  for (j in 1:(k-1)){
    l = l + 1
    ts.plot(Ps[,l],xlab="",ylab="",main=paste("P",i,j,sep=""))
    abline(h=P[i,j],col=2)
  }

par(mfrow=c(2,2)) 
for (i in 1:k){
  ts.plot(pars[,i],xlab="",ylab="",main=paste("alpha",i,sep=""),ylim=range(c(pars[,i],alpha[i])))
  abline(h=alpha[i],col=2) 
} 
ts.plot(pars[,k+1],xlab="",ylab="",main="phi",ylim=range(pars[,k+1],phi))
abline(h=phi,col=2)
ts.plot(pars[,k+2],xlab="",ylab="",main="sig2",ylim=range(pars[,k+2],sig2))
abline(h=sig2,col=2)
dev.off()