Merge branch 'main' into blume_easley

Author: HumphreyYang

.github/workflows/collab.yml

@@ -11,10 +11,11 @@ jobs:
         with:
           ref: ${{ github.event.pull_request.head.sha }}
       # Install build software
-      - name: Install Build Software & LaTeX (kalman_2)
+      - name: Install Build Software & LaTeX
         shell: bash -l {0}
         run: |
           pip install jupyter-book==1.0.3 quantecon-book-theme==0.8.2 sphinx-tojupyter==0.3.0 sphinxext-rediraffe==0.2.7 sphinxcontrib-youtube==1.3.0 sphinx-togglebutton==0.3.2 arviz sphinx-proof sphinx-exercise sphinx-reredirects
+          apt-get update
           apt-get install dvipng texlive texlive-latex-extra texlive-fonts-recommended cm-super    
       - name: Check nvidia drivers
         shell: bash -l {0}

lectures/aiyagari.md

@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -26,10 +28,9 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython
----
-tags: [hide-output]
----
+```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install quantecon
 ```
 
@@ -54,7 +55,7 @@ The Aiyagari model has been used to investigate many topics, including
 
 Let's start with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
 plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 import numpy as np
@@ -208,7 +209,7 @@ when the parameters change.
 
 The class also includes a default set of parameters that we'll adopt unless otherwise specified.
 
-```{code-cell} python3
+```{code-cell} ipython3
 class Household:
     """
     This class takes the parameters that define a household asset accumulation
@@ -318,7 +319,7 @@ def asset_marginal(s_probs, a_size, z_size):
 
 As a first example of what we can do, let's compute and plot an optimal accumulation policy at fixed prices.
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Example prices
 r = 0.03
 w = 0.956
@@ -366,7 +367,7 @@ The following code draws aggregate supply and demand curves.
 
 The intersection gives equilibrium interest rates and capital.
 
-```{code-cell} python3
+```{code-cell} ipython3
 A = 1.0
 N = 1.0
 α = 0.33
@@ -410,7 +411,7 @@ def prices_to_capital_stock(am, r):
     # Extract the marginal distribution for assets
     asset_probs = asset_marginal(stationary_probs, am.a_size, am.z_size)
     # Return K
-    return np.sum(asset_probs * am.a_vals)
+    return asset_probs @ am.a_vals
 
 
 # Create an instance of Household

lectures/cass_koopmans_1.md

@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.16.4
+    jupytext_version: 1.17.2
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -62,6 +62,13 @@ The lecture uses important ideas including
   long but finite-horizon economies.
 - A **stable manifold** and a **phase plane**
 
+In addition to what's in Anaconda, this lecture will need the following libraries:
+
+```{code-cell} ipython
+:tags: [hide-output]
+!pip install quantecon
+```
+
 Let's start with some standard imports:
 
 ```{code-cell} ipython3
@@ -614,7 +621,7 @@ tolerance bounds), we stop.
 def bisection(pp, c0, k0, T=10, tol=1e-4, max_iter=500, k_ter=0, verbose=True):
 
     # initial boundaries for guess c0
-    c0_upper = pp.f(k0)
+    c0_upper = pp.f(k0) + (1 - pp.δ) * k0  
     c0_lower = 0
 
     i = 0
@@ -648,7 +655,7 @@ def plot_paths(pp, c0, k0, T_arr, k_ter=0, k_ss=None, axs=None):
 
     if axs is None:
         fix, axs = plt.subplots(1, 3, figsize=(16, 4))
-    ylabels = ['$c_t$', '$k_t$', '$\mu_t$']
+    ylabels = ['$c_t$', '$k_t$', r'$\mu_t$']
     titles = ['Consumption', 'Capital', 'Lagrange Multiplier']
 
     c_paths = []
@@ -801,6 +808,76 @@ A rule of thumb for the planner is
 The planner accomplishes this by adjusting the saving rate $\frac{f(K_t) - C_t}{f(K_t)}$
 over time.
 
+```{exercise}
+:label: ck1_ex1
+
+The turnpike property is independent of the initial condition 
+$K_0$ provided that $T$ is sufficiently large.
+
+Expand the `plot_paths` function so that it plots trajectories for multiple initial points using `k0s = [k_ss*2, k_ss*3, k_ss/3]`.
+```
+
+```{solution-start} ck1_ex1
+:class: dropdown
+```
+
+Here is one solution
+
+```{code-cell} ipython3
+def plot_multiple_paths(pp, c0, k0s, T_arr, k_ter=0, k_ss=None, axs=None):
+    if axs is None:
+        fig, axs = plt.subplots(1, 3, figsize=(16, 4))
+        
+    ylabels = ['$c_t$', '$k_t$', r'$\mu_t$']
+    titles = ['Consumption', 'Capital', 'Lagrange Multiplier']
+
+    colors = plt.cm.viridis(np.linspace(0, 1, len(k0s)))
+    
+    all_c_paths = []
+    all_k_paths = []
+    
+    for i, k0 in enumerate(k0s):
+        k0_c_paths = []
+        k0_k_paths = []
+        
+        for T in T_arr:
+            c_vec, k_vec = bisection(pp, c0, k0, T, k_ter=k_ter, verbose=False)
+            k0_c_paths.append(c_vec)
+            k0_k_paths.append(k_vec)
+
+            μ_vec = pp.u_prime(c_vec)
+            paths = [c_vec, k_vec, μ_vec]
+
+            for j in range(3):
+                axs[j].plot(paths[j], color=colors[i], 
+                           label=f'$k_0 = {k0:.2f}$' if j == 0 and T == T_arr[0] else "", alpha=0.7)
+                axs[j].set(xlabel='t', ylabel=ylabels[j], title=titles[j])
+
+            if k_ss is not None and i == 0 and T == T_arr[0]:
+                axs[1].axhline(k_ss, c='k', ls='--', lw=1)
+
+            axs[1].axvline(T+1, c='k', ls='--', lw=1)
+            axs[1].scatter(T+1, paths[1][-1], s=80, color=colors[i])
+        
+        all_c_paths.append(k0_c_paths)
+        all_k_paths.append(k0_k_paths)
+    
+    # Add legend if multiple initial points
+    if len(k0s) > 1:
+        axs[0].legend()
+
+    return all_c_paths, all_k_paths
+```
+
+```{code-cell} ipython3
+_ = plot_multiple_paths(pp, 0.3, [k_ss*2, k_ss*3, k_ss/3], [250, 150, 75, 50], k_ss=k_ss)
+```
+
+We see that the turnpike property holds for various initial values of $K_0$.
+
+```{solution-end}
+```
+
 Let's calculate and  plot the saving rate.
 
 ```{code-cell} ipython3
@@ -1075,15 +1152,15 @@ studied in {doc}`Cass-Koopmans Competitive Equilibrium <cass_koopmans_2>` is a f
 ### Exercise
 
 ```{exercise}
-:label: ck1_ex1
+:label: ck1_ex2
 
 - Plot the optimal consumption, capital, and saving paths when the
   initial capital level begins at 1.5 times the steady state level
   as we shoot towards the steady state at $T=130$.
 - Why does the saving rate respond as it does?
 ```
 
-```{solution-start} ck1_ex1
+```{solution-start} ck1_ex2
 :class: dropdown
 ```
 

lectures/egm_policy_iter.md

@@ -35,7 +35,7 @@ We found time iteration to be significantly more accurate and efficient.
 
 In this lecture, we'll look at a clever twist on time iteration called the **endogenous grid method** (EGM).
 
-EGM is a numerical method for implementing policy iteration invented by [Chris Carroll](http://www.econ2.jhu.edu/people/ccarroll/).
+EGM is a numerical method for implementing policy iteration invented by [Chris Carroll](https://econ.jhu.edu/directory/christopher-carroll/).
 
 The original reference is {cite}`Carroll2006`.
 

lectures/ge_arrow.md

@@ -964,15 +964,15 @@ class RecurCompetitive:
     def price_risk_free_bond(self):
         "Give Q, compute price of one-period risk free bond"
 
-        PRF = np.sum(self.Q, 0)
+        PRF = np.sum(self.Q, axis=1)
         self.PRF = PRF
 
         return PRF
 
     def risk_free_rate(self):
         "Given Q, compute one-period gross risk-free interest rate R"
 
-        R = np.sum(self.Q, 0)
+        R = np.sum(self.Q, axis=1)
         R = np.reciprocal(R)
         self.R = R
 

lectures/kalman_2.md

@@ -447,9 +447,9 @@ def simulate_workers(worker, T, ax, mu_0=None, Sigma_0=None,
     A, C, G, R = worker.A, worker.C, worker.G, worker.R
     xhat_0, Σ_0 = worker.xhat_0, worker.Σ_0
     
-    if isinstance(mu_0, type(None)):
+    if mu_0 is None:
         mu_0 = xhat_0
-    if isinstance(Sigma_0, type(None)):
+    if Sigma_0 is None:
         Sigma_0 = worker.Σ_0
         
     ss = LinearStateSpace(A, C, G, np.sqrt(R), 
@@ -471,12 +471,12 @@ def simulate_workers(worker, T, ax, mu_0=None, Sigma_0=None,
         kalman.update(y[i])
         x_hat, Σ = kalman.x_hat, kalman.Sigma
         Σ_t.append(Σ)
-        [y_hat_t[i]] = worker.G @ x_hat
-        [u_hat_t[i]] = x_hat[1]
+        y_hat_t[i] = (worker.G @ x_hat).item()
+        u_hat_t[i] = x_hat[1].item()
 
-    if diff == True:
+    if diff :
         title = ('Difference between inferred and true work ethic over time' 
-                 if title == None else title)
+                 if title is None else title)
         
         ax.plot(u_hat_t - u_0, alpha=.5)
         ax.axhline(y=0, color='grey', linestyle='dashed')
@@ -485,10 +485,10 @@ def simulate_workers(worker, T, ax, mu_0=None, Sigma_0=None,
         ax.set_title(title)
         
     else:
-        label_line = (r'$E[u_t|y^{t-1}]$' if name == None 
+        label_line = (r'$E[u_t|y^{t-1}]$' if name is None 
                       else name)
         title = ('Inferred work ethic over time' 
-                if title == None else title)
+                if title is None else title)
         
         u_hat_plot = ax.plot(u_hat_t, label=label_line)
         ax.axhline(y=u_0, color=u_hat_plot[0].get_color(), 
@@ -525,9 +525,7 @@ for i, (uhat_0, α, β) in enumerate(zip(uhat_0s, αs, βs)):
     simulate_workers(worker, T, ax,
                     # By setting diff=False, it will give u_t
                     diff=False, name=r'$u_{{{}, t}}$'.format(i))
-    
-ax.axhline(y=u_0, xmin=0, xmax=0, color='grey', 
-           linestyle='dashed', label=r'$u_{i, 0}$')
+
 ax.legend(bbox_to_anchor=(1, 0.5))
 plt.show()
 ```
@@ -554,8 +552,6 @@ for i, (uhat_0, α, β) in enumerate(zip(uhat_0s, αs, βs)):
     
 # This controls the boundary of plots
 ax.set_ylim(ymin=-3, ymax=3)
-ax.axhline(y=u_0, xmin=0, xmax=0, color='grey', 
-           linestyle='dashed', label=r'$u_{i, 0}$')
 ax.legend(bbox_to_anchor=(1, 0.5))
 plt.show()
 ```

lectures/mccall_model.md

@@ -404,7 +404,7 @@ class McCallModel:
         # Evaluate value for each state-action pair
         # Consider action = accept or reject the current offer
         accept = w[i] / (1 - β)
-        reject = c + β * np.sum(v * q)
+        reject = c + β * (v @ q)
 
         return np.array([accept, reject])
 ```
@@ -488,7 +488,7 @@ def compute_reservation_wage(mcm,
 
     # == Now compute the reservation wage == #
 
-    return (1 - β) * (c + β * np.sum(v * q))
+    return (1 - β) * (c + β * (v @ q))
 ```
 
 The next line computes the reservation wage at  default parameters
@@ -622,13 +622,13 @@ def compute_reservation_wage_two(mcm,
 
     # == First compute h == #
 
-    h = np.sum(w * q) / (1 - β)
+    h = (w @ q) / (1 - β)
     i = 0
     error = tol + 1
     while i < max_iter and error > tol:
 
         s = np.maximum(w / (1 - β), h)
-        h_next = c + β * np.sum(s * q)
+        h_next = c + β * (s @ q)
 
         error = np.abs(h_next - h)
         i += 1

lectures/mccall_q.md

@@ -171,7 +171,7 @@ class McCallModel:
         # Evaluate value for each state-action pair
         # Consider action = accept or reject the current offer
         accept = w[i] / (1 - β)
-        reject = c + β * np.sum(v * q)
+        reject = c + β * (v @ q)
 
         return np.array([accept, reject])