############################################################################################
# Simple Markov chain example
############################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
############################################################################################# 
# Transition matrix
P=matrix(c(
0.3333,0.6667,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
0.0000,0.3283,0.5672,0.1045,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
0.0000,0.0372,0.4288,0.4921,0.0405,0.0014,0.0000,0.0000,0.0000,0.0000,
0.0000,0.0005,0.0796,0.5366,0.3653,0.0179,0.0001,0.0000,0.0000,0.0000,
0.0000,0.0001,0.0028,0.1481,0.5964,0.2445,0.0081,0.0000,0.0000,0.0000,
0.0000,0.0000,0.0000,0.0075,0.2463,0.6029,0.1409,0.0024,0.0000,0.0000,
0.0000,0.0000,0.0000,0.0000,0.0189,0.3567,0.5468,0.0770,0.0006,0.0000,
0.0000,0.0000,0.0000,0.0000,0.0005,0.0329,0.4891,0.4460,0.0315,0.0000,
0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0625,0.5938,0.3359,0.0078,
0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.6000,0.4000),
10,10,byrow=T)

# Computing the actual limiting distribution
P1 = P
for (i in 1:10000)
  P1 = P1%*%P

# Simulating from the Markov chain
x = rep(0,200000)
x[1]=1
for (i in 2:200000)
  x[i]=sample(1:10,size=1,P[x[i-1],],replace=T)

# Computing the observed limiting distribution
ps = rep(0,10)
for (i in 1:10)
  ps[i] = mean(x[100001:200000]==i)

# Comparing the actual limiting distribution with
# the observed one
cbind(c(1:10),round(P1[1,],4),ps)

# Graphical summary
pdf(file="transition.pdf",height=7,width=10)
par(mfrow=c(1,1))
PP = matrix(P,100,1)
xx = rep(1:10,10)
yy = rep(1,10)
for (i in 2:10)
  yy = c(yy,rep(i,10))
plot(yy,xx,col=0,xlab="state tomorrow",ylab="state today",axes=F)
axis(1,at=1:10)
axis(2,at=1:10)
text(yy,xx,round(PP,3),col=2)
dev.off()

pdf(file="equilibrium.pdf",height=7,width=10)
plot(1:10,round(P1[1,],4),type="h",xlab="state",ylab="",axes=F)
title("Equilibrium distribution")
axis(1,1:10)
axis(2)
lines(1:10,round(P1[1,],4),lwd=2,type="h")
lines(1:10+0.1,ps,lwd=2,lty=2,col=2,type="h")
legend(1,0.3,legend=c("True","Estimated"),lty=1:2,col=1:2,lwd=c(2,2),bty="n")
dev.off()