######################################################################
#
# 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 of two kinds:
#
#   Model 1: e(t) ~ N(0,sig2)
#   Model 2: e(t) ~ t(0,tau2,nu) where tau2 = ((nu-2)/nu)*sig2
#
# sig2 and nu are known quantities.
#
# Prior specification:
#
#        phi ~ MN(phi0,V0)  
#
######################################################################
library(mvtnorm)
#rmvnorm(n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)))
#rmvt(n, sigma = diag(2), df = 1, delta = rep(0, nrow(sigma))) 
#dmvnorm(x, mean = rep(0, p), sigma = diag(p), log = FALSE) 
#dmvt(x, delta = rep(0, p), sigma = diag(p), df=df)


# Simulating from an Gaussian AR(2) process
set.seed(1245)
gaussian.data = FALSE
p    = 2
sig  = 0.1
nu   = 4
tau  = sqrt((nu-2)/nu)*sig
phi0 = 0.00
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)
if (gaussian.data){
  for (t in 3:(n+100))
    y[t] = sum(c(1,y[t-1],y[t-2])*phi.true) + sig*rnorm(1)
}else{
  for (t in 3:(n+100))
    y[t] = sum(c(1,y[t-1],y[t-2])*phi.true) + tau*rt(1,df=nu)
}

y = y[101:(n+100)]
par(mfrow=c(1,2))
plot(y,type="b")
acf(y)

# Likelihood-based inference (Gaussian model)
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("Estimate","Std. Error","t value","Pr(>|t|)")
tab

# Bayesian inference via SIR 
# The proposal density is a Student's t with 
# df0 degrees of freedom, location phi.mle
# and scale fac0*var.phi.mle
loglike.g = function(phi){sum(dnorm((Y-X%*%phi)/sig,log=TRUE))}
loglike.t = function(phi){sum(dt((Y-X%*%phi)/tau,df=nu,log=T))-(n-p)*log(tau)}

# Prior hyperparameters
phi0 = rep(0,q)
V0   = diag(4,q)

# Implementing SIR
df0  = 5 
fac0 = 4 
W    = fac0*var.phi.mle
M    = 100000

# Step 1: Draws from proposal distribution
set.seed(1235)
phi.draw = rmvt(M,delta=phi.mle,sigma=W,df=df0)

# Step 2: Evaluating weights
logw.g = rep(0,M)
logw.t = rep(0,M)
for (i in 1:M){
  llike.g   = loglike.g(phi.draw[i,])
  llike.t   = loglike.t(phi.draw[i,])  
  lprop     = dmvt(phi.draw[i,],delta=phi.mle,sigma=W,df=df0,log=TRUE)
  lprior    = dmvnorm(phi.draw[i,],mean=phi0,sigma=V0,log=TRUE)
  logw.g[i] = llike.g+lprior-lprop
  logw.t[i] = llike.t+lprior-lprop  
}
w.g = exp(logw.g)
pred.g = sum(w.g)
w.t = exp(logw.t)
pred.t = sum(w.t)
w.g = exp(logw.g-max(logw.g))
w.t = exp(logw.t-max(logw.t))

# Step 3: Resampling
i1 = sample(1:M,size=M,replace=TRUE,prob=w.g)
i2 = sample(1:M,size=M,replace=TRUE,prob=w.t)
phi.draws = array(0,c(2,M,q))
phi.draws[1,,] = phi.draw[i1,]
phi.draws[2,,] = phi.draw[i2,]

# Posterior plots
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("Bayes (Normal model)","Bayes (t model)")
par(mfrow=c(2,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)
  points(phi.mle[1],phi.mle[2],col=3,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)
  points(phi.mle[1],phi.mle[3],col=3,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)
  points(phi.mle[2],phi.mle[3],col=3,pch=16,cex=2)
}

par(mfrow=c(2,2))
for (i in 1:3){
  range = range(phi.draws[,,i])
  plot(density(phi.draws[2,,i]),main="",xlab="",col=2,xlim=range)
  lines(density(phi.draws[1,,i]))
  points(phi.true[i],0,col=3,pch=16)
  points(phi.mle[i],0,col=4,pch=16)
  title(paste("phi(",i,")",sep=""))
  if (i==1){
    legend("topleft",legend=name.prior,col=1:2,lty=1,bty="n")
    legend("topright",legend=c("TRUE","MLE (Normal model)"),
           col=3:4,pch=16,bty="n")
  }
}
range = range(phi.draws[,,2]+phi.draws[,,3])
plot(density(phi.draws[2,,2]+phi.draws[2,,3]),main="",xlab="",col=2,xlim=range)
lines(density(phi.draws[1,,2]+phi.draws[1,,3]))
points(phi.true[2]+phi.true[3],0,col=3,pch=16)
points(phi.mle[2]+phi.mle[3],0,col=4,pch=16)
title("phi(2)+phi(3)")


# Model comparison
BF = pred.g/pred.t
c(pred.g,pred.t,round(BF,2))