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__

```
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
```

`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)

`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)

`1`

and `4`

with certain frequencies.

```
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()

```
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()

`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()

`[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`

.

`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

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

```
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 [ ]:
```