####################################################################
#  Simulated data
####################################################################

#pdf(file="AR1-with-changingpoint.pdf",width=10,height=8)
set.seed(31416)
n      = 200
tau   = 100
phi1 = 0.2
phi2 = 0.9
sig   = 0.2
sig2 = sig^2
y = rep(0,2*n)
for (t in 2:(n+tau))
  y[t] = rnorm(1,phi1*y[t-1],sig)
for (t in (n+tau+1):(2*n))
  y[t] = rnorm(1,phi2*y[t-1],sig)
y = y[(n+1):(2*n)]

ts.plot(y)
abline(v=tau,lty=2)


# Gibbs sampler
M0 = 10000
M = 10000
niter = M0+M
draws = matrix(0,niter,3)
taus = tau
for (iter in 1:niter){
y1 = y[1:taus]
y2 = y[(taus+1):n]
n1 = taus
n2 = n-taus
v1 = sig2/sum(y1[1:(n1-1)]^2)
v2 = sig2/sum(y2[1:(n2-1)]^2)
phi1hat = sum(y1[2:n1]*y1[1:(n1-1)])*v1/sig2
phi2hat = sum(y2[2:n2]*y2[1:(n2-1)])*v2/sig2
draw1 = rnorm(1,phi1hat,sqrt(v1))
draw2 = rnorm(1,phi2hat,sqrt(v2))
w = rep(0,n)
for (j in 2:(n-2)){
  n1 = j
  n2 = n-j
  y1 = y[1:n1]
  y2 = y[(n1+1):n]
  w[j] = sum(dnorm(y1[2:n1],draw1*y1[1:(n1-1)],sig,log=TRUE))+
             sum(dnorm(y2[2:n2],draw2*y2[1:(n2-1)],sig,log=TRUE))                   
}
w = w[2:(n-2)]
w = exp(w-max(w))
taus = sample(2:(n-2),size=1,prob=w)
draws[iter,] = c(draw1,draw2,taus)
}
draws = draws[(M0+1):niter,]
count = rep(0,n)
for (i in 1:M)
  count[draws[i,3]]=  count[draws[i,3]]+1

pdf(file="AR1-with-changingpoint-GibbsSampler.pdf",width=10,height=8)
par(mfrow=c(1,1))
ts.plot(y)
abline(v=tau,lty=2)

par(mfrow=c(3,3))
ts.plot(draws[,1],xlab="iterations",ylab="",main=expression(phi[1]))
acf(draws[,1],main="")
hist(draws[,1],prob=TRUE,xlab="",main="")
points(phi1,0,pch=16,col=2)
ts.plot(draws[,2],xlab="iterations",ylab="",main=expression(phi[2]))
acf(draws[,2],main="")
hist(draws[,2],prob=TRUE,xlab="",main="")
points(phi2,0,pch=16,col=2)
ts.plot(draws[,3],xlab="iterations",ylab="",main=expression(tau))
acf(draws[,1],main="")
#hist(draws[,3],xlab="",main="")
plot(1:n,count/M,type="h",xlab=expression(tau),ylab="Proportion")
points(tau,0,pch=16,col=2)
dev.off()