#simulate data set.seed(1) source("AR_PS.R") #simulate AR process n <- 500 t <- c(160,240,360,500) alpha<-matrix(0,nrow=4,ncol=3) alpha[1,]<-c(-0.6,0,0) alpha[2,] <- c(-1.2,-0.36,0) alpha[3,] <- c(-1.2,-0.36,0.25) alpha[4,] <- c(-0.9,0,0) sig<-c(1.2,0.8,1.3,1) y <- rep(0,n+3) y[1:3] <- rnorm(3,0,2) j <- 1 for(i in 1:n){ if(i>t[j]) j <- j+1 y[3+i] <- sum(y[2:0+i]*alpha[j,])+rnorm(1,0,sig[j]) } #load in functions from AR_PS.R #analyse under a model that allows AR models of up to log 3; uniform prior # on these lags For<-AR.forward(y,mod=3,delta=c(1,0.4,0.15),lambda=0.005,gamma=2,pi=rep(1/3,3)) #simulate 10 realisations of the changepoints cpts.sim <- ARcpts.sim(For,10,mod=3,pi=rep(1/3,3),lambda=0.005) #MAP estimate conditional on 3 segments cMAP <- AR.MAP(For,mod=3,pi=rep(1/3,3),lambda=0.005,mhat=5) #unconditional MAP estimate uMAP <- AR.MAP(For,mod=3,pi=rep(1/3,3),lambda=0.005)