diff --git a/lectures/markov_chains_I.md b/lectures/markov_chains_I.md index c50d472e..2c308ce2 100644 --- a/lectures/markov_chains_I.md +++ b/lectures/markov_chains_I.md @@ -61,6 +61,8 @@ import matplotlib as mpl from mpl_toolkits.mplot3d import Axes3D from matplotlib.animation import FuncAnimation from IPython.display import HTML +from matplotlib.patches import Polygon +from mpl_toolkits.mplot3d.art3d import Poly3DCollection ``` ## Definitions and examples @@ -743,6 +745,11 @@ This is, in some sense, a steady state probability of unemployment. Not surprisingly it tends to zero as $\beta \to 0$, and to one as $\alpha \to 0$. + + + + + ### Calculating stationary distributions A stable algorithm for computing stationary distributions is implemented in [QuantEcon.py](http://quantecon.org/quantecon-py). @@ -757,6 +764,11 @@ mc = qe.MarkovChain(P) mc.stationary_distributions # Show all stationary distributions ``` + + + + + ### Asymptotic stationarity Consider an everywhere positive stochastic matrix with unique stationary distribution $\psi^*$. @@ -767,17 +779,24 @@ For example, we have the following result (strict_stationary)= ```{prf:theorem} +:label: mc_gs_thm + If there exists an integer $m$ such that all entries of $P^m$ are -strictly positive, with unique stationary distribution $\psi^*$, then +strictly positive, then $$ \psi_0 P^t \to \psi^* \quad \text{ as } t \to \infty $$ + +where $\psi^*$ is the unique stationary distribution. ``` +This situation is often referred to as **asymptotic stationarity** or **global stability**. + +A proof of the theorem can be found in Chapter 4 of {cite}`sargent2023economic`, as well as many other sources. + -See, for example, {cite}`sargent2023economic` Chapter 4. @@ -848,103 +867,96 @@ Here You might like to try experimenting with different initial conditions. -#### An alternative illustration -We can show this in a slightly different way by focusing on the probability that $\psi_t$ puts on each state. -First, we write a function to draw initial distributions $\psi_0$ of size `num_distributions` - -```{code-cell} ipython3 -def generate_initial_values(num_distributions): - n = len(P) - - draws = np.random.randint(1, 10_000_000, size=(num_distributions,n)) - ψ_0s = draws/draws.sum(axis=1)[:, None] - - return ψ_0s -``` +#### Example: failure of convergence -We then write a function to plot the dynamics of $(\psi_0 P^t)(i)$ as $t$ gets large, for each state $i$ with different initial distributions -```{code-cell} ipython3 -def plot_distribution(P, ts_length, num_distributions): +Consider the periodic chain with stochastic matrix - # Get parameters of transition matrix - n = len(P) - mc = qe.MarkovChain(P) - ψ_star = mc.stationary_distributions[0] +$$ +P = +\begin{bmatrix} + 0 & 1 \\ + 1 & 0 \\ +\end{bmatrix} +$$ - ## Draw the plot - fig, axes = plt.subplots(nrows=1, ncols=n, figsize=[11, 5]) - plt.subplots_adjust(wspace=0.35) +This matrix does not satisfy the conditions of +{ref}`strict_stationary` because, as you can readily check, - ψ_0s = generate_initial_values(num_distributions) +* $P^m = P$ when $m$ is odd and +* $P^m = I$, the identity matrix, when $m$ is even. - # Get the path for each starting value - for ψ_0 in ψ_0s: - ψ_t = iterate_ψ(ψ_0, P, ts_length) +Hence there is no $m$ such that all elements of $P^m$ are strictly positive. - # Obtain and plot distributions at each state - for i in range(n): - axes[i].plot(range(0, ts_length), ψ_t[:,i], alpha=0.3) +Moreover, we can see that global stability does not hold. - # Add labels - for i in range(n): - axes[i].axhline(ψ_star[i], linestyle='dashed', lw=2, color = 'black', - label = fr'$\psi^*({i})$') - axes[i].set_xlabel('t') - axes[i].set_ylabel(fr'$\psi_t({i})$') - axes[i].legend() +For instance, if we start at $\psi_0 = (1,0)$, then $\psi_m = \psi_0 P^m$ is $(1, 0)$ when $m$ is even and $(0,1)$ when $m$ is odd. - plt.show() -``` +We can see similar phenomena in higher dimensions. -The following figure shows +The next figure illustrates this for a periodic Markov chain with three states. ```{code-cell} ipython3 -# Define the number of iterations -# and initial distributions -ts_length = 50 -num_distributions = 25 - -P = np.array([[0.971, 0.029, 0.000], - [0.145, 0.778, 0.077], - [0.000, 0.508, 0.492]]) - -plot_distribution(P, ts_length, num_distributions) -``` - -The convergence to $\psi^*$ holds for different initial distributions. +ψ_1 = (0.0, 0.0, 1.0) +ψ_2 = (0.5, 0.5, 0.0) +ψ_3 = (0.25, 0.25, 0.5) +ψ_4 = (1/3, 1/3, 1/3) +P = np.array([[0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + [1.0, 0.0, 0.0]]) +fig = plt.figure() +ax = fig.add_subplot(projection='3d') +colors = ['red','yellow', 'green', 'blue'] # Different colors for each initial point -#### Example: failure of convergence - +# Define the vertices of the unit simplex +v = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]]) -In the case of a periodic chain, with +# Define the faces of the unit simplex +faces = [ + [v[0], v[1], v[2]], + [v[0], v[1], v[3]], + [v[0], v[2], v[3]], + [v[1], v[2], v[3]] +] -$$ -P = -\begin{bmatrix} - 0 & 1 \\ - 1 & 0 \\ -\end{bmatrix} -$$ +def update(n): + ax.clear() + ax.set_xlim([0, 1]) + ax.set_ylim([0, 1]) + ax.set_zlim([0, 1]) + ax.view_init(45, 45) + + # Plot the 3D unit simplex as planes + simplex = Poly3DCollection(faces,alpha=0.05) + ax.add_collection3d(simplex) + + for idx, ψ_0 in enumerate([ψ_1, ψ_2, ψ_3, ψ_4]): + ψ_t = iterate_ψ(ψ_0, P, n+1) + + point = ψ_t[-1] + ax.scatter(point[0], point[1], point[2], color=colors[idx], s=60) + points = np.array(ψ_t) + ax.plot(points[:, 0], points[:, 1], points[:, 2], color=colors[idx],linewidth=0.75) + + return fig, -we find the distribution oscillates +anim = FuncAnimation(fig, update, frames=range(20), blit=False, repeat=False) +plt.close() +HTML(anim.to_jshtml()) +``` +This animation demonstrates the behavior of an irreducible and periodic stochastic matrix. -```{code-cell} ipython3 -P = np.array([[0, 1], - [1, 0]]) +The red, yellow, and green dots represent different initial probability distributions. -ts_length = 20 -num_distributions = 30 +The blue dot represents the unique stationary distribution. -plot_distribution(P, ts_length, num_distributions) -``` +Unlike Hamilton’s Markov chain, these initial distributions do not converge to the unique stationary distribution. -Indeed, this $P$ fails our asymptotic stationarity condition, since, as you can -verify, $P^t$ is not everywhere positive for any $t$. +Instead, they cycle periodically around the probability simplex, illustrating that asymptotic stability fails. (finite_mc_expec)= @@ -1105,42 +1117,6 @@ mc = qe.MarkovChain(P) ψ_star ``` -Solution 3: - -We find the distribution $\psi$ converges to the stationary distribution more quickly compared to the {ref}`hamilton's chain `. - -```{code-cell} ipython3 -ts_length = 10 -num_distributions = 25 -plot_distribution(P, ts_length, num_distributions) -``` - -In fact, the rate of convergence is governed by {ref}`eigenvalues` {cite}`sargent2023economic`. - -```{code-cell} ipython3 -P_eigenvals = np.linalg.eigvals(P) -P_eigenvals -``` - -```{code-cell} ipython3 -P_hamilton = np.array([[0.971, 0.029, 0.000], - [0.145, 0.778, 0.077], - [0.000, 0.508, 0.492]]) - -hamilton_eigenvals = np.linalg.eigvals(P_hamilton) -hamilton_eigenvals -``` - -More specifically, it is governed by the spectral gap, the difference between the largest and the second largest eigenvalue. - -```{code-cell} ipython3 -sp_gap_P = P_eigenvals[0] - np.diff(P_eigenvals)[0] -sp_gap_hamilton = hamilton_eigenvals[0] - np.diff(hamilton_eigenvals)[0] - -sp_gap_P > sp_gap_hamilton -``` - -We will come back to this when we discuss {ref}`spectral theory`. ```{solution-end} ```