Reading: J. P. Jarvis and D. R. Shier, "Graph-Theoretic Analysis of Finite Markov Chain."
In [1]:
from __future__ import division
import numpy as np
import quantecon as qe
Make sure that you have installed quantecon version 0.1.6 (or above):
In [2]:
qe.__version__
Out[2]:
Consider the Markov chain given by the following stochastic matrix (where the actual values of non-zero probabilities are not important):
In [3]:
P = np.zeros((6, 6))
P[0, 0] = 1
P[1, 4] = 1
P[2, [2, 3, 4]] = 1/3
P[3, [0, 5]] = 1/2
P[4, [1, 4]] = 1/2
P[5, [0, 3]] = 1/2
In [4]:
print P
Create a MarkovChain instance:
In [5]:
mc0 = qe.MarkovChain(P)
We call the states $0, \ldots, 5$, respectively, instead of $1, \ldots, 6$.
(a) Determine the communication classes.
In [6]:
mc0.communication_classes
Out[6]:
(b) Classify the states of this Markov chain.
In [7]:
mc0.recurrent_classes
Out[7]:
Obtain a list of the recurrent states.
In [8]:
# Write your own code
recurrent_states = np.concatenate(mc0.recurrent_classes)
print recurrent_states
Obtain a list of the transient states.
In [9]:
# Write your own code
transient_states = np.setdiff1d(np.arange(mc0.n), recurrent_states)
print transient_states
(c) Does the chain have a unique stationary distribution?
In [10]:
mc0.num_recurrent_classes
Out[10]:
Obtain the stationary distributions:
In [11]:
vecs = mc0.stationary_distributions
print vecs
Verify that the above vectors are indeed stationary distributions.
Hint: Use np.dot.
In [12]:
# Write your own code
print np.dot(vecs, mc0.P)
Let us simulate our Markov chain mc0.
The simualte method generates a sample path from an initial state as specified by the init argument,
of length specified by sample_size, which is set to 1000 by default when omitted.
A sample path from state 0:
In [13]:
mc0.simulate(init=0, sample_size=50)
Out[13]:
As is clear from the transition matrix P,
if it starts at state 0, the chain stays there forever,
i.e., 0 is an absorbing state, a state that constitutes a singleton recurrent class.
Start with state 1:
In [14]:
mc0.simulate(init=1, sample_size=50)
Out[14]:
You can observe that the chain stays in the recurrent class $\{1, 4\}$ and
visits states 1 and 4 with certain frequencies.
Let us compute the frequency distribution along a sample path. We will repeat this operation, so let us write a function for it.
In [15]:
def time_series_dist(mc, init, ts_length=100):
"""
Return a distribution of visits by a sample path of length ts_length
of mc with an initial state init.
"""
X = mc.simulate(init=init, sample_size=ts_length)
bins = np.arange(mc.n+1)
hist, bin_edges = np.histogram(X, bins=bins)
dist = hist/len(X)
return dist
Here is a frequency distribution along a sample path from initial state 1:
In [16]:
time_series_dist(mc0, init=1)
Out[16]:
Let us visualize the distribution.
In [17]:
def draw_histogram(distribution,
ax=None, title=None, xlabel=None, ylabel=None, ylim=(0, 1)):
"""
Plot the given distribution.
"""
if ax is None:
fig, ax = plt.subplots()
n = len(distribution)
ax.bar(np.arange(n), distribution, align='center')
ax.set_xlim(-0.5, (n-1)+0.5)
ax.set_ylim(*ylim)
if title:
ax.set_title(title)
if xlabel:
ax.set_xlabel(xlabel)
if ylabel:
ax.set_ylabel(ylabel)
if ax is None:
plt.show()
In [18]:
%matplotlib inline
import matplotlib.pyplot as plt
In [19]:
init = 1
draw_histogram(time_series_dist(mc0, init=init),
title='Time series distribution with init={0}'.format(init),
xlabel='States')
plt.show()
Observe that the distribution is close to the stationary distribution
[0, 1/3, 0, 0, 2/3, 0].
Start with state 2:
In [20]:
init = 2
draw_histogram(time_series_dist(mc0, init=init),
title='Time series distribution with init={0}'.format(init),
xlabel='States')
plt.show()
Run the above (or below) cell several times; you will observe that the distribution differs across sample paths. Sometimes the state is absorbed into the recurrent class $\{0\}$, while other times it is absorbed into the recurrent class $\{1, 4\}$.
In [21]:
init = 2
fig, axes = plt.subplots(1, 2, figsize=(8, 3))
titles = ['Some sample path', 'Another sample path']
titles = [title + ' init={0}'.format(init) for title in titles]
for ax, title in zip(axes, titles):
draw_histogram(time_series_dist(mc0, init=init), ax=ax, title=title, xlabel='States')
fig.suptitle('Time series distributions', y=-0.05, fontsize=12)
plt.show()
Let us repeat the simulation for many times (say, 100 times)
and obtain the distribution of visits to each states
at a given time period T.
That is, we want to simulate the marginal distribution at time T.
In [22]:
def cross_sectional_dist(mc, init, T, num_reps=100):
x = np.empty(num_reps, dtype=int)
for i in range(num_reps):
x[i] = mc.simulate(init=init, sample_size=T+1)[-1]
bins = np.arange(mc.n+1)
hist, bin_edges = np.histogram(x, bins=bins)
dist = hist/len(x)
return dist
In [23]:
init = 2
T = 100
draw_histogram(cross_sectional_dist(mc0, init=init, T=T),
title='Empirical marginal distribution ' + \
'at T={0} with init={1}'.format(T, init))
plt.show()
Observe that the distribution is close to a convex combination of
the stationary distributions [1, 0, 0, 0, 0, 0] and [0, 1/3, 0, 0, 2/3, 0].
Let us simulate with the remaining states, 3, 4, and 5.
In [24]:
inits = [3, 4, 5]
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
for init, ax in zip(inits, axes):
draw_histogram(time_series_dist(mc0, init=init), ax=ax,
title='Initial state = {0}'.format(init),
xlabel='States')
fig.suptitle('Time series distributions', y=-0.05, fontsize=12)
plt.show()
Plot empirical marginal distributions at T=100 with initial states 3, 4, and 5.
In [25]:
# Write your own code
inits = [3, 4, 5]
T = 100
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
for init, ax in zip(inits, axes):
draw_histogram(cross_sectional_dist(mc0, init=init, T=T), ax=ax,
title='Initial state = {0}'.format(init),
xlabel='States')
fig.suptitle('Empirical marginal distributions at T={0}'.format(T), y=-0.05, fontsize=12)
plt.show()
The marginal distrubtions at time $T$ are obtained by $P^T$.
In [26]:
np.set_printoptions(suppress=True) # Suppress printing with floating point notation
T = 10
print 'P^{0} ='.format(T)
print np.linalg.matrix_power(mc0.P, T)
In [27]:
T = 100
print 'P^{0} ='.format(T)
print np.linalg.matrix_power(mc0.P, T)
Compare the rows with the stationary distributions obtained by mc0.stationary_distributions.
Note that mc0 is aperiodic
(i.e., the least common multiple of the periods of the recurrent class is one),
so that $P^T$ converges as $T \to \infty$.
In [28]:
mc0.period
Out[28]:
Consider the Markov chain given by the following stochastic matrix (where the actual values of non-zero probabilities are not important):
In [29]:
P = np.zeros((10, 10))
P[0, 3] = 1
P[1, [0, 4]] = 1/2
P[2, 6] = 1
P[3, [1, 2, 7]] = 1/3
P[4, 3] = 1
P[5, 4] = 1
P[6, 3] = 1
P[7, [6, 8]] = 1/2
P[8, 9] = 1
P[9, 5] = 1
In [30]:
np.set_printoptions(precision=3) # Reduce the number of digits printed
print P
In [31]:
mc1 = qe.MarkovChain(P)
We call the states $0, \ldots, 9$, respectively, instead of $1, \ldots, 10$.
(a) Check that this Markov chain is irreducible.
In [32]:
mc1.is_irreducible
Out[32]:
(c) Determine the period of this Markov chain.
In [33]:
mc1.period
Out[33]:
(d) Identify the cyclic classes.
In [34]:
mc1.cyclic_classes
Out[34]:
Let us discuss this exercise using the Markov chain from Exercise 9.
In [35]:
print mc1.P
Denote the period of $P$ by $d$:
In [36]:
d = mc1.period
Reorder the states so that the transition matrix is in block form:
In [37]:
permutation = np.concatenate(mc1.cyclic_classes)
In [38]:
P = mc1.P[permutation, :][:, permutation]
print P
Let us create a new MarkovChain instance:
In [39]:
mc2 = qe.MarkovChain(P)
Obtain the block components $P_0, \ldots, P_d$:
In [40]:
P_blocks = []
for i in range(d):
P_blocks.append(mc2.P[mc2.cyclic_classes[i%d], :][:, mc2.cyclic_classes[(i+1)%d]])
print 'P_{0} ='.format(i)
print P_blocks[i]
(b) Show that $P^d$ is block diagonal.
Hint: You may use np.linalg.matrix_power.
In [41]:
# Write your own code
P_power_d = np.linalg.matrix_power(mc2.P, d)
print P_power_d
(c) Verify that the $i$th diagonal block of $P^d$ equals $P_i P_{i+1} \ldots P_{d-1} P_0 \ldots P_{i-1}$.
Store these diagonal blocks in a list called P_power_d_blocks.
In [42]:
# Write your own code
P_power_d_blocks = []
ordinals = ['0th', '1st', '2nd']
for i in range(d):
P_power_d_blocks.append(P_power_d[mc2.cyclic_classes[i], :][:, mc2.cyclic_classes[i]])
print '{0} diagonal block of P^d ='.format(ordinals[i])
print P_power_d_blocks[i]
Print $P_i P_{i+1} \ldots P_{d-1} P_0 \ldots P_{i-1}$ for each $i$.
In [43]:
for i in range(d):
R = np.eye(P_blocks[i].shape[0])
string = ''
for j in range(d):
R = R.dot(P_blocks[(i+j)%d])
string += 'P_{0} '.format((i+j)%d)
P_power_d_blocks.append(R)
print string + '='
print R
(d) Obtain the stationary distributions each associated with the diagonal blocks of $P^d$.
Store them in a list called pi_s.
In [44]:
# Write your own code
pi_s = []
for i in range(d):
pi_s.append(qe.MarkovChain(P_power_d_blocks[i]).stationary_distributions[0])
print 'pi^{0} ='.format(i)
print pi_s[i]
Verify that $\pi^{i+1} = \pi^i P_i$.
In [45]:
# Write your own code
for i in range(d):
print 'pi^{0} P_{0} ='.format(i)
print np.dot(pi_s[i], P_blocks[i])
Obtain the unique stationary distribution $\pi$ of the original Markov chain.
In [46]:
print mc2.stationary_distributions[0]
Verify that $\pi = (\pi^0, \ldots, \pi^d) / d$.
In [47]:
# Write your own code
rhs = np.zeros(mc2.n)
for i in range(d):
rhs[mc2.cyclic_classes[i]] = pi_s[i]
rhs /= d
print rhs
(e)
In [48]:
# Write your answer in a Markdown mode, with providing code if necessary.
Print $P^1, P^2, \ldots, P^d$.
In [49]:
# Write your own code
for i in range(1, d+1):
print 'P^{0} ='.format(i)
print np.linalg.matrix_power(mc2.P, i)
Print $P^{2d}$, $P^{4d}$, and $P^{6d}$.
In [50]:
# Write your own code
for i in [2, 4, 6]:
print 'P^{0}d ='.format(i)
print np.linalg.matrix_power(mc2.P, i*d)
Print $P^{10d + 1}, \ldots, P^{10d + d}$.
In [51]:
# Write your own code
for i in range(10*d+1, 10*d+1+d):
print 'P^{0}d ='.format(i)
print np.linalg.matrix_power(mc2.P, i)
Let us simulate the Markov chain mc2.
Plot the frequency distribution of visits to the states along a sample path starting at state 0:
In [52]:
init = 0
draw_histogram(time_series_dist(mc2, init=init),
title='Time series distribution with init={0}'.format(init),
xlabel='States', ylim=(0, 0.35))
plt.show()
Observe that the distribution is close to the (unique) stationary distribution.
Next, plot the simulated marginal distributions at $T = 10d + 1, \ldots, 11d, 11d + 1, \ldots, 12d$ with initial state 0:
In [53]:
init = 0
T = 10 * d + 1
fig, axes = plt.subplots(2, d, figsize=(12, 6))
for t, ax in enumerate(axes.flatten()):
draw_histogram(cross_sectional_dist(mc2, init=init, T=T+t), ax=ax,
title='T = {0}'.format(T+t))
fig.suptitle('Empirical marginal distributions with init={0}'.format(init), y=0.05, fontsize=12)
plt.show()
Compare the rows of $P^{10d + 1}, \ldots, P^{10d + d}$.
In [53]: