Reading: J. P. Jarvis and D. R. Shier, "Graph-Theoretic Analysis of Finite Markov Chain."
In [ ]:
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 [ ]:
qe.__version__
Consider the Markov chain given by the following stochastic matrix (where the actual values of non-zero probabilities are not important):
In [ ]:
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 [ ]:
print P
Create a MarkovChain instance:
In [ ]:
mc0 = qe.MarkovChain(P)
We call the states $0, \ldots, 5$, respectively, instead of $1, \ldots, 6$.
(a) Determine the communication classes.
In [ ]:
mc0.communication_classes
(b) Classify the states of this Markov chain.
In [ ]:
mc0.recurrent_classes
Obtain a list of the recurrent states.
In [ ]:
# Write your own code
Obtain a list of the transient states.
In [ ]:
# Write your own code
(c) Does the chain have a unique stationary distribution?
In [ ]:
mc0.num_recurrent_classes
Obtain the stationary distributions:
In [ ]:
vecs = mc0.stationary_distributions
print vecs
Verify that the above vectors are indeed stationary distributions.
Hint: Use np.dot.
In [ ]:
# Write your own code
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 [ ]:
mc0.simulate(init=0, sample_size=50)
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 [ ]:
mc0.simulate(init=1, sample_size=50)
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 [ ]:
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 [ ]:
time_series_dist(mc0, init=1)
Let us visualize the distribution.
In [ ]:
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 [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
In [ ]:
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 [ ]:
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 [ ]:
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 [ ]:
def cross_sectional_dist(mc, init, T, num_reps=100):
"""
Return a distribution of visits at time T across num_reps sample paths
of mc with an initial state init.
"""
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 [ ]:
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 [ ]:
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 [ ]:
# Write your own code
The marginal distrubtions at time $T$ are obtained by $P^T$.
In [ ]:
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 [ ]:
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 [ ]:
mc0.period
Consider the Markov chain given by the following stochastic matrix (where the actual values of non-zero probabilities are not important):
In [ ]:
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 [ ]:
np.set_printoptions(precision=3) # Reduce the number of digits printed
print P
In [ ]:
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 [ ]:
mc1.is_irreducible
(c) Determine the period of this Markov chain.
In [ ]:
mc1.period
(d) Identify the cyclic classes.
In [ ]:
mc1.cyclic_classes
Let us discuss this exercise using the Markov chain from Exercise 9.
In [ ]:
print mc1.P
Denote the period of $P$ by $d$:
In [ ]:
d = mc1.period
Reorder the states so that the transition matrix is in block form:
In [ ]:
permutation = np.concatenate(mc1.cyclic_classes)
In [ ]:
P = mc1.P[permutation, :][:, permutation]
print P
Let us create a new MarkovChain instance:
In [ ]:
mc2 = qe.MarkovChain(P)
Obtain the block components $P_0, \ldots, P_d$:
In [ ]:
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 [ ]:
# Write your own code
(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 [ ]:
# Write your own code
P_power_d_blocks = []
(d) Obtain the stationary distributions each associated with the diagonal blocks of $P^d$.
Store them in a list called pi_s.
In [ ]:
# Write your own code
pi_s = []
Verify that $\pi^{i+1} = \pi^i P_i$.
In [ ]:
# Write your own code
Obtain the unique stationary distribution $\pi$ of the original Markov chain.
In [ ]:
print mc2.stationary_distributions[0]
Verify that $\pi = (\pi^0, \ldots, \pi^d) / d$.
In [ ]:
# Write your own code
(e)
In [ ]:
# Write your answer in a Markdown mode, with providing code if necessary.
Print $P^1, P^2, \ldots, P^d$.
In [ ]:
# Write your own code
Print $P^{2d}$, $P^{4d}$, and $P^{6d}$.
In [ ]:
# Write your own code
Print $P^{10d + 1}, \ldots, P^{10d + d}$.
In [ ]:
# Write your own code
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 [ ]:
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 [ ]:
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 [ ]: