[markov_chain_I] latest edits (#510)

* latest edits to mc I lecture

* misc

* misc

* update the label

there are two labels for this theorem.

* Revert "update the label"

This reverts commit be885de.

* update label

update the ref to the theorem

---------

Co-authored-by: Longye Tian <longye.tian@outlook.com>

Author: jstac

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 <hamilton>`.
-
-```{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<eigen>` {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<spec_markov>`.
 
 ```{solution-end}
 ```