# Stat 753 - Markov Chain Simulation (part III)
# See https://cran.r-project.org/web/packages/markovchain/markovchain.pdf for details
# February 20, 2020

# Install the markovchain package in R
library(markovchain)


# First define the transition matrix
tmA <- matrix(c(0,0.5,0.5,0.5,0,0.5,0.5,0.5,0), nrow=3, byrow=TRUE);
tmA

# Create a discrete time Markov chain (DTMC)
dtmcA <- new("markovchain", transitionMatrix=tmA, states=c("a","b","c"), name="Markov Chain A");
dtmcA

# Shortcut way to create the same DTMC:
#dtmcA2 <- as(tmA, "markovchain");
#dtmcA2
#states(dtmcA2)  # default names assigned to states

# Display the graph (uses the igraph package)
plot(dtmcA)

# If desired, one can access the transition probabilities, for example
dtmcA[2,3]

# Or by calling...
transitionProbability(dtmcA,"b","c")


# Simulate the distribution after n-steps
initialState <- c(0,1,0)
steps <- 5
finalState <- initialState * dtmcA^steps
finalState

# Show a realization of the MC 
markovchainSequence(n=20, markovchain=dtmcA, include=TRUE)

# Compute the steady state distribution
steadyStates(dtmcA)

# Describe absorbing states (if any)
absorbingStates(dtmcA)

# Summary info - classification of states
summary(dtmcA)



### Another DTMC example
tmB <- matrix(c(0.2,0.5,0.3,0,1,0,0.1,0.8,0.1), nrow=3, byrow=TRUE);

dtmcB <-new("markovchain", transitionMatrix=tmB, states=c("a","b","c"), name="Markov Chain B");
dtmcB

# Display the graph (uses the igraph package)
plot(dtmcB)

# Simulate the distribution after n-steps
initialState <- c(1,0,0)
steps <- 5
finalState <- initialState * dtmcB^steps
finalState

# Show a realization of the MC 
markovchainSequence(n=20, markovchain=dtmcB, include=TRUE)

# Compute the steady state distribution
steadyStates(dtmcB)

# Describe absorbing states (if any)
absorbingStates(dtmcB)

# Summary info - classification of states
summary(dtmcB)



### In-class example 
tmC <- matrix(c(0,0,0.5,0.5,1,0,0,0,0,1,0,0,0,1,0,0), nrow=4, byrow=TRUE);

dtmcC <-new("markovchain", transitionMatrix=tmC, states=c("0","1","2","3"), name="Markov Chain Class Example");
dtmcC

plot(dtmcC)

# Compute the steady state distribution
steadyStates(dtmcC)

# Summary info - classification of states
summary(dtmcC)



### In-class example 2
tmD <- matrix(c(0.5,0.5,0,0,0.5,0.5,0,0,0.25,0.25,0.25,0.25,0,0,0,1), nrow=4, byrow=TRUE);

dtmcD <-new("markovchain", transitionMatrix=tmD, states=c("0","1","2","3"), name="Markov Chain Class Example 2");
dtmcD

plot(dtmcD)

# Compute the steady state distribution
steadyStates(dtmcD)

# Summary info - classification of states
summary(dtmcD)


# Extra things you can do
dtmcD^10

# Compare to raising the transition matrix to the 10th power
library(expm)
tmD %^% 10