update eigenvalue in mc and add example to mc2

Author: HumphreyYang

lectures/markov_chains_I.md

@@ -13,7 +13,7 @@ kernelspec:
 
 +++ {"user_expressions": []}
 
-# Markov Chains I
+# Markov Chains: Basic Concepts and Stationarity
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
@@ -270,12 +270,12 @@ $$
 
 ```{code-cell} ipython3
 nodes = ['DG', 'DC', 'NG', 'NC', 'AG', 'AC']
-trans_matrix = [[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],
-                [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],
-                [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],
-                [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],
-                [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],
-                [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]]
+P = [[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],
+     [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],
+     [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],
+     [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],
+     [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],
+     [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]]
 ```
 
 ```{code-cell} ipython3
@@ -285,7 +285,7 @@ label_dict = {}
 
 for start_idx, node_start in enumerate(nodes):
     for end_idx, node_end in enumerate(nodes):
-        value = trans_matrix[start_idx][end_idx]
+        value = P[start_idx][end_idx]
         if value != 0:
             G.add_edge(node_start,node_end, weight=value, len=100)
             
@@ -815,6 +815,7 @@ See, for example, {cite}`sargent2023economic` Chapter 4.
 
 +++ {"user_expressions": []}
 
+(hamilton)=
 #### Example: Hamilton's Chain
 
 Hamilton's chain satisfies the conditions of the theorem because $P^2$ is everywhere positive:
@@ -823,7 +824,6 @@ Hamilton's chain satisfies the conditions of the theorem because $P^2$ is everyw
 P = np.array([[0.971, 0.029, 0.000],
               [0.145, 0.778, 0.077],
               [0.000, 0.508, 0.492]])
-P = np.array(P)
 P @ P
 ```
 
@@ -1099,8 +1099,7 @@ Compare your solution to the paper.
 :class: dropdown
 ```
 
-1. 
-
+1.
 
 ```{code-cell} ipython3
 :tags: [hide-output]
@@ -1155,12 +1154,43 @@ mc = qe.MarkovChain(P)
 
 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 in 
+
+TODO: add a reference to eigen II
+
 ```{solution-end}
 ```
 
@@ -1192,7 +1222,7 @@ In this exercise,
 
 1. 
 
-Although $P$ is not every positive, $P^m$ when $m=3$ is everywhere positive. 
+Although $P$ is not every positive, $P^m$ when $m=3$ is everywhere positive.
 
 ```{code-cell} ipython3
 P = np.array([[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],
@@ -1209,18 +1239,12 @@ So it satisfies the requirement.
 
 2.
 
-We can see the distribution $\psi$ converges to the stationary distribution quickly regardless of the initial distributions
+We find the distribution $\psi$ converges to the stationary distribution quickly regardless of the initial distributions
 
 ```{code-cell} ipython3
 ts_length = 30
 num_distributions = 20
 nodes = ['DG', 'DC', 'NG', 'NC', 'AG', 'AC']
-P = [[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],
-     [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],
-     [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],
-     [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],
-     [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],
-     [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]]
 
 # Get parameters of transition matrix
 n = len(P)
@@ -1277,4 +1301,4 @@ Therefore $P^{k+1} \mathbf 1 = P P^k \mathbf 1 = P \mathbf 1 = \mathbf 1$
 The proof is done.
 
 ```{solution-end}
-```
+```

lectures/markov_chains_II.md

@@ -13,7 +13,7 @@ kernelspec:
 
 +++ {"user_expressions": []}
 
-# Markov Chains II
+# Markov Chains: Irreducibility and Ergodicity
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
@@ -262,7 +262,7 @@ This is one aspect of the concept  of ergodicity.
 ### Example 2
 
 
-Another example is Hamilton {cite}`Hamilton2005` dynamics {ref}`discussed above <mc_eg2>`.
+Another example is Hamilton {cite}`Hamilton2005` dynamics {ref}`discussed before <mc_eg2>`.
 
 The diagram of the Markov chain shows that it is **irreducible**.
 
@@ -298,6 +298,61 @@ plt.show()
 
 ### Example 3
 
+Let's look at one more example with six states {ref}`discussed before <mc_eg3>`.
+
+
+$$
+$$
+P :=
+\left(
+  \begin{array}{cccccc}
+0.86 & 0.11 & 0.03 & 0.00 & 0.00 & 0.00 \\
+0.52 & 0.33 & 0.13 & 0.02 & 0.00 & 0.00 \\
+0.12 & 0.03 & 0.70 & 0.11 & 0.03 & 0.01 \\
+0.13 & 0.02 & 0.35 & 0.36 & 0.10 & 0.04 \\
+0.00 & 0.00 & 0.09 & 0.11 & 0.55 & 0.25 \\
+0.00 & 0.00 & 0.09 & 0.15 & 0.26 & 0.50
+  \end{array}
+\right)
+$$
+$$
+
+
+The graph for the chain shows states are densely connected indicating that it is **irreducible**.
+
+Similar to previous examples, the sample path averages for each state converges to the stationary distribution
+
+```{code-cell} ipython3
+P = [[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],
+     [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],
+     [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],
+     [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],
+     [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],
+     [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]]
+
+ts_length = 10_000
+mc = qe.MarkovChain(P)
+ψ_star = mc.stationary_distributions[0]
+fig, ax = plt.subplots(figsize=(9, 6))
+X = mc.simulate(ts_length)
+# Center the plot at 0
+ax.set_ylim(-0.25, 0.25)
+ax.axhline(0, linestyle='dashed', lw=2, color = 'black', alpha=0.4)
+
+
+for x0 in range(6):
+    # Calculate the fraction of time for each state
+    X_bar = (X == x0).cumsum() / (1 + np.arange(ts_length, dtype=float))
+    ax.plot(X_bar - ψ_star[x0], label=f'$X = {x0+1} $')
+    ax.set_xlabel('t')
+    ax.set_ylabel(r'fraction of time spent in a state $- \psi^* (x)$')
+
+ax.legend()
+plt.show()
+```
+
+### Example 4
+
 Let's look at another example with two states: 0 and 1.
 
 
@@ -444,16 +499,14 @@ In this exercise,
 Use the technique we learnt before, we can take the power of the transition matrix
 
 ```{code-cell} ipython3
-P = [
-    [0.222, 0.222, 0.215, 0.187, 0.081, 0.038, 0.029, 0.006],
-    [0.221, 0.22,  0.215, 0.188, 0.082, 0.039, 0.029, 0.006],
-    [0.207, 0.209, 0.21,  0.194, 0.09,  0.046, 0.036, 0.008],
-    [0.198, 0.201, 0.207, 0.198, 0.095, 0.052, 0.04,  0.009],
-    [0.175, 0.178, 0.197, 0.207, 0.11,  0.067, 0.054, 0.012],
-    [0.182, 0.184, 0.2,   0.205, 0.106, 0.062, 0.05,  0.011],
-    [0.123, 0.125, 0.166, 0.216, 0.141, 0.114, 0.094, 0.021],
-    [0.084, 0.084, 0.142, 0.228, 0.17,  0.143, 0.121, 0.028]
-    ]
+P = [[0.222, 0.222, 0.215, 0.187, 0.081, 0.038, 0.029, 0.006],
+     [0.221, 0.22,  0.215, 0.188, 0.082, 0.039, 0.029, 0.006],
+     [0.207, 0.209, 0.21,  0.194, 0.09,  0.046, 0.036, 0.008],
+     [0.198, 0.201, 0.207, 0.198, 0.095, 0.052, 0.04,  0.009],
+     [0.175, 0.178, 0.197, 0.207, 0.11,  0.067, 0.054, 0.012],
+     [0.182, 0.184, 0.2,   0.205, 0.106, 0.062, 0.05,  0.011],
+     [0.123, 0.125, 0.166, 0.216, 0.141, 0.114, 0.094, 0.021],
+     [0.084, 0.084, 0.142, 0.228, 0.17,  0.143, 0.121, 0.028]]
 
 P = np.array(P)
 codes_B = ('1','2','3','4','5','6','7','8')
@@ -476,11 +529,9 @@ ts_length = 1000
 mc = qe.MarkovChain(P)
 fig, ax = plt.subplots(figsize=(9, 6))
 X = mc.simulate(ts_length)
-# Center the plot at 0
 ax.set_ylim(-0.25, 0.25)
 ax.axhline(0, linestyle='dashed', lw=2, color = 'black', alpha=0.4)
 
-
 for x0 in range(8):
     # Calculate the fraction of time for each worker
     X_bar = (X == x0).cumsum() / (1 + np.arange(ts_length, dtype=float))