<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Finite Markov Chains
Step1:  Example 1
Step2:  Create a MarkovChain instance
Step3:  Classification of states
This Markov chain is reducible
Step4:  Determine the communication classes
Step5:  Classify the states of this Markov chain
Step6:  Obtain a list of the recurrent states
Step7:  Obtain a list of the transient states
Step8:  A Markov chain is reducible (i.e., its directed graph is not strongly connected)
if and only if by symmetric permulations of rows and columns,
its transition probability matrix is written in the form ("canonical form")
$$
\begin{pmatrix}
U & 0 \
W & V
\end{pmatrix},
$$
where $U$ and $W$ are square matrices.
Such a form for mc1 is obtained by the following
Step9:  This Markov chain is aperiodic
(i.e., the least common multiple of the periods of the recurrent sub-chains is one)
Step10:  Indeed, each of the sub-chains corresponding to the recurrent classes has period $1$,
i.e., every recurrent state is aperiodic
Step11:  Stationary distributions
For each recurrent class $C$, there is a unique stationary distribution $\psi^C$
such that $\psi^C_i > 0$ for all $i \in C$ and $\psi^C_i = 0$ otherwise.
MarkovChain.stationary_distributions returns
these unique stationary distributions for the recurrent classes.
Any stationary distribution is written as a convex combination of these distributions.
Step12:  These are indeed stationary distributions
Step14:  Plot these distributions.
Step15:  Simulation
Let us simulate our Markov chain mc1.
The simualte method generates a sample path
of length given by the first argument, ts_length,
with an initial state as specified by an optional argument init;
if not specified, the initial state is randomly drawn. 
A sample path from state 0
Step16:  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
Step17:  You can observe that the chain stays in the recurrent class ${1,4}$
and visits states 1 and 4 with certain frequencies.
If init is not specified, the initial state is randomly chosen
Step18:  Note on reproducibility
Step20:  Time series averages
Now, let us compute the frequency distribution along a sample path, given by
$$
\frac{1}{t} \sum_{\tau=0}^{t-1} \mathbf{1}{X_{\tau} = s}
\quad (s \in S).
$$
Step21:  Here is a frequency distribution along a sample path,  of length 100,
from initial state 1, which is a recurrent state
Step22:  Length 10,000
Step23:  The distribution becomes close to the stationary distribution (0, 1/3, 0, 0, 2/3, 0).
Plot the frequency distributions for a couple of different time lengths
Step24:  Start with state 2,
which is a transient state
Step25:  Run the above cell several times;
you will observe that the limit 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}$.
Step26:  In fact,
for almost every sample path of a finite Markov chain ${X_t}$,
for some recurrent class $C$ we have
$$
\frac{1}{t} \sum_{\tau=0}^{t-1} \mathbf{1}{X_{\tau} = s} \to \psi^C[s]
\quad \text{as $t \to \infty$}
$$
for all states $s$,
where $\psi^C$ is the stationary distribution associated with the recurrent class $C$.
If the initial state $s_0$ is a recurrent state,
then the recurrent class $C$ above is the one that contains $s_0$,
while if it is a transient state,
then the recurrent class to which the convergence occurs depends on the sample path.
Let us simulate with the remaining states, 3, 4, and 5.
Step28:  Cross sectional averages
Next, let us repeat the simulation for many times (say, 10,000 times)
and obtain the distribution of visits to each state at a given time period T.
That is, we want to simulate the marginal distribution at time T.
Step29:  Start with state 1
Step30:  The distribution is close to the stationary distribution (0, 1/3, 0, 0, 2/3, 0).
Plot the simulated marginal distribution at T for some values of T.
Step31:  Starting with a transient state 2
Step32:  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),
which is a stationary distribution itself.
How the simulated marginal distribution evolves
Step33:  Since our Markov chain is aperiodic (i.e., every recurrent class is aperiodic),
the marginal disribution at time $T$ converges as $T \to \infty$
to some stationary distribution,
and the limit distribution depends on the initial state,
according to the probabilities that the state is absorbed into the recurrent classes.
For initial states 3, 4, and 5
Step34:  Powers of $P$
The marginal distributions at time $T$ are obtained by $P^T$.
Step35:  In the canonical form
Step36:  Observe that the first three rows, which correspond to the recurrent states,
are close to the stationary distributions associated with the corresponding recurrent classes.
Example 2
Step37:  This Markov chain is irreducible
Step38:  This Markov chain is periodic
Step39:  Its period, which we denote by $d$
Step40:  Identify the cyclic classes
Step41:  Cyclic normal form
If a Markov chain is periodic with period $d \geq 2$,
then its transition probability matrix is written in the form ("cyclic normal form")
$$
\begin{pmatrix}
0       & P_0    & 0      & 0      & \cdots & 0 \
0       & 0      & P_1    & 0      & \cdots & 0 \
0       & 0      & 0      & P_2    & \cdots & 0 \
\vdots  & \vdots & \vdots & \vdots & \ddots & \vdots \
0       & 0      & 0      & 0      & \cdots & P_{d-2} \
P_{d-1} & 0      & 0      & 0      & \cdots & 0
\end{pmatrix}.
$$
Represent our Markov chain in cyclic normal form
Step42:  Re-define the Markov chain with the above matrix Q
Step43:  Obtain the block components $P_0, \cdots, P_{d-1}$
Step44:  $P^d$ is block diagonal
Step45:  The $i$th diagonal block of $P^d$ equals $P_i P_{i+1} \cdots P_{d-1} P_0 \cdots P_{i-1}$
Step46:  Stationary distributions
The Markov chain mc2 has a unique stationary distribution,
which we denote by $\pi$
Step47:  Obtain the stationary distributions $\pi^0, \ldots, \pi^{d-1}$
each associated with the diagonal blocks of $P^d$
Step48:  Verify that $\pi^{i+1} = \pi^i P_i$
Step49:  Verify that $\pi = (\pi^0 + \cdots + \pi^{d-1})/d$
Step50:  Powers of $P$
Since the Markov chain in consideration is periodic,
the marginal distribution does not converge, but changes periodically.
Let us compute the powers of the transition probability matrix (in cyclic normal form)
Step51:  Print $P^1, P^2, \ldots, P^d$
Step52:  Print $P^{2d}$, $P^{4d}$, and $P^{6d}$
Step53:  $P^{kd}$ converges as $k \to \infty$ to a matrix that contains $\pi^0, \ldots, \pi^{d-1}$.
Print $P^{kd+1}, \ldots, P^{kd+d}$ with $k = 10$ for example
Step54:  But $P^i$ itself does not converge.
Simulation
Plot the frequency distribution of visits to the states
along a sample path starting at state 0
Step55:  Observe that the distribution is close to the (unique) stationary distribution $\pi$.
Step56:  Next, plot the simulated marginal distributions
at $T = 10d+1, \ldots, 11d, 11d+1, \ldots, 12d$ with initial state 0
Step57:  Compare these with the rows of $P^{10d+1}, \ldots, P^{10d+d}$.
Example 3
Step58:  If $\varepsilon = 0$,
then the Markovh chain is reducible into two recurrent classes, [0, 1] and [2]
Step59:  If $\varepsilon > 0$ but small, the chain is irreducible,
but transition within each of the subsets [0, 1] and [2] is much more likely
than that between these sets.
Step60:  Analytically, the unique stationary distribution of the chain with $\varepsilon > 0$
is (1/3, 1/3, 1/3), independent of the value of $\varepsilon$.
However,
for such matrices with small values of $\varepsilon > 0$,
general purpose eigenvalue solvers are numerically unstable.
For example, if we use numpy.linalg.eig
to compute the eigenvector that corresponds
to the dominant (i.e., largest in magnitude) eigenvalue
Step61:  The output becomes farther from the actual stationary distribution (1/3, 1/3, 1/3)
as $\varepsilon$ becomes smaller.
The same applies to scipy.linalg.eig
Step62:  MarkovChain in quantecon employs
the algorithm called the "GTH algorithm",
which is a numerically stable variant of Gaussian elimination,
specialized for Markov chains.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
from quantecon.markov import MarkovChain
Explanation: Finite Markov Chains: Examples
Daisuke Oyama
Faculty of Economics, University of Tokyo
This notebook demonstrates how to analyze finite-state Markov chains
with the MarkovChain class.
For basic concepts and properties on Markov chains, see
the lecture on finite Markov chains
  in Quantitative Economics, and
the documentation for MarkovChain.
For algorithmic issues in detecting reducibility and periodicity of a Markov chain,
see, for example,
J. P. Jarvis and D. R. Shier,
  "Graph-Theoretic Analysis of Finite Markov Chains,"
from which we draw some examples below.
End of explanation
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
print(P)
Explanation: Example 1: Reducible chain
Consider the Markov chain given by the following stochastic matrix,
taken from Exercise 3 in Jarvis and Shier
(where the actual values of non-zero probabilities are not important):
End of explanation
mc1 = MarkovChain(P)
Explanation: Create a MarkovChain instance:
End of explanation
mc1.is_irreducible
mc1.num_communication_classes
Explanation: Classification of states
This Markov chain is reducible:
End of explanation
mc1.communication_classes
Explanation: Determine the communication classes:
End of explanation
mc1.recurrent_classes
Explanation: Classify the states of this Markov chain:
End of explanation
recurrent_states = np.concatenate(mc1.recurrent_classes)
print(recurrent_states)
Explanation: Obtain a list of the recurrent states:
End of explanation
transient_states = np.setdiff1d(np.arange(mc1.n), recurrent_states)
print(transient_states)
Explanation: Obtain a list of the transient states:
End of explanation
permutation = np.concatenate([recurrent_states, transient_states])
print(mc1.P[permutation, :][:, permutation])
Explanation: A Markov chain is reducible (i.e., its directed graph is not strongly connected)
if and only if by symmetric permulations of rows and columns,
its transition probability matrix is written in the form ("canonical form")
$$
\begin{pmatrix}
U & 0 \
W & V
\end{pmatrix},
$$
where $U$ and $W$ are square matrices.
Such a form for mc1 is obtained by the following:
End of explanation
mc1.is_aperiodic
Explanation: This Markov chain is aperiodic
(i.e., the least common multiple of the periods of the recurrent sub-chains is one):
End of explanation
for recurrent_class in mc1.recurrent_classes:
    sub_matrix = P[recurrent_class, :][:, recurrent_class]
    d = MarkovChain(sub_matrix).period
    print('Period of the sub-chain\n{0}\n = {1}'.format(sub_matrix, d))
Explanation: Indeed, each of the sub-chains corresponding to the recurrent classes has period $1$,
i.e., every recurrent state is aperiodic:
End of explanation
print(mc1.stationary_distributions)
Explanation: Stationary distributions
For each recurrent class $C$, there is a unique stationary distribution $\psi^C$
such that $\psi^C_i > 0$ for all $i \in C$ and $\psi^C_i = 0$ otherwise.
MarkovChain.stationary_distributions returns
these unique stationary distributions for the recurrent classes.
Any stationary distribution is written as a convex combination of these distributions.
End of explanation
print(mc1.stationary_distributions.dot(mc1.P))
Explanation: These are indeed stationary distributions:
End of explanation
def draw_histogram(distribution, ax=None, figsize=None,
                   title=None, xlabel=None, ylabel=None, ylim=(0, 1)):
    
    Plot the given distribution.
    
    
    if ax is None:
        fig, ax = plt.subplots(figsize=figsize)
    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()
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
titles = ['Stationary distribution for the recurrent class {0}'.format(recurrent_class)
          for recurrent_class in mc1.recurrent_classes]
for ax, title, dist in zip(axes, titles, mc1.stationary_distributions):
    draw_histogram(dist, ax=ax, title=title, xlabel='States')
fig.suptitle('Stationary distributions', y=-0.05, fontsize=12)
plt.show()
Explanation: Plot these distributions.
End of explanation
mc1.simulate(50, init=0)
Explanation: Simulation
Let us simulate our Markov chain mc1.
The simualte method generates a sample path
of length given by the first argument, ts_length,
with an initial state as specified by an optional argument init;
if not specified, the initial state is randomly drawn. 
A sample path from state 0:
End of explanation
mc1.simulate(50, init=1)
Explanation: 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:
End of explanation
mc1.simulate(50)
Explanation: You can observe that the chain stays in the recurrent class ${1,4}$
and visits states 1 and 4 with certain frequencies.
If init is not specified, the initial state is randomly chosen:
End of explanation
mc1.simulate(50, random_state=12345)
Explanation: Note on reproducibility:
The simulate method offers an option random_state to set a random seed
to initialize the pseudo-random number generator.
As you provide the same random seed value,
simulate returns the same outcome.
For example, the following will always give the same sequence:
End of explanation
def time_series_dist(mc, t, init=None, random_state=None):
    
    Return the distribution of visits by a sample path of length t
    of mc with an initial state init.
    
    
    t_max = np.max(t)
    dim = 1
    try:
        ts_size = len(t)  # t is an array
        ts_array = t
        dim = 2
    except:  # t is an int
        ts_size = 1
        ts_array = [t]
        
    X = mc.simulate(ts_length=t_max, init=init, random_state=random_state)
    dists = np.empty((ts_size, mc.n))
    bins = np.arange(mc.n+1)
    for i, length in enumerate(ts_array):
        hist, bin_edges = np.histogram(X[:length], bins=bins)
        dists[i, :] = hist / length
    if dim == 1:
        return dists[0]
    else:
        return dists
Explanation: Time series averages
Now, let us compute the frequency distribution along a sample path, given by
$$
\frac{1}{t} \sum_{\tau=0}^{t-1} \mathbf{1}{X_{\tau} = s}
\quad (s \in S).
$$
End of explanation
time_series_dist(mc1, t=100, init=1)
Explanation: Here is a frequency distribution along a sample path,  of length 100,
from initial state 1, which is a recurrent state:
End of explanation
time_series_dist(mc1, t=10**4, init=1)
Explanation: Length 10,000:
End of explanation
def plot_time_series_dists(mc, init, ts, seed=None, figsize=(12,4)):
    dists = time_series_dist(mc, t=ts, init=init, random_state=seed)
    fig, axes = plt.subplots(1, len(ts), figsize=figsize)
    titles = ['t={0}'.format(t) for t in ts]
    for ax, title, dist in zip(axes, titles, dists):
        draw_histogram(dist, ax=ax, title=title, xlabel='States')
    fig.suptitle('Time series distributions with init={0}'.format(init),
                 y=-0.05, fontsize=12)
    plt.show()
init = 1
ts = [5, 10, 50, 100]
plot_time_series_dists(mc1, init, ts)
Explanation: The distribution becomes close to the stationary distribution (0, 1/3, 0, 0, 2/3, 0).
Plot the frequency distributions for a couple of different time lengths:
End of explanation
init = 2
ts = [5, 10, 50, 100]
plot_time_series_dists(mc1, init, ts)
Explanation: Start with state 2,
which is a transient state:
End of explanation
init = 2
ts = [5, 10, 50, 100]
seeds = [222, 2222]
descriptions = ['{0} sample path with init={1}'.format(adjective, init) for adjective
                in ['Some'] + ['Another'] + ['Yet another']*(len(seeds)-1)]
for seed, description in zip(seeds, descriptions):
    print(description)
    plot_time_series_dists(mc1, init, ts, seed=seed)
Explanation: Run the above cell several times;
you will observe that the limit 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}$.
End of explanation
inits = [3, 4, 5]
t = 100
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
for init, ax in zip(inits, axes):
    draw_histogram(time_series_dist(mc1, t=t, init=init), ax=ax,
                   title='Initial state = {0}'.format(init),
                   xlabel='States')
fig.suptitle('Time series distributions for t={0}'.format(t),
             y=-0.05, fontsize=12)
plt.show()
Explanation: In fact,
for almost every sample path of a finite Markov chain ${X_t}$,
for some recurrent class $C$ we have
$$
\frac{1}{t} \sum_{\tau=0}^{t-1} \mathbf{1}{X_{\tau} = s} \to \psi^C[s]
\quad \text{as $t \to \infty$}
$$
for all states $s$,
where $\psi^C$ is the stationary distribution associated with the recurrent class $C$.
If the initial state $s_0$ is a recurrent state,
then the recurrent class $C$ above is the one that contains $s_0$,
while if it is a transient state,
then the recurrent class to which the convergence occurs depends on the sample path.
Let us simulate with the remaining states, 3, 4, and 5.
End of explanation
def cross_sectional_dist(mc, T, init=None, num_reps=10**4, random_state=None):
    
    Return the distribution of visits at time T by num_reps times of simulation
    of mc with an initial state init.
    
    
    T_max = np.max(T)
    dim = 1
    try:
        Ts_size = len(T)  # T is an array
        Ts_array = T
        dim = 2
    except:  # T is an int
        Ts_size = 1
        Ts_array = [T]
    
    x = mc.simulate(ts_length=T_max+1, init=init, num_reps=num_reps,
                    random_state=random_state)[:, Ts_array]
    dists = np.empty((x.shape[-1], mc.n))
    bins = np.arange(mc.n+1)
    for i in range(x.shape[-1]):
        hist, bin_edges = np.histogram(x[:, i], bins=bins)
        dists[i, :] = hist / num_reps
    if dim == 1:
        return dists[0]
    else:
        return dists
Explanation: Cross sectional averages
Next, let us repeat the simulation for many times (say, 10,000 times)
and obtain the distribution of visits to each state at a given time period T.
That is, we want to simulate the marginal distribution at time T.
End of explanation
init = 1
T = 10
cross_sectional_dist(mc1, init=init, T=T)
T = 100
cross_sectional_dist(mc1, init=init, T=T)
Explanation: Start with state 1:
End of explanation
def plot_cross_sectional_dists(mc, init, Ts, num_reps=10**4, seed=None,
                               figsize=(12,4)):
    dists = cross_sectional_dist(mc, T=Ts, init=init, num_reps=num_reps,
                                 random_state=seed)
    fig, axes = plt.subplots(1, len(Ts), figsize=figsize)
    titles = ['T={0}'.format(T) for T in Ts]
    for ax, title, dist in zip(axes, titles, dists):
        draw_histogram(dist, ax=ax, title=title, xlabel='States')
    fig.suptitle('Cross sectional distributions with init={0}'.format(init),
                 y=-0.05, fontsize=12)
    plt.show()
init = 1
Ts = [2, 3, 5, 10]
plot_cross_sectional_dists(mc1, init, Ts)
Explanation: The distribution is close to the stationary distribution (0, 1/3, 0, 0, 2/3, 0).
Plot the simulated marginal distribution at T for some values of T.
End of explanation
init = 2
T = 10
cross_sectional_dist(mc1, init=init, T=T)
T = 100
dist = cross_sectional_dist(mc1, init=init, T=T)
dist
draw_histogram(dist,
               title='Cross sectional distribution at T={T} with init={init}'
                     .format(T=T, init=init),
               xlabel='States')
Explanation: Starting with a transient state 2:
End of explanation
init = 2
Ts = [2, 3, 5, 10]
plot_cross_sectional_dists(mc1, init, Ts)
Explanation: 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),
which is a stationary distribution itself.
How the simulated marginal distribution evolves:
End of explanation
inits = [3, 4, 5]
T = 10
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
for init, ax in zip(inits, axes):
    draw_histogram(cross_sectional_dist(mc1, T=T, init=init), ax=ax,
                   title='Initial state = {0}'.format(init),
                   xlabel='States')
fig.suptitle('Cross sectional distribution at T={0}'.format(T),
             y=-0.05, fontsize=12)
plt.show()
Explanation: Since our Markov chain is aperiodic (i.e., every recurrent class is aperiodic),
the marginal disribution at time $T$ converges as $T \to \infty$
to some stationary distribution,
and the limit distribution depends on the initial state,
according to the probabilities that the state is absorbed into the recurrent classes.
For initial states 3, 4, and 5:
End of explanation
np.set_printoptions(suppress=True)  # Suppress printing with floating point notation
Ts = [10, 20, 30]
for T in Ts:
    print('P^{T} =\n{P_T}'.format(T=T, P_T=np.linalg.matrix_power(mc1.P, T)))
Explanation: Powers of $P$
The marginal distributions at time $T$ are obtained by $P^T$.
End of explanation
Q = mc1.P[permutation, :][:, permutation]
print('Q =\n{Q}'.format(Q=Q))
for T in Ts:
    print('Q^{T} =\n{Q_T}'.format(T=T, Q_T=np.linalg.matrix_power(Q, T)))
Explanation: In the canonical form:
End of explanation
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
np.set_printoptions(precision=3)  # Reduce the number of digits printed
print(P)
mc2 = MarkovChain(P)
Explanation: Observe that the first three rows, which correspond to the recurrent states,
are close to the stationary distributions associated with the corresponding recurrent classes.
Example 2: Periodic chain
Consider the Markov chain given by the following stochastic matrix,
taken from Exercise 9 (see also Exercise 11) in Jarvis and Shier
(where the actual values of non-zero probabilities are not important):
End of explanation
mc2.is_irreducible
Explanation: This Markov chain is irreducible:
End of explanation
mc2.is_aperiodic
Explanation: This Markov chain is periodic:
End of explanation
d = mc2.period
print(d)
Explanation: Its period, which we denote by $d$:
End of explanation
mc2.cyclic_classes
Explanation: Identify the cyclic classes:
End of explanation
permutation = np.concatenate(mc2.cyclic_classes)
Q = mc2.P[permutation, :][:, permutation]
print(Q)
Explanation: Cyclic normal form
If a Markov chain is periodic with period $d \geq 2$,
then its transition probability matrix is written in the form ("cyclic normal form")
$$
\begin{pmatrix}
0       & P_0    & 0      & 0      & \cdots & 0 \
0       & 0      & P_1    & 0      & \cdots & 0 \
0       & 0      & 0      & P_2    & \cdots & 0 \
\vdots  & \vdots & \vdots & \vdots & \ddots & \vdots \
0       & 0      & 0      & 0      & \cdots & P_{d-2} \
P_{d-1} & 0      & 0      & 0      & \cdots & 0
\end{pmatrix}.
$$
Represent our Markov chain in cyclic normal form:
End of explanation
mc2 = MarkovChain(Q)
Explanation: Re-define the Markov chain with the above matrix Q:
End of explanation
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_{i} =\n{P_block}'.format(i=i, P_block=P_blocks[i]))
Explanation: Obtain the block components $P_0, \cdots, P_{d-1}$:
End of explanation
P_power_d = np.linalg.matrix_power(mc2.P, d)
print(P_power_d)
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('{ordinal} diagonal block of P^d =\n{P_power_d_block}'
          .format(ordinal=ordinals[i], P_power_d_block=P_power_d_blocks[i]))
Explanation: $P^d$ is block diagonal:
End of explanation
products = []
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)
    products.append(R)
    print(string + '=\n{R}'.format(R=R))
for matrix0, matrix1 in zip(P_power_d_blocks, products):
    print(np.array_equal(matrix0, matrix1))
Explanation: The $i$th diagonal block of $P^d$ equals $P_i P_{i+1} \cdots P_{d-1} P_0 \cdots P_{i-1}$:
End of explanation
len(mc2.stationary_distributions)
pi = mc2.stationary_distributions[0]
print(pi)
draw_histogram(pi,
               title='Stationary distribution', xlabel='States', ylim=(0, 0.35))
Explanation: Stationary distributions
The Markov chain mc2 has a unique stationary distribution,
which we denote by $\pi$:
End of explanation
pi_s = []
for i in range(d):
    pi_s.append(MarkovChain(P_power_d_blocks[i]).stationary_distributions[0])
    print('pi^{i} =\n{pi_i}'.format(i=i, pi_i=pi_s[i]))
fig, axes = plt.subplots(1, d, figsize=(12, 3))
for i, ax in enumerate(axes):
    pi_i_full_dim = np.zeros(mc2.n)
    pi_i_full_dim[mc2.cyclic_classes[i]] = pi_s[i]
    draw_histogram(pi_i_full_dim,
                   ax=ax, title='$\pi^{i}$'.format(i=i), xlabel='States')
fig.suptitle('Stationary distributions for the diagonal blocks',
             y=-0.05, fontsize=12)
plt.show()
Explanation: Obtain the stationary distributions $\pi^0, \ldots, \pi^{d-1}$
each associated with the diagonal blocks of $P^d$:
End of explanation
for i in range(d):
    print('pi^{i} P_{i} =\n{dot}'.format(i=i, dot=np.dot(pi_s[i], P_blocks[i])))
Explanation: Verify that $\pi^{i+1} = \pi^i P_i$:
End of explanation
# Right hand side of the above identity
rhs = np.zeros(mc2.n)
for i in range(d):
    rhs[mc2.cyclic_classes[i]] = pi_s[i]
rhs /= d
print(rhs)
np.allclose(pi, rhs)
Explanation: Verify that $\pi = (\pi^0 + \cdots + \pi^{d-1})/d$:
End of explanation
np.set_printoptions(suppress=True)  # Suppress printing with floating point notation
Explanation: Powers of $P$
Since the Markov chain in consideration is periodic,
the marginal distribution does not converge, but changes periodically.
Let us compute the powers of the transition probability matrix (in cyclic normal form):
End of explanation
for i in range(1, d+1):
    print('P^{i} =\n{P_i}'.format(i=i, P_i=np.linalg.matrix_power(mc2.P, i)))
Explanation: Print $P^1, P^2, \ldots, P^d$:
End of explanation
for i in [k*d for k in [2, 4, 6]]:
    print('P^{i} =\n{P_i}'.format(i=i, P_i=np.linalg.matrix_power(mc2.P, i)))
Explanation: Print $P^{2d}$, $P^{4d}$, and $P^{6d}$:
End of explanation
for i in range(10*d+1, 10*d+1+d):
    print('P^{i} =\n{P_i}'.format(i=i, P_i=np.linalg.matrix_power(mc2.P, i)))
Explanation: $P^{kd}$ converges as $k \to \infty$ to a matrix that contains $\pi^0, \ldots, \pi^{d-1}$.
Print $P^{kd+1}, \ldots, P^{kd+d}$ with $k = 10$ for example:
End of explanation
init = 0
dist = time_series_dist(mc2, init=init, t=10**4)
print(dist)
draw_histogram(dist,
               title='Time series distribution with init={0}'.format(init),
               xlabel='States', ylim=(0, 0.35))
plt.show()
Explanation: But $P^i$ itself does not converge.
Simulation
Plot the frequency distribution of visits to the states
along a sample path starting at state 0:
End of explanation
print(pi)
Explanation: Observe that the distribution is close to the (unique) stationary distribution $\pi$.
End of explanation
init = 0
k = 10
Ts = [k*d + 1 + i for i in range(2*d)]
num_reps = 10**2
dists = cross_sectional_dist(mc2, T=Ts, init=init, num_reps=num_reps)
fig, axes = plt.subplots(2, d, figsize=(12, 6))
for dist, T, ax in zip(dists, Ts, axes.flatten()):
    draw_histogram(dist, ax=ax,
                   title='T = {T}'.format(T=T))
fig.suptitle('Cross sectional distributions with init={init}'.format(init=init),
             y=0.05, fontsize=12)
plt.show()
Explanation: Next, plot the simulated marginal distributions
at $T = 10d+1, \ldots, 11d, 11d+1, \ldots, 12d$ with initial state 0:
End of explanation
def P_epsilon(eps, p=0.5):
    P = np.array([[1-(p+eps), p,         eps],
                  [p,         1-(p+eps), eps],
                  [eps,       eps,       1-2*eps]])
    return P
Explanation: Compare these with the rows of $P^{10d+1}, \ldots, P^{10d+d}$.
Example 3: Nearly completely decomposable chain
Consider the Markov chain given by the following stochastic matrix $P^{\varepsilon}$,
parameterized by $\varepsilon$:
End of explanation
P_epsilon(0)
MarkovChain(P_epsilon(0)).recurrent_classes
Explanation: If $\varepsilon = 0$,
then the Markovh chain is reducible into two recurrent classes, [0, 1] and [2]:
End of explanation
P_epsilon(0.001)
MarkovChain(P_epsilon(0.001)).recurrent_classes
Explanation: If $\varepsilon > 0$ but small, the chain is irreducible,
but transition within each of the subsets [0, 1] and [2] is much more likely
than that between these sets.
End of explanation
np.set_printoptions(precision=15)  # Increase the number of digits printed
epsilons = [float('1e-{i}'.format(i=i)) for i in range(12, 18)]
for eps in epsilons:
    print('epsilon = {eps}'.format(eps=eps))
    w, v = np.linalg.eig(P_epsilon(eps).T)
    i = w.argmax()
    print(v[:, i]/v[:, i].sum())
Explanation: Analytically, the unique stationary distribution of the chain with $\varepsilon > 0$
is (1/3, 1/3, 1/3), independent of the value of $\varepsilon$.
However,
for such matrices with small values of $\varepsilon > 0$,
general purpose eigenvalue solvers are numerically unstable.
For example, if we use numpy.linalg.eig
to compute the eigenvector that corresponds
to the dominant (i.e., largest in magnitude) eigenvalue:
End of explanation
import scipy.linalg
for eps in epsilons:
    print('epsilon = {eps}'.format(eps=eps))
    w, v = scipy.linalg.eig(P_epsilon(eps), left=True, right=False)
    i = w.argmax()
    print(v[:, i]/v[:, i].sum())
Explanation: The output becomes farther from the actual stationary distribution (1/3, 1/3, 1/3)
as $\varepsilon$ becomes smaller.
The same applies to scipy.linalg.eig:
End of explanation
for eps in epsilons + [1e-100]:
    print('epsilon = {eps}'.format(eps=eps))
    print(MarkovChain(P_epsilon(eps)).stationary_distributions[0])
Explanation: MarkovChain in quantecon employs
the algorithm called the "GTH algorithm",
which is a numerically stable variant of Gaussian elimination,
specialized for Markov chains.
End of explanation
<END_TASK>