Contents

This computational worksheet is an opportunity to learn more about Markov chains using R . This practical doubles as Assessment  2, and is an assessed part of the course.  This assessment counts for 3% of your final module grade.

You should write a report on your work (see Question 0 below), answering all the questions with explanations and figures where necessary, and showing any R code you used.  Good formats for your report include PDF, Word document, or HTML . You work must be a report: do not just submit a list of R commands.

The deadline for submission is Monday  27 March  at  2pm. Late submissions up to Thursday 3 April at 2pm will still be marked, but the total mark will be reduced by 10% for each day or part-day it is late.

Your report will by marked out of 8, with 2 marks available for each of Questions 1 –4 .

Question  0:  Answer the questions on this problem sheet.   Write up a short report on your findings, with answers to the questions. Make sure to include relevant graphs produced and the R code you used.

Remember that this is assessed work, and your report must be submitted by Monday 27 March at 2pm. You might like to look at the example report for Computational Worksheet 1 to see what a good report looks like.

We start by reading three functions into R by entering the entire text of the following code chunk.  If you are using the Rmd version of this question sheet in RStudio, you can simply the green “play”button next to the following code chunk. If you are reading the HTML version, you should carefully copy and paste the following code chunk into R.

# Produces a transition matrix for a Markov chain based on a parameter eps

TransMat <- function(eps) {

matrix(c(0.1, 0.9,  0,      0,      0,  0,

0.4, 0.2, 0.4,      0,      0,  0,

0, 0.5,  0,    0.5,      0,  0,

0,  0, 0.5, 0.5-eps,    eps,  0,

0,  0,  0,    eps, 0.3-eps, 0.7,

0,  0,  0,      0,    0.5, 0.5),

nrow = 6, ncol = 6, byrow = TRUE)

}

# Simulates a Markov chain with transition matrix P for N steps from X0

MarkovChain <- function(N, X0, P) {

X <- rep(0, N)    # empty vector of length N

space <- nrow(P)  # number of points in sample space

now <- X0

for (n in 1:N) {

now <- sample(space, 1, prob = P[now, ])

X[n] <- now

}

X

}

# Finds the stationary distribution of a Markov chain

StatDist <- function(P) {

LeftEigen <- eigen(t(P))

statmeas <- LeftEigen$vectors[, 1]  # Unnormalised stationary measure

statmeas / sum(statmeas)            # Normalised stationary distribution

}

It is not important that you understand how these functions work, but you do need to know what they do. The first function TransMat() produces a transition matrix

  ⎞⎟

⎜                                                                                                                                                                                         ⎟

⎜  0        0.5        0            0.5                 0              0     ⎟

⎜                                                                                                                                                                                         ⎟

⎜                                                                                                                                                    ⎟

⎝  0        0        0            0              0.5        0.5 ⎠

for some value of the parameter  .

The second function, MarkovChain(), simulates a Markov chain.  For example, MarkovChain(10, 2, P) simulates a Markov chain with transition matrix P, starting from state 2 for 10 steps.

The third function, StatDist(), calculates a stationary distribution for a given transition matrix.  This function works for any irreducible Markov chain on a finite state space – it may behave strangely for non-

irreducible Markov chains.

We will be investigating the Markov chain with  = 0 .2 .

Question 1.  Using the function TransMat(), produce the transition matrix for  = 0 .2 .  Using the function MarkovChain(), simulate N  = 50  steps of this Markov chain starting from x0   =  1 .  Check it has worked, perhaps by plotting a graph or examining the vector produced.  Comment on what you find.

You will presumably want to try something like this:

eps <- 0.2

P   <- TransMat(eps)

N   <- 50