# Simulate a switching model # rbinom(n, size = 1, prob = p) # Rather use a Markov process for switching between states # First generate Markov chain with transition matrix(c(p1,1-p1,q1,1-q1),ncol=2) # but since we are only in two states this can be avoided markov.chain1 = function(n,p1,q1){ # Burnin is huge, it would have been smarter to find its stationary distribution # I think using at least 1000 burn in is necessary. m = (n+1000) # The Transition Matrix is matrix(c(1-p1,p1,q1,1-q1),ncol=2,byrow=TRUE) # Initial value randomized ini = rbinom(1,size=1,prob=0.5) S = rep(ini,m) # we override most of S for(j in c(2:m)){ if(S[j-1]==0) S[j] <- rbinom(1,size=1,prob = (1-p1)) if(S[j-1]==1) S[j] <- (1-rbinom(1,size=1,prob = (1-q1))) } S = S[-c(1:1000)] return(S) } markov.switch.ar = function(a0,a1,p1,q1,n0,n){ m = (n+n0) y <- rep(0,m) S <- markov.chain1(n=m,p1=p1,q1=q1) innov = rnorm(m) for(i in c(2:m)){ if(S[i]==0){ y[i] = a0*y[i-1] + innov[i] } if(S[i]==1){ y[i] = a1*y[i-1] + innov[i] } } yS = cbind(y[-c(1:n0)],S[-c(1:n0)]) return(yS) } # Model 1 (both are in the stationary region) test = markov.switch.ar(a0=0.9,a1=-0.9,p1=0.9,q1=0.9,n0=500,n=300) # Model 2 (one is in the stationary region the other is not). But collectively # it is stationary test = markov.switch.ar(a0=1.5,a1=0.5,p1=0.2,q1=0.2,n0=1000,n=1000) # However if p1 is increased it will go outside the stationary region #pdf(file = '/my_plot.pdf') par(mfrow = c(2,1)) plot.ts(test[,1],col="blue") plot.ts(test[,2],col="blue") #dev.off()