Finding steady state probabilities by solving equation system

Question

(I know that there are numerous questions on this, but my problem is in actually solving the equations, which isn't the problem in other questions.)

I'm trying to figure out the steady state probabilities for a Markov Chain, but I'm having problems with actually solving the equations that arise. So,

$$ V = V\times \left[ { \begin{array}{ccc} 0 & 0.4 & 0 & 0.6\\ 0 & 0.5 & 0.5 & 0 \\ 0.4 & 0 & 0.2 & 0.4\\ 0 & 0 & 1 & 0 \end{array} }\right]$$ Where $V = \left[ \begin{array}{cccc}a & b & c & d \end{array} \right]$

So typing out the equations that I get from multiplying V with the transition matrix.

\begin{align} &0a &+ &0b &+ &0.4c &+ &0d&= a\\ &0.4a &+ &0.5b &+ &0c &+ &0d &= b\\ &0a &+ &0.5b &+ &0.2c &+ &d &=c\\ &0.6a &+ &0b &+ &0.4c &+ &0d &=d \end{align}

How do I solve these equations?

I know the answer for a, b, c, d but I get lost in trying to get to the answer from the equations.

Answer 1

Start with the good equations:

\begin{align} &0a &+ &0b &+ &0.4c &+ &0d&= a\\ &0.4a &+ &0.5b &+ &0c &+ &0d &= b\\ &0a &+ &0.5b &+ &0.2c &+ &d &=c\\ &0.6a &+ &0b &+ &0.4c &+ &0d &=d \end{align}

---

now you can just substitute:

\begin{align} &-a & & & &0.4c & &&= 0\\ &0.4a &- &0.5b & & & &&= 0\\ & & &0.5b &- &0.8c &+ &d &=0\\ &0.6a & & &+ &0.4c &- &d &=0 \end{align}

substitute $a = 0.4c$: \begin{align} -&0.5b &+ &0.16c & &&= 0\\ &0.5b &- &0.8c &+ &d &=0\\ && &0.64c &- &d &=0 \end{align}

substitute $d = 0.64c$: \begin{align} -&0.5b&+&0.16c =0\\ &0.5b&-&0.16c =0 \end{align} you get a consistent system.

This confirms that the system has solutions. The general solution is

$$(0.4c, 0.32c, c, 0.64c)$$

Now choose $c$ such as $a+b+c+d = 1$ and check that $a>0, b>0, c>.0, d>0$.

Answer 2

The equation you want to solve is $[a\ b\ c\ d](V-I)=0$, where $I$ is the identity matrix which has $1$'s on the diagonal and $0$'s elsewhere. To solve the equation, use Gaussian elimination.

$$V-I=\left[ { \begin{array}{ccc} -1 & 0.4 & 0 & 0.6\\ 0 & -0.5 & 0.5 & 0 \\ 0.4 & 0 & -0.8 & 0.4\\ 0 & 0 & 1 & -1 \end{array} }\right]$$

Answer 3

Another thing you can do, which makes the procedure completely general - in addition to the four equations you get from

$$v(P-I) = 0$$

Where $P$ is your transition probabilities matrix above, you also need -

$$v_1 + v_2 +v_3 + \dots = 1$$ So, just add a column to your $P-I$ matrix with all ones. This new matrix has n rows and n+1 columns. Let's call it $Q$. And your right hand side vector becomes $b = [0, 0, \dots, 1]$, meaning repeat zeros until n and then add a 1. Now you get -

$$vQ = b$$ Post multiply both sides by $Q^T$

$$vQQ^T = bQ^T$$

Which leads to -

$$v = bQ^T(QQ^T)^{-1}$$

Here, $v$ is a column vector. Most linear algebra solvers expect the inverse matrix first multiplied by the vector, which is not the form of the result above. To convert it into this format (so it can be easily passed to these solvers), simply take the transpose.

$$v^T = (Q^T Q)^{-1}(Qb^T)$$

However, note that $Qb^T$ is simply a vector of ones (since all but the last column of $b$ is zero and the last row of $Q^T$ is all ones). So, no need to explicitly calculate $Qb^T$ it is simply a vector of ones. Here is a python method that does this for you:

import numpy as np

def steady_state_prop(
  p=np.matrix([
          [0,.2,.4,.4],
          [.2,0,.4,.4],
          [.2,.3,0,.5],
          [.3,.4,.3,0]
          ])):
    dim = p.shape[0]
    q = (p-np.eye(dim))
    ones = np.ones(dim)
    q = np.c_[q,ones]
    QTQ = np.dot(q, q.T)
    bQT = np.ones(dim)
    return np.linalg.solve(QTQ,bQT)