Plot how the probability P(Xn=s) changes as a function of time?

Question

Entire problem: We have a Markov chain model with with 5 states: s, t, m, f, r

TPM follows:

      P <- matrix(c(.84,.03,.01,.03,.03,
                    .11,.80,.15,.19,.09,
                    .01,.04,.70,.02,.05,
                    .04,.10,.07,.75,.00,
                    .00,.03,.07,.01,.83),
                  nrow=5

                   )

With matrix multiplication, the limiting distribution comes out to:

(0.1478365,0.4149259,0.09555939,0.2163813,0.1252968)

I am attempting to plot how P(Xn = s) changes as a function of time.

Given the initial distribution is P(X0 = i) = 1/5 i.e:

                  s   t   m   f   r  
           α = ( 1/5 1/5 1/5 1/5 1/5 )

I need to plot P(Xn = s) (on the y-axis) against n = 0, 1, 2, 3, 4, 5 (x-axis).

Answer 1

The expm package has a matrix-power function, %^%. I first checked the 10th state, but it has not reach stability so I went with the 20th state

Each successive state is in the progression a %*% P, a %*% (P%^%2), ...., a %*% (P%^%20)

library(expm)

#Attaching package: ‘expm’

#The following object is masked from ‘package:Matrix’:

 #   expm

So I guess I had it already available, no matter.

evolve <- sapply(1:20, function(n){ a%*% (P %^% n)})  # 5 x 20 matrix
png(); matplot( 1:20, t(evolve) );dev.off()  # matplot needs data in rows

Answer 2

Here is a base R solution which enables the matrix-power operation and view the progress of evolution

n <- 20
res <- do.call(rbind,Reduce(`%*%`,c(list(a),replicate(n,P,simplify = FALSE)),accumulate = T))

then, similar to the plotting approach by @42-, you can use to see the evolution w.r.t. 1:n

matplot(seq(nrow(r)),r)