#################################################################### # Simulated data #################################################################### set.seed(154321) n = 500 tau = 200 phi1 = 0.0 phi2 = 1.0 phi = 0.90 sig = 1 sig2 = sig^2 y = rep(0,2*n) for (t in 2:(n+tau)) y[t] = rnorm(1,phi1+phi*y[t-1],sig) for (t in (n+tau+1):(2*n)) y[t] = rnorm(1,phi2+phi*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 z1 = y1[2:n1]-phi*y1[1:(n1-1)] z2 = y2[2:n2]-phi*y2[1:(n2-1)] y1bar = mean(z1) y2bar = mean(z2) draw1 = rnorm(1,z1,sqrt(sig2/n1)) draw2 = rnorm(1,z2,sqrt(sig2/n2)) 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]-phi*y1[1:(n1-1)],draw1,sig,log=TRUE))+ sum(dnorm(y2[2:n2]-phi*y2[1:(n2-1)],draw2,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-intercept-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()