from __future__ import division
import numpy as np
import quantecon as qe
from mc_tools import mc_compute_stationary, mc_sample_path
from quantecon.mc_tools import MarkovChain
def ScaledGE(A):
A = np.array(A, dtype=float)
n = A.shape[0]
x = np.zeros(n)
#Reduction stage
# 1st: scaling column vector for i-th column
# 2nd: reduction.Obtain a submatrix for each i
for i in range(n-1):
A[i+1:n, i] /= -A[i,i]
A[i+1:n, i+1:n] += np.dot(A[i+1:n, i:i+1], A[i:i+1, i+1:n])
#Backward Substitution
x[n-1]=1
for i in range(n-2, -1, -1):
x[i] = np.sum((x[j] * A[j, i] for j in range(i+1, n)))
x /= np.sum(x)
return x
def GTH(A):
A = np.array(A, dtype=float)
n = A.shape[0]
x = np.zeros(n)
#Reduction stage
# 1st: scaling column vector for i-th column.Scaling is conducted by sum of other components in each row.
# 2nd: reduction.Obtain a submatrix for each i
# Key point is sum of row vector equals to 0. Keep an attention to this point.
for i in range(n-1):
scale = np.sum(A[i, i+1:n])
A[i+1:n, i] /= scale
A[i+1:n, i+1:n] += np.dot(A[i+1:n, i:i+1], A[i:i+1, i+1:n])
x[n-1] = 1
for i in range(n-2, -1, -1):
x[i] = np.sum((x[j] * A[j, i] for j in range(i+1, n)))
x /= np.sum(x)
return x
q=0.5
for i in range(8,18):
e=1e-1**i
P =np.array([[-(q+e), q , e ],
[q , -(q+e), e ],
[e , e , -2*e]])
print('e = 1e-{0}'.format(i))
print(ScaledGE(P))
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 ] e = 1e-17 [ nan nan -0.]
q=0.5
for i in range(8,18):
e=1e-1**i
P =np.array([[-(q+e), q , e ],
[q , -(q+e), e ],
[e , e , -2*e]])
print('e = 1e-{0}'.format(i))
print(GTH(P))
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-17 [ 0.33333333 0.33333333 0.33333333]
#Parameter
p = 1/3
N = 3
def KMR_2x2_P_sequential(N, p, epsilon):
P = np.zeros((N+1, N+1), dtype=float)
P[0, 0], P[0, 1] = 1 - epsilon * (1/2), epsilon * (1/2)
P[N, N-1], P[N, N] = epsilon * (1/2), 1 - epsilon * (1/2)
for x in range(1, N):
if x/N >p:
P[x, x-1] = (x/N)*epsilon * (1/2)
P[x, x+1] = (1-(x/N))*(1 - epsilon*(1/2))
P[x, x] = 1 - P[x, x-1] - P[x, x+1]
elif x/N<p:
P[x, x-1] =(x/N)*(1 - epsilon*(1/2))
P[x, x+1] =(1-(x/N))*epsilon * (1/2)
P[x, x] = 1 - P[x, x-1] - P[x, x+1]
else:
P[x, x-1] = (x/N)*(1/2)
P[x, x+1] = (1-(x/N))*(1/2)
P[x, x] = 1 - P[x, x-1] - P[x, x+1]
return P
for i in range(2,17):
epsilon=1e-1**i
kmr_seq=KMR_2x2_P_sequential(N, p, epsilon)
print('epsilon = 1-{0}'.format(i))
print(kmr_seq)
epsilon = 1-2 [[ 0.995 0.005 0. 0. ] [ 0.16666667 0.5 0.33333333 0. ] [ 0. 0.00333333 0.665 0.33166667] [ 0. 0. 0.005 0.995 ]] epsilon = 1-3 [[ 9.99500000e-01 5.00000000e-04 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-04 6.66500000e-01 3.33166667e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-04 9.99500000e-01]] epsilon = 1-4 [[ 9.99950000e-01 5.00000000e-05 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-05 6.66650000e-01 3.33316667e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-05 9.99950000e-01]] epsilon = 1-5 [[ 9.99995000e-01 5.00000000e-06 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-06 6.66665000e-01 3.33331667e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-06 9.99995000e-01]] epsilon = 1-6 [[ 9.99999500e-01 5.00000000e-07 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-07 6.66666500e-01 3.33333167e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-07 9.99999500e-01]] epsilon = 1-7 [[ 9.99999950e-01 5.00000000e-08 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-08 6.66666650e-01 3.33333317e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-08 9.99999950e-01]] epsilon = 1-8 [[ 9.99999995e-01 5.00000000e-09 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-09 6.66666665e-01 3.33333332e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-09 9.99999995e-01]] epsilon = 1-9 [[ 9.99999999e-01 5.00000000e-10 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-10 6.66666666e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-10 9.99999999e-01]] epsilon = 1-10 [[ 1.00000000e+00 5.00000000e-11 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-11 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-11 1.00000000e+00]] epsilon = 1-11 [[ 1.00000000e+00 5.00000000e-12 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-12 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-12 1.00000000e+00]] epsilon = 1-12 [[ 1.00000000e+00 5.00000000e-13 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-13 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-13 1.00000000e+00]] epsilon = 1-13 [[ 1.00000000e+00 5.00000000e-14 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-14 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-14 1.00000000e+00]] epsilon = 1-14 [[ 1.00000000e+00 5.00000000e-15 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-15 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-15 1.00000000e+00]] epsilon = 1-15 [[ 1.00000000e+00 5.00000000e-16 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-16 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-16 1.00000000e+00]] epsilon = 1-16 [[ 1.00000000e+00 5.00000000e-17 0.00000000e+00 0.00000000e+00] [ 1.66666667e-01 5.00000000e-01 3.33333333e-01 0.00000000e+00] [ 0.00000000e+00 3.33333333e-17 6.66666667e-01 3.33333333e-01] [ 0.00000000e+00 0.00000000e+00 5.00000000e-17 1.00000000e+00]]
for i in range(2,17):
epsilon=1e-1**i
kmr_seq=KMR_2x2_P_sequential(N, p, epsilon)
kmr_seq_m=qe.MarkovChain(kmr_seq)
kmr_seq_s=kmr_seq_m.stationary_distributions
print('epsilon = 1-{0}'.format(i))
print(kmr_seq_s)
epsilon = 1-2 [[ 4.92538049e-03 1.47761415e-04 1.47761415e-02 9.80150717e-01]] epsilon = 1-3 [[ 4.99250376e-04 1.49775113e-06 1.49775113e-03 9.98001501e-01]] epsilon = 1-4 [[ 4.99925004e-05 1.49977501e-08 1.49977501e-04 9.99800015e-01]] epsilon = 1-5 [[ 4.99992500e-06 1.49997750e-10 1.49997750e-05 9.99980000e-01]] epsilon = 1-6 [[ 4.99999250e-07 1.49999775e-12 1.49999775e-06 9.99998000e-01]] epsilon = 1-7 [[ 4.99999925e-08 1.49999978e-14 1.49999978e-07 9.99999800e-01]] epsilon = 1-8 [[ 4.99999993e-09 1.49999998e-16 1.49999998e-08 9.99999980e-01]] epsilon = 1-9 [[ 4.99999999e-10 1.50000000e-18 1.50000000e-09 9.99999998e-01]] epsilon = 1-10 [[ 5.00000000e-11 1.50000000e-20 1.50000000e-10 1.00000000e+00]] epsilon = 1-11 [[ 5.00000000e-12 1.50000000e-22 1.50000000e-11 1.00000000e+00]] epsilon = 1-12 [[ 5.00000000e-13 1.50000000e-24 1.50000000e-12 1.00000000e+00]] epsilon = 1-13 [[ 5.00000000e-14 1.50000000e-26 1.50000000e-13 1.00000000e+00]] epsilon = 1-14 [[ 5.00000000e-15 1.50000000e-28 1.50000000e-14 1.00000000e+00]] epsilon = 1-15 [[ 5.00000000e-16 1.50000000e-30 1.50000000e-15 1.00000000e+00]] epsilon = 1-16 [[ 5.00000000e-17 1.50000000e-32 1.50000000e-16 1.00000000e+00]]
for i in range(2,17):
epsilon=1e-1**i
kmr_seq=KMR_2x2_P_sequential(N, p, epsilon)
I=np.identity(kmr_seq.shape[0]) # used for SGE algorithm.
kmr_seq_SGE=ScaledGE(I-kmr_seq)
kmr_seq_GTH=GTH(I-kmr_seq)
print('epsilon = 1-{0}'.format(i))
print(kmr_seq_SGE)
print(kmr_seq_GTH)
epsilon = 1-2 [ 4.92538049e-03 1.47761415e-04 1.47761415e-02 9.80150717e-01] [ 4.92538049e-03 1.47761415e-04 1.47761415e-02 9.80150717e-01] epsilon = 1-3 [ 4.99250376e-04 1.49775113e-06 1.49775113e-03 9.98001501e-01] [ 4.99250376e-04 1.49775113e-06 1.49775113e-03 9.98001501e-01] epsilon = 1-4 [ 4.99925004e-05 1.49977501e-08 1.49977501e-04 9.99800015e-01] [ 4.99925004e-05 1.49977501e-08 1.49977501e-04 9.99800015e-01] epsilon = 1-5 [ 4.99992500e-06 1.49997750e-10 1.49997750e-05 9.99980000e-01] [ 4.99992500e-06 1.49997750e-10 1.49997750e-05 9.99980000e-01] epsilon = 1-6 [ 4.99999250e-07 1.49999775e-12 1.49999775e-06 9.99998000e-01] [ 4.99999250e-07 1.49999775e-12 1.49999775e-06 9.99998000e-01] epsilon = 1-7 [ 4.99999925e-08 1.49999977e-14 1.49999978e-07 9.99999800e-01] [ 4.99999925e-08 1.49999978e-14 1.49999978e-07 9.99999800e-01] epsilon = 1-8 [ 4.99999997e-09 1.49999998e-16 1.49999998e-08 9.99999980e-01] [ 4.99999993e-09 1.49999998e-16 1.49999998e-08 9.99999980e-01] epsilon = 1-9 [ 4.99999937e-10 1.49999994e-18 1.50000000e-09 9.99999998e-01] [ 4.99999999e-10 1.50000000e-18 1.50000000e-09 9.99999998e-01] epsilon = 1-10 [ 4.99999938e-11 1.49999994e-20 1.50000000e-10 1.00000000e+00] [ 5.00000000e-11 1.50000000e-20 1.50000000e-10 1.00000000e+00] epsilon = 1-11 [ 4.99999938e-12 1.49999994e-22 1.50000000e-11 1.00000000e+00] [ 5.00000000e-12 1.50000000e-22 1.50000000e-11 1.00000000e+00] epsilon = 1-12 [ 4.99933333e-13 1.49993333e-24 1.50000000e-12 1.00000000e+00] [ 5.00000000e-13 1.50000000e-24 1.50000000e-12 1.00000000e+00] epsilon = 1-13 [ 5.00600178e-14 1.50060018e-26 1.50000000e-13 1.00000000e+00] [ 5.00000000e-14 1.50000000e-26 1.50000000e-13 1.00000000e+00] epsilon = 1-14 [ 5.00600178e-15 1.50060018e-28 1.50000000e-14 1.00000000e+00] [ 5.00000000e-15 1.50000000e-28 1.50000000e-14 1.00000000e+00] epsilon = 1-15 [ 4.29061342e-16 1.42906134e-30 1.50000000e-15 1.00000000e+00] [ 5.00000000e-16 1.50000000e-30 1.50000000e-15 1.00000000e+00] epsilon = 1-16 [ nan nan nan nan] [ 5.00000000e-17 1.50000000e-32 1.50000000e-16 1.00000000e+00]