A Simulation of Markov Chain

Post date: Feb 4, 2011 7:03:26 AM

mchainLen = 100000 # length of Markov chain

mchain = rep(NA, mchainLen) # initialize a vector of length mchainLen filled with "NA" (no values/missing values)

initvalue = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(0, 1/4, 3/4,0,0,0)) # initial value pmf = (0, 1/4, 3/4,0,0,0) on state space 1,2,3,4,5,6

initvalue

mchain[1] = initvalue

for (i in 2:100000)

{

## using transition probability below (t.p.m. specifies )

if (mchain[i-1] == 1) # conditional pmf if current state = 1

mchain[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(1/3, 0, 1/3,0,0,1/3)) # 1 row of P

if (mchain[i-1] == 2) # conditional pmf if current state = 2

mchain[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(1/2, 1/4, 1/4,0,0,0)) # 2 row of P

if (mchain[i-1] == 3) # conditional pmf if current state = 3

mchain[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(0, 0, 0,0,1,0)) # 3 row of P

if (mchain[i-1] == 4) # conditional pmf if current state = 4

mchain[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(1/4,1/4,1/4,0,0,1/4)) # 4 row of P

if (mchain[i-1] == 5) # conditional pmf if current state = 5

mchain[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(0,0,1,0,0,0)) # 5 row of P

if (mchain[i-1] == 6) # conditional pmf if current state = 6

mchain[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE, prob = c(0,0,0,0,0,1)) # 6 row of P

}

## calculate proportion of time spent in each of the states 1, 2, 3,4,5,6

sum(mchain==1)

sum(mchain==1)/mchainLen

sum(mchain==2)

sum(mchain==2)/mchainLen

sum(mchain==3)

sum(mchain==3)/mchainLen

sum(mchain==4)

sum(mchain==4)/mchainLen

sum(mchain==5)

sum(mchain==5)/mchainLen

sum(mchain==6)

sum(mchain==6)/mchainLen