Finding stationary distribution of a markov process given a transition probability matrix

Question

There has been two threads related to this issue on Stack Overflow:

• How can I obtain stationary distribution of a Markov Chain given a transition probability matrix describes what a transition probability matrix is, and demonstrate how a stationary distribution is reached by taking powers of this matrix;
• How to find when a matrix converges with a loop uses an R loop to determine when the matrix power converges.

The above is straightforward, but very expensive. If we have a transition matrix of order n, then at each iteration we compute a matrix-matrix multiplication at costs O(n ^ 3).

Is there a more efficient way to do this? One thing that occurs to me is to use Eigen decomposition. A Markov matrix is known to:

• be diagonalizable in complex domain: A = E * D * E^{-1};
• have a real Eigen value of 1, and other (complex) Eigen values with length smaller than 1.

The stationary distribution is the Eigen vector associated with the Eigen value of 1, i.e., the first Eigen vector.

Well, the theory is nice, but I can't get it work. Taking the matrix P in the first linked question:

P <- structure(c(0, 0.1, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0.2, 
0.3, 0, 0, 0.5, 0.4, 0.3, 0.5, 0.4, 0, 0, 0, 0, 0, 0.6, 0.4, 
0.5, 0.4, 0.3, 0.2, 0, 0.6), .Dim = c(6L, 6L))

If I do:

Re(eigen(P)$vectors[, 1])
# [1] 0.4082483 0.4082483 0.4082483 0.4082483 0.4082483 0.4082483

What's going on? According to previous questions, the stationary distribution is:

# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

Answer 1

Your vector y = Re(eigen(P)$vectors[, 1]) is not a distribution (since it doesn't add up to one) and solves P'y = y, not x'P = x. The solution from your linked Q&A does approximately solve the latter:

x = c(0.00259067357512953, 0.0259067357512953, 0.116580310880829, 
0.310880829015544, 0.272020725388601, 0.272020725388601)
all(abs(x %*% P - x) < 1e-10) # TRUE

By transposing P, you can use your eigenvalue approach:

x2 = Re(eigen(t(P))$vectors[, 1])
x2 <- x2/sum(x2) 
(function(x) all(abs(x %*% P - x) < 1e-10))(
  x2
) # TRUE

It's finding a different stationary vector in this instance, though.