######################################################################
#
# Bayesian inference in the AR(2) process with known variance:
#
# For t=1,...,n
#    
#     y(t) = phi0+phi1*y(t-1)+phi2*y(t-2) + e(t) 
#
# with iid errors e(t) ~ N(0,sig2)
#
# We entertain two prior specifications for phi=(phi0,phi1,phi2)
#
# Prior 1: phi ~ MN(phi0,V0)  
# Prior 2: phi ~ Mt(nu0,phi0,C0), where nu0>2 and C0=(nu0-2)/nu0*V0
#
# There phi's under both priors have the same first two moments, 
# but Prior 2 has heavier tails than Prior1
# 
######################################################################
library(mvtnorm)
#dmvnorm(x, mean = rep(0, p), sigma = diag(p), log = FALSE, checkSymmetry = TRUE) 
#rmvnorm(n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)), method=c("eigen", "svd", "chol"), pre0.9_9994 = FALSE, checkSymmetry = TRUE)
#rmvt(n, sigma = diag(2), df = 1, delta = rep(0, nrow(sigma)), type = c("shifted", "Kshirsagar"), ...) 
#dmvt(x, delta = rep(0, p), sigma = diag(p), df


# Simulating from an AR(2) process
set.seed(1245)
p    = 2
sig  = 1
phi0 = 0.00
#phi1 = 1.78
#phi2 = -0.8
phi1 = 0.60
phi2 = 0.38
phi.true = c(phi0,phi1,phi2)
n    = 100
sig2 = sig^2
q = p + 1
y = rep(0,n+100)
for (t in 3:(n+100))
  y[t] = sum(c(1,y[t-1],y[t-2])*phi.true)+rnorm(1,0,sig)
y = y[101:(n+100)]
par(mfrow=c(1,2))
plot(y,type="b")
acf(y)

# Likelihood-based inference
Y           = y[3:n]
X           = cbind(1,y[2:(n-1)],y[1:(n-2)])
phi.mle     = solve(t(X)%*%X)%*%t(X)%*%Y
var.phi.mle = sig2*solve(t(X)%*%X)
error.mle   = sqrt(diag(var.phi.mle))
t.stat      = phi.mle/error.mle
pvalue      = 2*(1-pnorm(t.stat))

# Estimation summaries
tab = cbind(round(phi.mle,5),round(error.mle,5),
            round(t.stat,3),round(pvalue,6))
rownames(tab) = c("phi0","phi1","phi2")
colnames(tab) = c("Estimat","Std. Error","t value","Pr(>|t|)")
tab

# Compare to the output from R function "lm"
summary(lm(Y~X-1))

###################################
# Prior 1) phi(i) ~ N(phi0,V0)
# Prior 2) phi(i) ~ t(nu0,phi0,C0)
###################################
phi0 = rep(0,q)
V0   = diag(4,q)
nu0  = 4 
C0   = (nu0-2)/nu0*V0



# Posterior inference via SIR 
# The proposal density is Student's t with 
# df0 degrees of freedom, location phi.mle
# and scale fac0*var.phi.mle
loglike = function(phi){
  -sum((Y-X%*%phi)^2)/(2*sig2)
}
df0  = 5 
fac0 = 4 
W    = fac0*var.phi.mle
M    = 10000

set.seed(1235)
phi.draw = rmvt(M,delta=phi.mle,sigma=W,df=df0)
logw1 = rep(0,M)
logw2 = rep(0,M)
for (i in 1:M){
  llike   = loglike(phi.draw[i,])
  lprop   = dmvt(phi.draw[i,],delta=phi.mle,sigma=W,df=nu0,log=TRUE)
  lprior1 = dmvnorm(phi.draw[i,],mean=phi0,sigma=V0,log=TRUE)
  lprior2 = dmvt(phi.draw[i,],delta=phi0,sigma=C0,df=nu0,log=TRUE)
  logw1[i] = llike+lprior1-lprop
  logw2[i] = llike+lprior2-lprop  
}
w1 = exp(logw1-max(logw1))
w2 = exp(logw2-max(logw2))
i1 = sample(1:M,size=M,replace=TRUE,prob=w1)
i2 = sample(1:M,size=M,replace=TRUE,prob=w2)
phi.draws = array(0,c(2,M,q))
phi.draws[1,,] = phi.draw[i1,]
phi.draws[2,,] = phi.draw[i2,]


min = max = rep(0,q)
for (i in 1:q){
  min[i] = min(apply(phi.draws[,,i],1,min))
  max[i] = max(apply(phi.draws[,,i],1,max))
}
name.prior = c("Normal","Student's t")
par(mfrow=c(3,q))
for (i in 1:2){
  plot(phi.draws[i,,1],phi.draws[i,,2],col=gray(0.75),pch=16,
  xlim=c(min[1],max[1]),ylim=c(min[2],max[2]),xlab="phi(1)",ylab="phi(2)")
  points(phi.true[1],phi.true[2],col=2,pch=16,cex=2)
  title(name.prior[i])
  plot(phi.draws[i,,1],phi.draws[i,,3],col=gray(0.75),pch=16,
  xlim=c(min[1],max[1]),ylim=c(min[3],max[3]),xlab="phi(1)",ylab="phi(3)")
  points(phi.true[1],phi.true[3],col=2,pch=16,cex=2)
  plot(phi.draws[i,,2],phi.draws[i,,3],col=gray(0.75),pch=16,
  xlim=c(min[2],max[2]),ylim=c(min[3],max[3]),xlab="phi(2)",ylab="phi(3)")
  points(phi.true[2],phi.true[3],col=2,pch=16,cex=2)
}
for (i in 1:3){
  plot(density(phi.draws[1,,i]),main="",xlab="")
  lines(density(phi.draws[2,,i]),col=2)
  title(paste("phi(",i,")",sep=""))
  if (i==2){
    legend("topleft",legend=name.prior,col=1:2,lty=1,bty="n")
  }
}