## Functions to check properties of Markov chains by simulation
##
## Author: Andrea Riebler, 15.09.2014
##
## Comment: Compared to the lecture demonstration I added more 
## comments and splitted the function we wrote in two functions 
## one to simulate a Marko chain, called simMC, and one to get the
## relative frequencies of all states in the state space as time 
## goes on, called genFreqTable.

## 1.) Function to simulate a Markov chain
# What do we need:
# q : vector of initial probabilities
# P : matrix with transition probabilities
# n : length of the Markov chain
# states : labels of the states

simMC <- function(q, P, n, states){
  
  # generate a vector of the desired length
  myMC <- rep(NA, n)
  # sample the initial state from the states vector using
  # the initial probabilities.
  myMC[1] <- sample(states, 1, prob=q)
  # conditioning on where we at time i-1 we sample the next state 
  # at time i.
  for(i in 2:n){
    # we know the state at time i-1, which is myMC[i-1], to get the
    # corresponding line in P we use the function match. 
    # The function math returns the position of myMC[i-1] in the 
    # states vector, i.e the line in P that belongs to state myMC[i-1].
    myMC[i] <- sample(states, 1, prob=P[match(myMC[i-1], states),])
  }
  # return the Markov chain
  return(myMC)
}

## Try the function
# desired length of the Markov chain
n <- 100
# state space for the inventory models
states <- c(-1, 0, 1, 2)
# vector of initial probabilities.  Here, we say
# to start in state 2.
q <- c(0, 0, 0, 1)
# transition matrix from the inventory model
P <- matrix(c(0, 0.1, 0.4, 0.5,
             0, 0.1, 0.4, 0.5,
             0.1, 0.4, 0.5, 0,
             0, 0.1, 0.4, 0.5), 4,4, byrow=T)

MC <- simMC(q=q, P=P, n=n, states=states)
plot(MC, type="o", ylab="States", xlab="time")


## 2. Function to return the cumulative frequency distribution
## of the states in a Markov chain
# Arguments:
# myMC : vector representing the Markov chain 
#       (returned for example by simMC)
# states : vector of states that are potentially possible
genFreqTable <- function(myMC, states){

  # length of the Markov chain
  n <- length(myMC)
  myMC = factor(myMC, levels=states)
  # generate a matrix to save the frequencies of all states
  # as the MC increases.
  freqStates <- matrix(NA, ncol=length(states), nrow=n)
  colnames(freqStates) = levels(myMC)
  for(i in 1:n){
    freqStates[i,] <- table(myMC[1:i])
  }
  # transform the absolute frequencies to relative frequencies
  return(freqStates/1:n)
}
MC <- simMC(q=q, P=P, n=10, states=c(2,4,1, 9))
genFreqTable(MC, c(2,4,1, 9))

## Try the function
# Generate a longer chain to check the long-time behaviour
MC <- simMC(q=q, P=P, n=10000, states=states)
# get a table on the relative frequencies of the 
# state distribution after each time step.
fT <- genFreqTable(MC, states)
# use matplot to plot the relative frequencies
matplot(fT, type="l", lwd=2, xlab="time", ylab="Relative frequencies")
# add a grid in the background 
grid()
# add a legend
legend("topright", legend=states, col=1:4, lty=1)