Compute the equilibrium probabilities of the Markov chain using Jacobi Iteration in Python

Question

I am trying to compute a function that calculates the equilibrium probabilities of a Markov chain. For this problem, is have already got my transition matrix.

Now I am trying to define a function called Jacobi but confused about the most effective way to do so. Any suggestions on how to do this?

So far I have tried setting it us like a system of equations and solving x=a^(-1)*b but fail to implement it correctly due to the transition matrix being singular.

I know I need to multiple the transition matrix by a variable matrix to get 7 separate equations. Then I need to add the equation x0 + x1 + x2 + x3 + x4 + x5 + x6 = 1. After I have all 8 equations I can solve for x0 through x6 to get my equilibrium probabilities. Do you know how I can implement this process in python code?

Answer 1

I'm not sure the Jacobi method works for a Markov transition matrix. There are other, easier techniques for finding the stationary distribution though. First by solving the system of equations like you described:

import numpy as np


>>> M = np.array([[0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.5768, 0.224, 0.1494, 0.0498, 0.0, 0.0, 0.0], 
                  [0.3528, 0.224, 0.224, 0.1494, 0.0498, 0.0, 0.0], 
                  [0.1848, 0.168, 0.224, 0.224, 0.1494, 0.0498, 0.0], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498]])
>>>
>>> I = np.eye(len(M)) # identity matrix
>>>
>>> A = M.T - I  # rearrange each "equation" in the system
>>>
>>> A = np.vstack([A, np.ones((1, len(A)))]) # add row of ones
>>>
>>> b = np.zeros(len(A)) # prepare solution vector
>>> b[-1] = 1.0 # add a one to the last entry
>>>
>>> np.linalg.solve(A.T @ A, A.T @ b) # solve the normal equation
array([0.22288974, 0.14776044, 0.1781388 , 0.181211  , 0.1504296 ,
       0.09080838, 0.02876204])

You can also repeatedly multiply M by itself until it converges:

>>> np.linalg.matrix_power(M, 25)[0]
array([0.22288974, 0.14776044, 0.1781388 , 0.181211  , 0.1504296 ,
       0.09080838, 0.02876204])