##### BEGIN USER INPUT ##### nIter <- 10000 # Number of simulated steps of the Markov chain P = matrix(c(0.8,0.1,0.1,0.2,0.6,0.2,0.3,0.3,0.4),3,3, byrow=TRUE) # Define a transition matrix ##### END USER INPUT ##### #install.packages("expm", repos="http://R-Forge.R-project.org") # To get the matrix power library(expm) lineColors = c("red","blue","green","black","yellow") nStates <- dim(P)[1] print(P) # Stationary distribution is normalized left eigenvector of P corresponding to eigenvalue of one. eigObj <- eigen(t(P)) eigenVect <- eigObj\$vector[,eigObj\$values == 1] piStat <- eigenVect/sum(eigenVect) print(piStat) print(piStat - piStat%*%P) # Check 1: should be the zeros vector if piStat is the stationary distribution print(P%^%100) # Check 2: each row should be approximately piStat # Simulating a Markov chain x <- matrix(NA,nIter,1) x[1] <- 3 # Initial value for (i in 2:nIter){ x[i] <- sample(1:nStates,1, prob = P[x[i-1],]) } # Graphics print(x[1:50]) # Printing the first 50 simulations approxStat <- list(nStates,1) for (j in 1:nStates){ approxStat[[j]] <- cumsum(x == j)/seq(1,nIter) } plot(0, type = "n", xlim = c(0,nIter), ylim = c(0,1), ylab = "Marginal probability", xlab = "Iteration number", lwd = 2) legendText <- c() for (j in 1:nStates){ lines(1:nIter,approxStat[[j]], col = lineColors[j], lwd = 2) legendText <- c(legendText, paste('Cumulative mean - state',j), paste('True stationary - state',j)) abline(piStat[j], 0, col = lineColors[j], lty = 2) } legend("topright", inset = 0.05, box.lty = 1, legend = legendText, col = rep(lineColors[1:nStates], each=2), lty = rep(c(1,2),nStates), lwd = rep(c(2,1),nStates), cex = 0.5)