

In [1]:

    
from __future__ import division, print_function
import numpy as np
import quantecon as qe

Computing Stationary Distributions

(Scaled) Gaussian Elimination



In [2]:

    
def scaled_ge_solve(A):
    """
    Solve x A = 0 for a nonzero (normalized) solution x,
    where A is an irreducible matrix with nonnegative off-diagonals
    and zero row sums, by (scalled) Gaussian elimination.

    """
    A = np.array(A, dtype=float)  # Make a copy of A
    n = A.shape[0]
    x = np.zeros(n)

    # === Reduction === #
    for k in range(n-1):
        for i in range(k+1, n):
            A[i, k] /= -A[k, k]

            for j in range(k+1, n):
                A[i, j] += A[i, k] * A[k, j]


    # === Backward substitution === #
    x[n-1] = 1
    for k in range(n-2, -1, -1):
        for i in range(k+1, n):
            x[k] += x[i] * A[i, k]

    # === Normalization === #
    x /= np.sum(x)

    return x

Let us apply the above to Hamilton's Markov chain:



In [3]:

    
P = [[0.971, 0.029, 0    ],
     [0.145, 0.778, 0.077],
     [0    , 0.508, 0.492]]



In [4]:

    
scaled_ge_solve(P - np.identity(3))









    Out[4]:





array([ 0.8128 ,  0.16256,  0.02464])

Check that the output is correct:



In [5]:

    
qe.gth_solve(P)









    Out[5]:





array([ 0.8128 ,  0.16256,  0.02464])

Consider the matrix \[ A = \begin{pmatrix}

(q + \varepsilon) & q & \varepsilon \ q & - (q + \varepsilon) & \varepsilon \ \varepsilon & \varepsilon & - 2\varepsilon \end{pmatrix}, \] where $0 < q < 1$, and $\varepsilon > 0$ is sufficiently small.

Observe that this has a (unique) stationary disribution $(\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$ for any values of $q$ and $\varepsilon$.

Apply scaled_ge_solve to $A$ with $q = 0.5$ and $\varepsilon = 10^{-8}, \ldots, 10^{-16}$:



In [6]:

    
q = 0.5

for i in range(8, 17):
    e = 1e-1 ** i
    A = [[-(q+e), q     , e   ],
         [q     , -(q+e), e   ],
         [e     , e     , -2*e]]
    print('e = 1e-{0}'.format(i))
    print(scaled_ge_solve(A))









    



e = 1e-8
[ 0.33333333  0.33333333  0.33333333]
e = 1e-9
[ 0.33333334  0.33333334  0.33333333]
e = 1e-10
[ 0.33333332  0.33333332  0.33333335]
e = 1e-11
[ 0.33333332  0.33333332  0.33333335]
e = 1e-12
[ 0.33333579  0.33333579  0.33332842]
e = 1e-13
[ 0.33329879  0.33329879  0.33340243]
e = 1e-14
[ 0.33342217  0.33342217  0.33315567]
e = 1e-15
[ 0.33342217  0.33342217  0.33315567]
e = 1e-16
[ 0.32152035  0.32152035  0.3569593 ]

You can see that the output becomes far from $(\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$ as $\varepsilon$ becomes small.

GTH Algorithm



In [7]:

    
def gth_solve(A):
    """
    Solve x A = 0 for a nonzero (normalized) solution x,
    where A is an irreducible matrix with nonnegative off-diagonals
    and zero row sums, by the GTH algorithm.

    """
    A = np.array(A, dtype=float)  # Make a copy of A
    n = A.shape[0]
    x = np.zeros(n)

    # === Reduction === #
    for k in range(n-1):
        scale = np.sum(A[k, k+1:n])
        
        for i in range(k+1, n):
            A[i, k] /= scale

            for j in range(k+1, n):
                A[i, j] += A[i, k] * A[k, j]


    # === Backward substitution === #
    x[n-1] = 1
    for k in range(n-2, -1, -1):
        for i in range(k+1, n):
            x[k] += x[i] * A[i, k]

    # === Normalization === #
    x /= np.sum(x)

    return x

Verify that gth_solve always returns the correct answer:



In [8]:

    
for i in list(range(8, 17)) + [100]:
    e = 1e-1 ** i
    A = [[-(q+e), q     , e   ],
         [q     , -(q+e), e   ],
         [e     , e     , -2*e]]
    print('e = 1e-{0}'.format(i))
    print(gth_solve(A))









    



e = 1e-8
[ 0.33333333  0.33333333  0.33333333]
e = 1e-9
[ 0.33333333  0.33333333  0.33333333]
e = 1e-10
[ 0.33333333  0.33333333  0.33333333]
e = 1e-11
[ 0.33333333  0.33333333  0.33333333]
e = 1e-12
[ 0.33333333  0.33333333  0.33333333]
e = 1e-13
[ 0.33333333  0.33333333  0.33333333]
e = 1e-14
[ 0.33333333  0.33333333  0.33333333]
e = 1e-15
[ 0.33333333  0.33333333  0.33333333]
e = 1e-16
[ 0.33333333  0.33333333  0.33333333]
e = 1e-100
[ 0.33333333  0.33333333  0.33333333]

For a vectorized version of gth_solve, see gth_solve.py in quantecon.



In [8]: