# Stat 753 - Markov Chain Simulation (part I)
# February 13, 2020

# Simulate a Markov chain with transition probability P and initial condition x1.

# Function to simulate a Markov chain. 
# Requires the following input: n = number of time steps in sim, P = transition matrix, x1 = initial condition
MC.sim <- function(n,P,x1) {
  sim <- as.numeric(n)
  m <- ncol(P)
  if (missing(x1)) {
    sim[1] <- sample(1:m,1)  # random initial condition
  } 
  else { 
    sim[1] <- x1 
  }
  for (i in 2:n) {
    newstate <- sample(1:m,1,prob=P[sim[i-1],])
    sim[i] <- newstate
    }
  sim
}



#### Example 1 ####
###################

# Number of time steps to run the Markov chain
n = 100;
# n = 1000;

# Transition probability matrix P
P = matrix(c(0,0.5,0,1,0,1,0,0.5,0), nrow=3); P

# Double check that the rows sum to 1!
rowSums(P)

# Set an initial condition (or leave this input blank to generate a random initial condition)
nstates = ncol(P);
x1 = 1;  # start in state 1

### Generate a realization of the Markov chain - uses the initial condition specified above
MC.sim(n,P,x1)

# Plot the sample path
run1 <- MC.sim(n,P,x1)
X11(); 
plot(run1, type="b", pch=19, main="Sample Path of 3-State Markov Chain", xlab="Time", ylab="State")

# Make a histogram - shows the relative amount of time spent in each state
X11();
hist(run1, col="wheat", main="Histogram of Example 1", xlab="State")

### Generate a realization of the Markov chain - uses a random initial condition
run1b <- MC.sim(n,P)
plot(run1b, type="b", pch=19, main="Sample Path of 3-State Markov Chain", xlab="Time", ylab="State")
hist(run1b, col="wheat", main="Histogram of Example 1b", xlab="State")



#### Example 2 #### 
###################

# This example includes a self-loop. What do you notice in the sample path?

# Number of time steps to run the MC
n = 100;

# Transition probability matrix P
P = matrix(c(0,0,1,0.5,0.5,0,0.5,0.5,0), nrow=3); P

# Double check that the rows sum to 1!
rowSums(P)

# Set an initial condition (or leave this input blank to generate a random initial condition)
nstates = ncol(P);
x1 = 1;  # start in state 1

### Generate a realization of the Markov chain 
run2 <- MC.sim(n,P)
plot(run2, type="b", pch=19, main="Sample Path of 3-State Markov Chain", xlab="Time", ylab="State", col="blue")
hist(run2, col="light blue", xlab="State")




#### Example 3 #### 
###################

# This is a modified verson of Example 2. What do you observe in this case?

# Number of time steps to run the MC
n = 100;

# Transition probability matrix P
P = matrix(c(0,0,1,0.5,1,0,0.5,0,0), nrow=3); P

# Double check that the rows sum to 1!
rowSums(P)

# Set an initial condition (or leave this input blank to generate a random initial condition)
nstates = ncol(P);
x1 = 1;  # start in state 1

### Generate a realization of the Markov chain 
run3 <- MC.sim(n,P)
plot(run3, type="b", pch=19, main="Sample Path of 3-State Markov Chain", xlab="Time", ylab="State", col="blue")
hist(run3, col="light blue", xlab="State")

# State 2 in this example is called an absorbing state. Why is this the case?