# Markov Chain code

pp<-matrix(c(.1,.2,.7,.2,.2,.6,.6,.1,.3),ncol=3,byrow=T)

z<-numeric(100000)
z[1]<-1
for(i in 2:100000){
	z[i] <- sample(c(1,2,3),prob=pp[z[i-1],],size=1)
}

z.converged <- z[11:100000]
length(z.converged[z.converged==1])/length(z.converged)
# or mean(z.converged==1)
length(z.converged[z.converged==2])/length(z.converged)
# or mean(z.converged==2)
length(z.converged[z.converged==3])/length(z.converged)
# or mean(z.converged==3)

hist(z.converged)

tpp <- matrix(c(-.9,.2,.6,.2,-.8,.1,1,1,1),ncol=3,byrow=T)
solve(tpp,c(0,0,1))

pp2 <- pp%*%pp
pp4 <- pp2%*%pp2
pp8 <- pp4%*%pp4
pp16 <- pp8%*%pp8
pp32 <- pp16%*%pp16
pp64 <- pp32%*%pp32

pp64