fix

Author: Smit-create

lectures/olg.md

@@ -728,44 +728,38 @@ def simulate_ts(m, x0_values, ts_length):
 simulate_ts(m_crra, x0, ts_length)
 ```
 
-## Exercises
 
-TODO Add exercise-solution environment
 
+## Exercises
 
 
-### Exercise 1
+```{exercise}
+:label: olg_ex1
 
-To solve the equilibrium of the CRRA case numerically using [](equilibrium_crra).
+Solve the equilibrium of the CRRA case numerically using [](equilibrium_crra).
 
 Plot the equilibrium quantity and equilibrium price in the equilibrium plot with the CRRA utility (TODO label and refer the plot generated by the code).
+```
 
 
-
-### Solution to exercise 1
-
-To solve the equation we need to turn to Newton's method.
-
-
-
-```{math}
-    L (1-\alpha)(K_t / L)^{\alpha} \left [ 1 + \beta^{-1/\gamma} R_{t+1}^{(\gamma-1)/\gamma} \right ]^{-1} = L\left (\frac{R_{t+1}}{\alpha} \right )^{1/(\alpha - 1)}
+```{solution-start} olg_ex1
+:class: dropdown
 ```
 
+To solve the equation we need to turn to Newton's method. `find_Rstar` is used to find $R^*_{t+1}$ by finding the zero of equation [](equilibrium_crra) using the helper function `find_Rstar_newton` for a given value of $K_t$. Similary, `find_Kstar` finds the equilibrium quantity $K^*_{t+1}$ using the value of $R^*_{t+1}$.
+
 ```{code-cell} ipython3
 def find_Rstar_newton(x, K_prev, model):
     α, β, γ, L = model.α, model.β, model.u.__defaults__[0], model.L
     lhs = L * (1-α) * (K_prev / L)**α
     lhs /= (1 + β**(-1/γ) * x**((γ-1)/γ))
-    # To should always non-negative for a real solution.
-    assert x>=0, x
     rhs = L * (x / α)**(1/(α-1))
     return lhs - rhs
 ```
 
 ```{code-cell} ipython3
 def find_Rstar(K_prev, model):
-    return optimize.newton(find_Rstar_newton, K_prev*100, args=(K_prev, model))
+    return optimize.newton(find_Rstar_newton, 0.5, args=(K_prev, model))
 
 def find_Kstar(R_star, model):
     return model.L * (R_star / model.α)**(1/(model.α-1))
@@ -792,25 +786,27 @@ def plot_ks_rs(K_t_vals, model):
 ```
 
 ```{code-cell} ipython3
-K_t_vals = np.linspace(500, 1000, 50)
+K_t_vals = np.linspace(0.1, 50, 50)
 m_crra = create_olg_model(u=crra)
 plot_ks_rs(K_t_vals, m_crra)
 ```
 
+```{solution-end}
+```
 
-
-### Exercise 2
-
+```{exercise}
+:label: olg_ex2
 
 Let's keep the model the same except for replacing the utility function $u(c)$ in {eq}`eq_crra`  with a quasilinear form $u(c)=c + c^{\theta}$.
 
-Like what we did in the crra case we don't have an analytical solution.
+Like what we did in the CRRA case we don't have an analytical solution.
 
 Try to compute the time path capital $\{k_t\}$ in this case.
+```
 
-+++
-
-### Solution to exercise 2
+```{solution-start} olg_ex2
+:class: dropdown
+```
 
 To get the time path capital $\{k_t\}$ first we need to solve the household's utility maximization problem for the optimal consumption and optimal saving.
 
@@ -825,74 +821,62 @@ Let $k_t := K_t / L$.
 
 Since [](aggregate_supply_capital_log_olg), [](wage_2) and [](interest_rate_2) the Euler equation becomes
 ```{math}
-:label:
+:label: euler_quasilinear1
     1 + \theta ((1-\alpha)k^{\alpha}_t - k_{t+1})^{\theta-1} = \beta \alpha k^{\alpha - 1}_t + \beta (\alpha k^{\alpha - 1}_t)^{\theta} \theta k_{t+1}^{\theta - 1}
 ```
 
 Obviously $k_{t+1}$ cannot be solved by pencil and paper.
 
 To solve for $k_{t+1}$ we need to turn to the newton's method.
 
+Let's start by defining the utility function.
+
 ```{code-cell} ipython3
-# Hi Smit please fill in the trial answers here
+def u_quasilinear(c, θ=4):
+    return c + c**θ
 ```
 
 ```{code-cell} ipython3
-
+def solve_for_k_t1(x, k_t, model):
+    α, β, L, θ = model.α, model.β, model.L, model.u.__defaults__[0]
+    l = 1 + θ * ((1 - α) * k_t**α - x)**(θ - 1)
+    r = β * α * k_t**(α - 1)
+    r += β * (α * k_t**(α - 1))**θ * θ * x**(θ - 1)
+    return l - r
 ```
 
 ```{code-cell} ipython3
-# Smit's previous code
-
-# def u_prime(c, θ):
-#     return 1 + θ * c ** (θ - 1)
-
-# def solve_for_c(x, θ, model, k_t, k_t_1):
-#     l = u_prime(x, θ)
-#     R = model.α * k_t_1**(model.α-1)
-#     r = model.β * R
-#     y = R * ((1 - model.α) * k_t**model.α - x)
-#     r = r * u_prime(y, θ)
-#     return l - r
-
-# def solve_for_k(x, θ, model, k_t):
-#     c = optimize.newton(solve_for_c, k_t, args=(θ, model, k_t, x))
-#     num = (1 - model.α) * k_t**model.α - c
-#     den = (1 + model.n)
-#     return x - num/den
-
-# def k_next(model, θ, k_t):
-#     return optimize.newton(solve_for_k, k_t, args=(θ, model, k_t), maxiter=200)
-
-# ts_length = 10
-# θ = 1.5
-# x0 = np.array([1.2, 2.6])
-
-# def simulate_ts(olg, θ, x0_values, ts_length):
+def find_k_next(k_t, model):
+    return optimize.newton(solve_for_k_t1, k_t, args=(k_t, model))
+```
 
-#     fig, ax = plt.subplots(figsize=(10, 5))
+```{code-cell} ipython3
+def simulate_ts(k0_values, model, ts_length=10):
 
-#     ts = np.zeros(ts_length)
+    fig, ax = plt.subplots(figsize=(10, 5))
 
-#     # simulate and plot time series
-#     for x_init in x0_values:
-#         ts[0] = x_init
-#         for t in range(1, ts_length):
-#             ts[t] = k_next(olg, θ, ts[t-1])
-#         ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6,
-#                 label=r'$k_0=%g$' %x_init)
-#     ax.plot(np.arange(ts_length), np.full(ts_length,k_star),
-#             alpha=0.6, color='red', label=r'$k^*$')
-#     ax.legend(fontsize=10)
+    ts = np.zeros(ts_length)
 
-#     ax.set_xlabel(r'$t$', fontsize=14)
-#     ax.set_ylabel(r'$k_t$', fontsize=14)
+    # simulate and plot time series
+    for x_init in k0_values:
+        ts[0] = x_init
+        for t in range(1, ts_length):
+            ts[t] = find_k_next(ts[t-1], model)
+        ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6,
+                label=r'$k_0=%g$' %x_init)
+    ax.legend(fontsize=10)
 
-#     plt.show()
+    ax.set_xlabel(r'$t$', fontsize=14)
+    ax.set_ylabel(r'$k_t$', fontsize=14)
 
-# simulate_ts(olg, θ, x0, ts_length)
+    plt.show()
 ```
 
 ```{code-cell} ipython3
+k0_values = [0.2, 10, 50, 100]
+m_quasilinear = create_olg_model(u=u_quasilinear)
+simulate_ts(k0_values, m_quasilinear)
+```
 
+```{solution-end}
 ```