Calculate Expected States of Markov Chain via Simulation

Post date: Feb 4, 2011 7:06:23 AM

ssize = 1000 # sample size of simulation

x0 = rep(NA, ssize) # initialize a vector of x0 filled with "NA" (no values/missing values)

x1 = rep(NA, ssize) # initialize a vector of x1 filled with "NA" (no values/missing values)

x2 = rep(NA, ssize) # initialize a vector of x2 filled with "NA" (no values/missing values)

x3 = rep(NA, ssize) # initialize a vector of x3 filled with "NA" (no values/missing values)

for (i in 1:1000)

{

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

## realization of x0

x0[i]= 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

## realization of x1

if (x0[i] == 2) # conditional pmf if current state = 2

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

if (x0[i] == 3) # conditional pmf if current state = 3

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

## realization of x2

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

x2[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 (x1[i] == 2) # conditional pmf if current state = 2

x2[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 (x1[i] == 3) # conditional pmf if current state = 3

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

if (x1[i] == 5) # conditional pmf if current state = 5

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

## realization of x3

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

x3[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 (x2[i] == 2) # conditional pmf if current state = 2

x3[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 (x2[i] == 3) # conditional pmf if current state = 3

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

if (x2[i] == 5) # conditional pmf if current state = 5

x3[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 (x2[i] == 6) # conditional pmf if current state = 6

x3[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 expected state of x0,x1,x2,x3:

sum(x0)/ssize

sum(x1)/ssize

sum(x2)/ssize

sum(x3)/ssize