MAINT: transfer np.sum(a * b) to a @ b

Author: mmcky

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,9 @@ The following code draws aggregate supply and demand curves.
 
 The intersection gives equilibrium interest rates and capital.
 
-```{code-cell} python3
+```{code-cell} ipython3
+:tags: hide-input
+
 A = 1.0
 N = 1.0
 α = 0.33
@@ -439,3 +442,85 @@ ax.legend(loc='upper right')
 
 plt.show()
 ```
+
+```{code-cell} ipython3
+k_vals_orig = k_vals
+```
+
+```{code-cell} ipython3
+A = 1.0
+N = 1.0
+α = 0.33
+β = 0.96
+δ = 0.05
+
+
+def r_to_w(r):
+    """
+    Equilibrium wages associated with a given interest rate r.
+    """
+    return A * (1 - α) * (A * α / (r + δ))**(α / (1 - α))
+
+def rd(K):
+    """
+    Inverse demand curve for capital.  The interest rate associated with a
+    given demand for capital K.
+    """
+    return A * α * (N / K)**(1 - α) - δ
+
+
+def prices_to_capital_stock(am, r):
+    """
+    Map prices to the induced level of capital stock.
+
+    Parameters:
+    ----------
+
+    am : Household
+        An instance of an aiyagari_household.Household
+    r : float
+        The interest rate
+    """
+    w = r_to_w(r)
+    am.set_prices(r, w)
+    aiyagari_ddp = DiscreteDP(am.R, am.Q, β)
+    # Compute the optimal policy
+    results = aiyagari_ddp.solve(method='policy_iteration')
+    # Compute the stationary distribution
+    stationary_probs = results.mc.stationary_distributions[0]
+    # Extract the marginal distribution for assets
+    asset_probs = asset_marginal(stationary_probs, am.a_size, am.z_size)
+    # Return K
+    return asset_probs @ am.a_vals
+
+
+# Create an instance of Household
+am = Household(a_max=20)
+
+# Use the instance to build a discrete dynamic program
+am_ddp = DiscreteDP(am.R, am.Q, am.β)
+
+# Create a grid of r values at which to compute demand and supply of capital
+num_points = 20
+r_vals = np.linspace(0.005, 0.04, num_points)
+
+# Compute supply of capital
+k_vals = np.empty(num_points)
+for i, r in enumerate(r_vals):
+    k_vals[i] = prices_to_capital_stock(am, r)
+
+# Plot against demand for capital by firms
+fig, ax = plt.subplots(figsize=(11, 8))
+ax.plot(k_vals, r_vals, lw=2, alpha=0.6, label='supply of capital')
+ax.plot(k_vals, rd(k_vals), lw=2, alpha=0.6, label='demand for capital')
+ax.grid()
+ax.set_xlabel('capital')
+ax.set_ylabel('interest rate')
+ax.legend(loc='upper right')
+
+plt.show()
+```
+
+```{code-cell} ipython3
+np.allclose(k_vals_orig, k_vals)
+```