Finding steady state probabilities by solving equation system

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)

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$.

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]$$