# 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()