# simulate discrete Markov chains according to transition matrix P
run.mc.sim <- function(P, init.state, num.iters) {
  
  # number of possible states
  num.states <- nrow(P)
  
  # stores the states X_t through time
  states     <- numeric(num.iters)
  
  # initialize variable for first state 
  states[1]    <- init.state
  
  for(t in 2:num.iters) {
    
    # probability vector to simulate next state X_{t+1}
    p  <- P[states[t-1], ]
    
    ## draw from multinomial and determine state
    states[t] <-  which(rmultinom(1, 1, p) == 1)
  }
  return(states)
}

P.ex <- rbind(c(0.7,0.25,0.05),c(0.6,0.3,0.1),c(0,0,1))
output1 <- run.mc.sim(P.ex, 1, 50)
output2 <- run.mc.sim(P.ex, 1, 50)
plot(31:80,output1, xlab="State", ylab="Time", col="red", type="p")
lines(31:80,output2, xlab="State", ylab="Time", col="blue", type="p", pch=0)