########################################################################
library("mvtnorm")
set.seed(12345)
phi  = c(0.95,0.5)
phi  = c(0.99,0.8)
sig  = 3
n    = 100
tau0 = 50
y    = rep(0,n)
for (t in 2:tau0)
  y[t] = phi[1]*y[t-1] + rnorm(1,0,sig)	
mu = y[tau0]*(1-phi[2])  
for (t in (tau0+1):n)
  y[t] = mu + phi[2]*y[t-1] + rnorm(1,0,sig)	  
y = y - mu/(1-phi[2])
########################################################################

y = c(-25.9,-24.2,-22.0,-22.4,-23.8,-22.0,-27.5,-25.6,-26.4,-27.3,-30.0,-30.3,
-24.8,-23.7,-22.2,-24.5,-22.0,-24.7,-25.7,-22.4,-21.5,-19.2,-14.9,-17.0,-21.7,
-26.6,-21.1,-22.6,-20.8,-19.0,-19.6,-17.2,-10.7,-4.7,0.0,0.5,1.7,0.5,-4.8,0.3,
0.1,3.2,-4.2,-7.6,-5.0,-2.6,1.5,-3.0,-1.5,0.0,-3.9,-4.8,2.0,1.8,2.5,0.0,0.8,
2.7,4.7,10.2,1.1,1.3,-3.0,-0.7,4.2,1.6,-4.2,-0.7,4.2,4.9,0.0,0.2,-2.2,-4.9,3.1,
6.7,8.2,9.0,4.8,5.2,7.3,7.7,9.3,6.5,12.7,13.0,16.0,14.8,11.0,10.4,10.8,5.7,2.0,
7.3,4.6,0.8,2.7,0.6,7.0,3.8)

n = length(y)

par(mfrow=c(1,1))
ts.plot(y)

#######################################################
# (1) MLE for various values of tau
#######################################################
taus = 25:45
ntau = length(taus)
phis = matrix(0,ntau,2)
sigs = rep(0,ntau)
for (i in 1:ntau){
  y1 = y[2:taus[i]]
  X1 = y[1:(taus[i]-1)]
  y2 = y[(taus[i]+1):n]
  X2 = y[taus[i]:(n-1)]
  phis[i,1] = sum(X1*y1)/sum(X1^2)
  phis[i,2] = sum(X2*y2)/sum(X2^2)
  sigs[i]   = sqrt(mean(c((y1-phis[i,1]*X1)^2,(y2-phis[i,2]*X2)^2)))
}

par(mfrow=c(1,2))
plot(taus,phis[,1],type="b",ylim=range(phis),xlab=expression(tau),ylab=expression(phi))
points(taus,phis[,2],type="b",col=2)
legend("right",legend=c(expression(phi[1]),expression(phi[2])),col=1:2,lty=1)
plot(taus,sigs,type="h",xlab=expression(tau),ylab=expression(sigma))


#######################################################
# (2)-(3) Gibbs sampler and posterior inference
#######################################################
# Prior hyperparameters
m0 = 0
C0 = 1
a0 = 5
b0 = 6
tl = 15
tu = 85


# Initial value
sig2 = 1
t0   = 35


# Gibbs sampler
M0    = 10000
M     = 2000
L     = 10
niter = M0+L*M
draws = matrix(0,niter,4)
for (iter in 1:niter){
  # sampling phi1
  var  = 1/(1/C0+sum(y[1:(t0-1)]^2)/sig2)
  mean = var*(m0/C0+sum(y[1:(t0-1)]*y[2:t0])/sig2)
  phi1 = rnorm(1,mean,sqrt(var))

  # sampling phi2
  var  = 1/(1/C0+sum(y[t0:(n-1)]^2)/sig2)
  mean = var*(m0/C0+sum(y[t0:(n-1)]*y[(t0+1):n])/sig2)
  phi2 = rnorm(1,mean,sqrt(var))

  # Sampling sig2
  phis = c(rep(phi1,t0-1),rep(phi2,n-t0))
  par1 = a0+(n-1)/2
  par2 = b0+sum((y[2:n]-phis*y[1:(n-1)])^2)/2
  sig2 = 1/rgamma(1,par1,par2)

  # sampling t0
  w = rep(0,tu-tl+1)
  for (tt in tl:tu){
    phis = c(rep(phi1,tt-1),rep(phi2,n-tt))
    mean = y[2:n]-phis*y[1:(n-1)]
    w[tt-tl+1] = sum(dnorm(n-1,mean,sqrt(sig2),log=TRUE))
  }
  w = exp(w-max(w))
  w = w/sum(w)
  t0 = sample(tl:tu,size=1,prob=w)
  
  # storing the draws
  draws[iter,] = c(phi1,phi2,sig2,t0)
}
draws = draws[seq(M0+1,niter,by=L),]

# Posterior summaries
par(mfrow=c(2,3))
ts.plot(draws[,1])
acf(draws[,1])
hist(draws[,1])
ts.plot(draws[,2])
acf(draws[,2])
hist(draws[,2])

par(mfrow=c(2,3))
ts.plot(draws[,3])
acf(draws[,2])
hist(draws[,2])
plot(draws[,4])
acf(draws[,4])
plot(table(draws[,4])/M)

par(mfrow=c(1,1))
ts.plot(y)
abline(v=mean(draws[,4]),col=4)

tab = apply(draws,2,summary)
colnames(tab) = c("phi1","phi2","sig2","tau")
round(tab,3)