######################################################################
#
# Bayesian inference in the AR(1) process with known variance:
#
# For t=1,...,n
#    
#              y(t)|y(t-1),phi ~ N(phi*y(t-1),sig2)
#
# We entertain two prior specifications:
#
# Prior 1: phi ~ N(phi0,V0)  - Closed-form posterior inference
# Prior 2: phi ~ t(4)        - Approximate posterior inference via MC 
#
######################################################################

# Simulating from an AR(1) process
set.seed(1245)
sig=1
phi=0.8
n = 60
y = rep(0,n)
for (t in 2:n)
  y[t] = phi*y[t-1]+rnorm(1,0,sig)
ts.plot(y)

# Above this line we only know the observed y1,...,yn and sig

# Likelihood function resembles a N(phi.hat,tau2)
tau2    = sig^2/sum(y[1:(n-1)]^2)
phi.hat = sum(y[2:n]*y[1:(n-1)])/sum(y[1:(n-1)]^2)
c(phi.hat,sqrt(tau2))


############################
# Prior 1) phi ~ N(phi0,V0)
############################
phi0 = 0.7
V0   = 0.025

# Closed form posterior updates
# Combining the prior and the likelihood
V1 = 1/(1/V0+1/tau2)
phi1 = V1*(phi0/V0+phi.hat/tau2)
c(phi1,sqrt(V1))

# Plotting prior, likelihood and posterior
phis = seq(0.3,1.2,length=1000)
plot(phis,dnorm(phis,phi0,sqrt(V0)),type="l",
     xlab="phi",ylab="Density",
ylim=c(0,10))
lines(phis,dnorm(phis,phi.hat,sqrt(tau2)),col=2)
lines(phis,dnorm(phis,phi1,sqrt(V1)),col=3)
abline(v=-1,lty=2)
abline(v=1,lty=2)
legend("topleft",legend=c("Prior","Likelihood",
       "Posterior"),col=1:3,lty=1)
tail = 1-pnorm(1,phi1,sqrt(V1))
title(paste("Pr(phi>1|data)=",round(100*tail,2),"%",sep=""))


############################
# Prior 2: phi ~ t(4)
############################
g = function(phi){
  dnorm(phi,phi.hat,sqrt(tau2))*dt(phi,df=nu)
}
nu = 4

# Implementing sampling importance resampling (SIR)
# -------------------------------------------------

# Step 1) sample M from proposal distribution U(-10,10)
M = 1000000
phi.draw = runif(M,-10,10)
hist(phi.draw)

# Step 2) compute reweights
w = g(phi.draw)
w = w/sum(w)

# Step 3) resampling
phi.post = sample(phi.draw,size=M,replace=TRUE,prob=w)

# Prior, likelihood and posterior 
hist(phi.post,prob=TRUE,xlab="phi",main="")
lines(phis,dnorm(phis,phi.hat,sqrt(tau2)),col=2)
lines(phis,dt(phis,df=nu),col=3)

# Comparing both posteriors
plot(density(phi.post),xlab="phi",main="")
lines(phis,dnorm(phis,phi1,sqrt(V1)),col=2)
legend("topleft",legend=c("Posterior (prior Gaussian)","Posterior (prior Student's t(4)"),col=2:1,lty=1)

# Comparing Pr(phi>1|data) under both priors
tail.gp = 1-pnorm(1,phi1,sqrt(V1))
tail.tp = mean(phi.post>1)
c(tail.gp,tail.tp,tail.tp/tail.gp)