Fix wealth dynamics lecture (#833)

* misc

* finished updates

Author: jstac

source/rst/wealth_dynamics.rst

@@ -221,15 +221,7 @@ normal in :math:`\mathbb R^3`.
 The value of :math:`c_r` should be close to zero, since rates of return
 on assets do not exhibit large trends.
 
-When we simulate a population of households, we will take
-
--  :math:`\{z_t\}` to be an **aggregate shock** that is common to all
-   households
--  :math:`\{\xi_t\}` and :math:`\{\zeta_t\}` to be **idiosyncratic
-   shocks**
-
-Idiosyncratic shocks are specific to individual households and
-independent across them.
+When we simulate a population of households, we will assume all shocks are idiosyncratic (i.e.,  specific to individual households and independent across them).
 
 Regarding the savings function :math:`s`, our default model will be
 
@@ -247,6 +239,7 @@ their wealth.
 We are using something akin to a fixed savings rate model, while
 acknowledging that low wealth households tend to save very little.
 
+
 Implementation
 ==============
 
@@ -267,6 +260,7 @@ Here's some type information to help Numba.
         ('b',      float64),    # aggregate shock parameter
         ('σ_z',    float64),    # aggregate shock parameter
         ('z_mean', float64),    # mean of z process
+        ('z_var', float64),     # variance of z process
         ('y_mean', float64),    # mean of y process
         ('R_mean', float64)     # mean of R process
     ]
@@ -280,15 +274,15 @@ the aggregate state and household wealth.
     class WealthDynamics:
         
         def __init__(self,
-                     w_hat=10.0,
+                     w_hat=1.0,
                      s_0=0.75, 
                      c_y=1.0, 
-                     μ_y=1.5,
+                     μ_y=1.0,
                      σ_y=0.2, 
                      c_r=0.05,
-                     μ_r=0.0, 
-                     σ_r=0.6, 
-                     a=0.85,  
+                     μ_r=0.1, 
+                     σ_r=0.5, 
+                     a=0.5,  
                      b=0.0, 
                      σ_z=0.1):
         
@@ -299,8 +293,8 @@ the aggregate state and household wealth.
             
             # Record stationary moments
             self.z_mean = b / (1 - a)
-            z_var = σ_z**2 / (1 - a**2)
-            exp_z_mean = np.exp(self.z_mean + z_var / 2)
+            self.z_var = σ_z**2 / (1 - a**2)
+            exp_z_mean = np.exp(self.z_mean + self.z_var / 2)
             self.R_mean = c_r * exp_z_mean + np.exp(μ_r + σ_r**2 / 2)
             self.y_mean = c_y * exp_z_mean + np.exp(μ_y + σ_y**2 / 2)
             
@@ -320,82 +314,54 @@ the aggregate state and household wealth.
                           self.a, self.b, self.σ_z)
             return parameters
             
-        def update_z(self, z):
+        def update_states(self, w, z):
             """
-            Update exogenous state one period, given current state z.
+            Update one period, given current wealth w and persistent 
+            state z.
             """   
+
             # Simplify names
             params = self.parameters()
             w_hat, s_0, c_y, μ_y, σ_y, c_r, μ_r, σ_r, a, b, σ_z = params
-            
-            # Update z
             zp = a * z + b + σ_z * np.random.randn()
-            return zp
-            
-        def update_wealth(self, w, zp):
-            """
-            Update one period, given current wealth w and next period 
-            exogenous state zp.
-            """   
-            # Simplify names
-            params = self.parameters()
-            w_hat, s_0, c_y, μ_y, σ_y, c_r, μ_r, σ_r, a, b, σ_z = params
 
             # Update wealth
             y = c_y * np.exp(zp) + np.exp(μ_y + σ_y * np.random.randn())   
             wp = y
             if w >= w_hat:
                 R = c_r * np.exp(zp) + np.exp(μ_r + σ_r * np.random.randn())
                 wp += R * s_0 * w 
-            return wp
+            return wp, zp
             
-        
-Here’s a function to simulate the aggregate state:
+
+Here's function to simulate the time series of wealth for in individual households.
 
 .. code:: ipython3
 
     @njit
-    def z_time_series(wdy, ts_length=10_000):
+    def wealth_time_series(wdy, w_0, n):
         """
-        Generate a single time series of length ts_length for the 
-        aggregate shock process {z_t}.
+        Generate a single time series of length n for wealth given 
+        initial value w_0.
 
-        * wdy is an instance of WealthDynamics
-
-        """
-        z = np.empty(ts_length)
-        z[0] = wdy.z_mean
-        for t in range(ts_length-1):
-            z[t+1] = wdy.update_z(z[t])
-        return z
+        The initial persistent state z_0 for each household is drawn from
+        the stationary distribution of the AR(1) process.
 
-Here’s a function to simulate the wealth of a single household, taking the
-time path of the aggregate state as an argument.
+            * wdy: an instance of WealthDynamics
+            * w_0: scalar
+            * n: int
 
-.. code:: ipython3
-
-    @njit
-    def wealth_time_series(wdy, w_0, z):
-        """
-        Generate a single time series for wealth given aggregate shock 
-        sequence z.
 
-        * wdy: an instance of WealthDynamics
-        * w_0: a scalar representing initial wealth.
-        * z: array_like
-        
-        The returned time series w satisfies len(w) = len(z).
         """
-        n = len(z)
+        z = wdy.z_mean + np.sqrt(wdy.z_var) * np.random.randn()
         w = np.empty(n)
         w[0] = w_0
         for t in range(n-1):
-            w[t+1] = wdy.update_wealth(w[t], z[t+1])
+            w[t+1], z = wdy.update_states(w[t], z)
         return w
 
 
-Finally, here's function to simulate a cross section of households forward in
-time.
+Now here's function to simulate a cross section of households forward in time.
 
 Note the use of parallelization to speed up computation.
 
@@ -414,16 +380,16 @@ Note the use of parallelization to speed up computation.
         j = shift_length.  
         
         Returns the new distribution.
+
         """
         new_distribution = np.empty_like(w_distribution)
     
-        z = z_time_series(wdy, ts_length=shift_length)
-    
         # Update each household
         for i in prange(len(new_distribution)):
+            z = wdy.z_mean + np.sqrt(wdy.z_var) * np.random.randn()
             w = w_distribution[i]
             for t in range(shift_length-1):
-                w = wdy.update_wealth(w, z[t+1])
+                w, z = wdy.update_states(w, z)
             new_distribution[i] = w
         return new_distribution
     
@@ -450,8 +416,8 @@ Let's look at the wealth dynamics of an individual household.
 
     wdy = WealthDynamics()
     
-    z = z_time_series(wdy, ts_length=10_000)
-    w = wealth_time_series(wdy, wdy.y_mean, z)
+    ts_length = 200
+    w = wealth_time_series(wdy, wdy.y_mean, ts_length)
     
     fig, ax = plt.subplots()
     ax.plot(w)
@@ -474,7 +440,7 @@ curve and Gini coefficient.
 
 .. code:: ipython3
 
-    def generate_lorenz_and_gini(wdy, num_households=250_000, T=500):
+    def generate_lorenz_and_gini(wdy, num_households=100_000, T=500):
         """
         Generate the Lorenz curve data and gini coefficient corresponding to a 
         WealthDynamics mode by simulating num_households forward to time T.
@@ -485,8 +451,7 @@ curve and Gini coefficient.
         ψ_star = update_cross_section(wdy, ψ_0, shift_length=T)
         return qe.gini_coefficient(ψ_star), qe.lorenz_curve(ψ_star)
 
-Now we investigate how the Lorenz curves associated with the wealth distribution change
-as return to savings increase.
+Now we investigate how the Lorenz curves associated with the wealth distribution change as return to savings varies.
 
 The code below plots Lorenz curves for three different values of :math:`\mu_r`. 
 
@@ -499,7 +464,7 @@ In fact the code, which is JIT compiled and parallelized, runs extremely fast re
 .. code:: ipython3
 
     fig, ax = plt.subplots()
-    μ_r_vals = np.linspace(0.0, 0.05, 3)
+    μ_r_vals = (0.0, 0.025, 0.05)
     gini_vals = []
     
     for μ_r in μ_r_vals:
@@ -544,30 +509,49 @@ Now let's check the Gini coefficient.
 Once again, we see that inequality increases as returns on financial income
 rise.
 
-This makes sense because higher returns reinforce the effect of existing wealth.
+Let's finish this section by investigating what happens when we change the
+volatility term :math:`\sigma_r` in financial returns.
 
 
+.. code:: ipython3
+
+    fig, ax = plt.subplots()
+    σ_r_vals = (0.35, 0.45, 0.52)
+    gini_vals = []
+    
+    for σ_r in σ_r_vals:
+        wdy = WealthDynamics(σ_r=σ_r)
+        gv, (f_vals, l_vals) = generate_lorenz_and_gini(wdy)
+        ax.plot(f_vals, l_vals, label=f'$\psi^*$ at $\mu_r = {σ_r:0.2}$')
+        gini_vals.append(gv)
+        
+    ax.plot(f_vals, f_vals, label='equality')
+    ax.legend(loc="upper left")
+    plt.show()
+
+
+We see that greater volatility has the effect of increasing inequality in this model.
+
 
 Exercises
 =========
 
 Exercise 1
 ----------
 
-For a wealth or income distribution with Pareto tail, a higher tail index suggests lower
-In fact, it is possible to prove that the Gini coefficient of the Pareto
+For a wealth or income distribution with Pareto tail, a higher tail index suggests lower inequality.
+
+Indeed, it is possible to prove that the Gini coefficient of the Pareto
 distribution with tail index :math:`a` is :math:`1/(2a - 1)`.
 
 To the extent that you can, confirm this by simulation.
 
 In particular, generate a plot of the Gini coefficient against the tail index
-using both the theoretical value and the value computed from a sample via
-``qe.gini_coefficient``.
+using both the theoretical value just given and the value computed from a sample via ``qe.gini_coefficient``.
 
 For the values of the tail index, use ``a_vals = np.linspace(1, 10, 25)``.
 
-Use sample of size 1,000 for each :math:`a` and the sampling method for generating
-Pareto draws employed in the discussion of Lorenz curves for the Pareto distribution.
+Use sample of size 1,000 for each :math:`a` and the sampling method for generating Pareto draws employed in the discussion of Lorenz curves for the Pareto distribution.
 
 To the extend that you can, interpret the monotone relationship between the
 Gini index and :math:`a`.