From 1099889642de41385b51aed69d9e04ae948f7b98 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 9 May 2023 22:45:03 +1000
Subject: [PATCH 01/23] update newton_method

---
 lectures/newton_method.md | 243 +++++++++++++++++---------------------
 1 file changed, 107 insertions(+), 136 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 141161695..0d386370a 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -65,23 +67,21 @@ approximation and Newton's method.
 Then we apply Newton's method to multi-dimensional settings to solve
 market for equilibria with multiple goods.
 
-At the end of the lecture we leverage the power of JAX and automatic
-differentiation to solve a very high-dimensional equilibrium problem.
+At the end of the lecture we leverage the power of automatic
+differentiation in [`autograd`](https://github.com/HIPS/autograd) to solve a very high-dimensional equilibrium problem.
 
 We use the following imports in this lecture
 
-```{code-cell} python3
-import numpy as np
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
 from collections import namedtuple
 from scipy.optimize import root
-import jax
-import jax.numpy as jnp
+from autograd import jacobian
+import autograd.numpy as np
 
 plt.rcParams["figure.figsize"] = (10, 5.7)
 ```
 
-
 ## Fixed Point Computation Using Newton's Method
 
 In this section we solve the fixed point of the law of motion for capital in
@@ -126,13 +126,13 @@ $$ k^* = \left(\frac{s A}{δ}\right)^{1/(1 - α)}  $$
 
 Let's store our parameters in [`namedtuple`](https://docs.python.org/3/library/collections.html#collections.namedtuple) to help us keep our code clean and concise.
 
-```{code-cell} python3
+```{code-cell} ipython3
 SolowParameters = namedtuple("SolowParameters", ('A', 's', 'α', 'δ'))
 ```
 
 This function creates a suitable `namedtuple` with default parameter values.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def create_solow_params(A=2.0, s=0.3, α=0.3, δ=0.4):
     "Creates a Solow model parameterization with default values."
     return SolowParameters(A=A, s=s, α=α, δ=δ)
@@ -140,7 +140,7 @@ def create_solow_params(A=2.0, s=0.3, α=0.3, δ=0.4):
 
 The next two functions implement the law of motion [](motion_law) and store the true fixed point $k^*$.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def g(k, params):
     A, s, α, δ = params
     return A * s * k**α + (1 - δ) * k
@@ -149,9 +149,10 @@ def exact_fixed_point(params):
     A, s, α, δ = params
     return ((s * A) / δ)**(1/(1 - α))
 ```
+
 Here is a function to provide a 45 degree plot of the dynamics.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def plot_45(params, ax, fontsize=14):
     
     k_min, k_max = 0.0, 3.0
@@ -184,14 +185,14 @@ def plot_45(params, ax, fontsize=14):
 
 Let's look at the 45 degree diagram for two parameterizations.
 
-```{code-cell} python3
+```{code-cell} ipython3
 params = create_solow_params()
 fig, ax = plt.subplots(figsize=(8, 8))
 plot_45(params, ax)
 plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 params = create_solow_params(α=0.05, δ=0.5)
 fig, ax = plt.subplots(figsize=(8, 8))
 plot_45(params, ax)
@@ -210,8 +211,7 @@ from some initial state $k_0$ using the law of motion.
 
 Here's a time series from a particular choice of $k_0$.
 
-
-```{code-cell} python3
+```{code-cell} ipython3
 def compute_iterates(k_0, f, params, n=25):
     "Compute time series of length n generated by arbitrary function f."
     k = k_0
@@ -222,7 +222,7 @@ def compute_iterates(k_0, f, params, n=25):
     return k_iterates
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 params = create_solow_params()
 k_0 = 0.25
 k_series = compute_iterates(k_0, g, params)
@@ -237,7 +237,7 @@ plt.show()
 
 Let's see the output for a long time series.
 
-```{code-cell} python3
+```{code-cell} ipython3
 k_series = compute_iterates(k_0, g, params, n=10_000)
 k_star_approx = k_series[-1]
 k_star_approx
@@ -246,11 +246,11 @@ k_star_approx
 This is close to the true value.
 
 (solved_k)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 k_star
 ```
 
-
 #### Newton's Method 
 
 In general, when applying Newton's fixed point method to some function $g$, 
@@ -295,7 +295,7 @@ g'(k) = \alpha s A k^{\alpha-1} + (1-\delta)
 
 Let's define this:
 
-```{code-cell} python3
+```{code-cell} ipython3
 def Dg(k, params):
     A, s, α, δ = params
     return α * A * s * k**(α-1) + (1 - δ)
@@ -303,14 +303,14 @@ def Dg(k, params):
 
 Here's a function $q$ representing [](newtons_method).
 
-```{code-cell} python3
+```{code-cell} ipython3
 def q(k, params):
     return (g(k, params) - Dg(k, params) * k) / (1 - Dg(k, params))
 ```
 
 Now let's plot some trajectories.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def plot_trajectories(params, 
                       k0_a=0.8,  # first initial condition
                       k0_b=3.1,  # second initial condition
@@ -343,7 +343,7 @@ def plot_trajectories(params,
     plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 params = create_solow_params()
 plot_trajectories(params)
 ```
@@ -396,7 +396,8 @@ x_{t+1} = x_t - \frac{ f(x_t) }{ f'(x_t) },
 The following code implements the iteration [](oneD-newton)
 
 (first_newton_attempt)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 def newton(f, Df, x_0, tol=1e-7, max_iter=100_000):
     x = x_0
 
@@ -436,14 +437,14 @@ Any zero of $f$ is clearly a fixed point of $g$.
 
 Let's apply this idea to the Solow problem
 
-```{code-cell} python3
+```{code-cell} ipython3
 params = create_solow_params()
 k_star_approx_newton = newton(f=lambda x: g(x, params) - x,
                               Df=lambda x: Dg(x, params) - 1,
                               x_0=0.8)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 k_star_approx_newton
 ```
 
@@ -534,12 +535,11 @@ $$
 
 The function below calculates the excess demand for given parameters
 
-```{code-cell} python3
+```{code-cell} ipython3
 def e(p, A, b, c):
     return np.exp(- A @ p) + c - b * np.sqrt(p)
 ```
 
-
 Our default parameter values will be
 
 
@@ -560,8 +560,7 @@ A = \begin{pmatrix}
         \end{pmatrix}
 $$
 
-
-```{code-cell} python3
+```{code-cell} ipython3
 A = np.array([
     [0.5, 0.4],
     [0.8, 0.2]
@@ -570,23 +569,20 @@ b = np.ones(2)
 c = np.ones(2)
 ```
 
-At a price level of $p = (1, 0.5)$, the excess demand is 
-
+At a price level of $p = (1, 0.5)$, the excess demand is
 
-```{code-cell} python3
+```{code-cell} ipython3
 ex_demand = e((1.0, 0.5), A, b, c)
 
 print(f'The excess demand for good 0 is {ex_demand[0]:.3f} \n'
       f'The excess demand for good 1 is {ex_demand[1]:.3f}')
 ```
 
-
 Next we plot the two functions $e_0$ and $e_1$ on a grid of $(p_0, p_1)$ values, using contour surfaces and lines.
 
 We will use the following function to build the contour plots
 
-```{code-cell} python3
-
+```{code-cell} ipython3
 def plot_excess_demand(ax, good=0, grid_size=100, grid_max=4, surface=True):
 
     # Create a 100x100 grid
@@ -609,14 +605,16 @@ def plot_excess_demand(ax, good=0, grid_size=100, grid_max=4, surface=True):
 ```
 
 Here's our plot of $e_0$:
-```{code-cell} python3
+
+```{code-cell} ipython3
 fig, ax = plt.subplots()
 plot_excess_demand(ax, good=0)
 plt.show()
 ```
 
 Here's our plot of $e_1$:
-```{code-cell} python3
+
+```{code-cell} ipython3
 fig, ax = plt.subplots()
 plot_excess_demand(ax, good=1)
 plt.show()
@@ -628,8 +626,7 @@ For a price vector $p$ such that $e_i(p)=0$ we know that good $i$ is in equilibr
 
 If these two contour lines cross at some price vector $p^*$, then $p^*$ is an equilibrium price vector.
 
-
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(10, 5.7))
 for good in (0, 1):
     plot_excess_demand(ax, good=good, surface=False)
@@ -645,27 +642,27 @@ To solve for $p^*$ more precisely, we use a zero-finding algorithm from `scipy.o
 
 We supply $p = (1, 1)$ as our initial guess.
 
-```{code-cell} python3
+```{code-cell} ipython3
 init_p = np.ones(2)
 ```
 
 This uses the [modified Powell method](https://docs.scipy.org/doc/scipy/reference/optimize.root-hybr.html#optimize-root-hybr) to find the zero
 
-```{code-cell} python3
+```{code-cell} ipython3
 %%time
 solution = root(lambda p: e(p, A, b, c), init_p, method='hybr')
 ```
 
 Here's the resulting value:
 
-```{code-cell} python3
+```{code-cell} ipython3
 p = solution.x
 p
 ```
 
 This looks close to our guess from observing the figure. We can plug it back into $e$ to test that $e(p) \approx 0$:
 
-```{code-cell} python3
+```{code-cell} ipython3
 np.max(np.abs(e(p, A, b, c)))
 ```
 
@@ -686,8 +683,8 @@ J(p) =
     \end{pmatrix}
 $$
 
-```{code-cell} python3
-def jacobian(p, A, b, c):
+```{code-cell} ipython3
+def jacobian_e(p, A, b, c):
     p_0, p_1 = p
     a_00, a_01 = A[0, :]
     a_10, a_11 = A[1, :]
@@ -700,22 +697,21 @@ def jacobian(p, A, b, c):
     return np.array(J)
 ```
 
-```{code-cell} python3 
+```{code-cell} ipython3
 %%time
 solution = root(lambda p: e(p, A, b, c),
                 init_p, 
-                jac=lambda p: jacobian(p, A, b, c), 
+                jac=lambda p: jacobian_e(p, A, b, c), 
                 method='hybr')
 ```
 
 Now the solution is even more accurate (although, in this low-dimensional problem, the difference is quite small):
 
-```{code-cell} python3
+```{code-cell} ipython3
 p = solution.x
 np.max(np.abs(e(p, A, b, c)))
 ```
 
-
 #### Using Newton's Method
 
 Now let's use Newton's method to compute the equilibrium price using the multivariate version of Newton's method
@@ -732,14 +728,14 @@ This is a multivariate version of [](oneD-newton)
 
 The iteration starts from some initial guess of the price vector $p_0$. 
 
-Here, instead of coding Jacobian by hand, We use the `jax.jacobian()` function to auto-differentiate and calculate the Jacobian.
+Here, instead of coding Jacobian by hand, We use the `jacobian()` function in the `autograd` library to auto-differentiate and calculate the Jacobian.
 
 With only slight modification, we can generalize [our previous attempt](first_newton_attempt) to multi-dimensional problems
 
-```{code-cell} python3
+```{code-cell} ipython3
 def newton(f, x_0, tol=1e-5, max_iter=10):
     x = x_0
-    q = jax.jit(lambda x: x - jnp.linalg.solve(jax.jacobian(f)(x), f(x)))
+    q = lambda x: x - np.linalg.solve(jacobian(f)(x), f(x))
     error = tol + 1
     n = 0
     while error > tol:
@@ -747,28 +743,28 @@ def newton(f, x_0, tol=1e-5, max_iter=10):
         if(n > max_iter):
             raise Exception('Max iteration reached without convergence')
         y = q(x)
-        if(any(jnp.isnan(y))):
+        if(any(np.isnan(y))):
             raise Exception('Solution not found with NaN generated')
-        error = jnp.linalg.norm(x - y)
+        error = np.linalg.norm(x - y)
         x = y
         print(f'iteration {n}, error = {error:.5f}')
     print('\n' + f'Result = {x} \n')
     return x
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 def e(p, A, b, c):
-    return jnp.exp(- jnp.dot(A, p)) + c - b * jnp.sqrt(p)
+    return np.exp(- np.dot(A, p)) + c - b * np.sqrt(p)
 ```
 
 We find the algorithm terminates in 4 steps
 
-```{code-cell} python3
+```{code-cell} ipython3
 %%time
-p = newton(lambda p: e(p, A, b, c), init_p).block_until_ready()
+p = newton(lambda p: e(p, A, b, c), init_p)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 np.max(np.abs(e(p, A, b, c)))
 ```
 
@@ -784,66 +780,57 @@ However, things will change when we move to higher dimensional problems.
 
 Our next step is to investigate a large market with 5,000 goods.
 
-To handle this large problem we will use Google JAX.
+This section will be very slow to run. 
 
-The excess demand function is essentially the same, but now the matrix $A$ is $5000 \times 5000$ and the parameter vectors $b$ and $c$ are $5000 \times 1$.
+A JAX version of this lecture section using JAX accelerated linear algebra,
+automatic differentiation, and a GPU is available here TODO: link to JAX lecture
 
+(Spoiler: the JAX version can compute the result in just a few seconds)
+
+The excess demand function is essentially the same, but now the matrix $A$ is $5000 \times 5000$ and the parameter vectors $b$ and $c$ are $5000 \times 1$.
 
-```{code-cell} python3
-dim = 5_000
+```{code-cell} ipython3
+dim = 5000
 np.random.seed(123)
 
 # Create a random matrix A and normalize the rows to sum to one
 A = np.random.rand(dim, dim)
-A = jnp.asarray(A)
-s = jnp.sum(A, axis=0)
+A = np.asarray(A)
+s = np.sum(A, axis=0)
 A = A / s
 
 # Set up b and c
-b = jnp.ones(dim)
-c = jnp.ones(dim)
-```
-
-Here is essentially the same demand function we applied before, but now using `jax.numpy` for the calculations.
-
-```{code-cell} python3
-def e(p, A, b, c):
-    return jnp.exp(- jnp.dot(A, p)) + c - b * jnp.sqrt(p)
+b = np.ones(dim)
+c = np.ones(dim)
 ```
 
 Here's our initial condition
 
-```{code-cell} python3
-init_p = jnp.ones(dim)
+```{code-cell} ipython3
+init_p = np.ones(dim)
 ```
 
-By leveraging the power of Newton's method, JAX accelerated linear algebra,
-automatic differentiation, and a GPU, we obtain a relatively small error for
-this very large problem in just a few seconds:
-
-```{code-cell} python3
+```{code-cell} ipython3
 %%time
-p = newton(lambda p: e(p, A, b, c), init_p).block_until_ready()
+p = newton(lambda p: e(p, A, b, c), init_p)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 np.max(np.abs(e(p, A, b, c)))
 ```
 
-With the same tolerance, SciPy's `root` function takes much longer to run,
-even with the Jacobian supplied.
+With the same tolerance, SciPy's `root` function takes about the same time, but is less accurate
 
-
-```{code-cell} python3
+```{code-cell} ipython3
 %%time
 solution = root(lambda p: e(p, A, b, c),
                 init_p, 
-                jac=lambda p: jax.jacobian(e)(p, A, b, c), 
+                jac=lambda p: jacobian(e)(p, A, b, c), 
                 method='hybr',
                 tol=1e-5)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 p = solution.x
 np.max(np.abs(e(p, A, b, c)))
 ```
@@ -924,36 +911,36 @@ The result should converge to the [analytical solution](solved_k).
 
 Let's first define the parameters for this problem
 
-```{code-cell} python3
-A = jnp.array([[2.0, 3.0, 3.0],
-               [2.0, 4.0, 2.0],
-               [1.0, 5.0, 1.0]])
+```{code-cell} ipython3
+A = np.array([[2.0, 3.0, 3.0],
+              [2.0, 4.0, 2.0],
+              [1.0, 5.0, 1.0]])
 
 s = 0.2
 α = 0.5
 δ = 0.8
 
-initLs = [jnp.ones(3),
-          jnp.array([3.0, 5.0, 5.0]),
-          jnp.repeat(50.0, 3)]
+initLs = [np.ones(3),
+          np.array([3.0, 5.0, 5.0]),
+          np.repeat(50.0, 3)]
 ```
 
 Then define the multivariate version of the formula for the [law of motion of captial](motion_law)
 
-```{code-cell} python3
+```{code-cell} ipython3
 def multivariate_solow(k, A=A, s=s, α=α, δ=δ):
-    return (s * jnp.dot(A, k**α) + (1 - δ) * k)
+    return (s * np.dot(A, k**α) + (1 - δ) * k)
 ```
 
 Let's run through each starting value and see the output
 
-```{code-cell} python3
+```{code-cell} ipython3
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
 
     %time k = newton(lambda k: multivariate_solow(k) - k, \
-                     init).block_until_ready()
+                     init)
     print('-'*64)
     attempt +=1
 ```
@@ -964,7 +951,7 @@ But the number of iterations it takes to converge is dependent on the starting v
 
 Let substitute the output back to the formulate to check our last result
 
-```{code-cell} python3
+```{code-cell} ipython3
 multivariate_solow(k) - k
 ```
 
@@ -972,8 +959,8 @@ Note the error is very small.
 
 We can also test our results on the known solution
 
-```{code-cell} python3
-A = jnp.array([[2.0, 0.0, 0.0],
+```{code-cell} ipython3
+A = np.array([[2.0, 0.0, 0.0],
                [0.0, 2.0, 0.0],
                [0.0, 0.0, 2.0]])
 
@@ -981,28 +968,20 @@ s = 0.3
 α = 0.3
 δ = 0.4
 
-init = jnp.repeat(1.0, 3)
+init = np.repeat(1.0, 3)
 
 
 %time k = newton(lambda k: multivariate_solow(k, A=A, s=s, α=α, δ=δ) - k, \
-                 init).block_until_ready()
+                 init)
 ```
 
 The result is very close to the ground truth but still slightly different.
 
-We can increase the precision of the floating point numbers and restrict the tolerance to obtain a more accurate approximation (see detailed discussion in the [lecture on JAX](https://python-programming.quantecon.org/jax_intro.html#differences))
-
-```{code-cell} python3
-from jax.config import config
-
-config.update("jax_enable_x64", True)
-
-init = init.astype('float64')
-
+```{code-cell} ipython3
 %time k = newton(lambda k: multivariate_solow(k, A=A, s=s, α=α, δ=δ) - k, \
                  init,\
-                 tol=1e-7).block_until_ready()
-```   
+                 tol=1e-7)
+```
 
 We can see it steps towards a more accurate solution.
 
@@ -1051,14 +1030,6 @@ $$
 
 Set the tolerance to $0.0$ for more accurate output.
 
-
-```{hint} 
-:class: dropdown
-
-Similar to [exercise 1](newton_ex1), enabling `float64` for JAX can improve the precision of our results.
-```
-
-
 ```{exercise-end}
 ```
 
@@ -1068,24 +1039,24 @@ Similar to [exercise 1](newton_ex1), enabling `float64` for JAX can improve the
 
 Define parameters and initial values
 
-```{code-cell} python3
-A = jnp.array([
+```{code-cell} ipython3
+A = np.array([
     [0.2, 0.1, 0.7],
     [0.3, 0.2, 0.5],
     [0.1, 0.8, 0.1]
 ])
 
-b = jnp.array([1.0, 1.0, 1.0])
-c = jnp.array([1.0, 1.0, 1.0])
+b = np.array([1.0, 1.0, 1.0])
+c = np.array([1.0, 1.0, 1.0])
 
-initLs = [jnp.repeat(5.0, 3),
-          jnp.ones(3),
-          jnp.array([4.5, 0.1, 4.0])] 
+initLs = [np.repeat(5.0, 3),
+          np.ones(3),
+          np.array([4.5, 0.1, 4.0])] 
 ```
 
 Let’s run through each initial guess and check the output
 
-```{code-cell} python3
+```{code-cell} ipython3
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
@@ -1094,7 +1065,7 @@ for init in initLs:
 
     %time p = newton(lambda p: e(p, A, b, c), \
                  init, \
-                 tol=0.0).block_until_ready()
+                 tol=0.0)
     print('-'*64)
     attempt +=1
 ```
@@ -1105,7 +1076,7 @@ Sometimes it may take a few initial guesses to achieve convergence.
 
 Substitute the result back to the formula to check our result
 
-```{code-cell} python3
+```{code-cell} ipython3
 e(p, A, b, c)
 ```
 

From 2b08451543db1c539cc2dc5037f9e35dd2c735a7 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 9 May 2023 22:49:11 +1000
Subject: [PATCH 02/23] add a comment on import

---
 lectures/newton_method.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 0d386370a..9926c91fe 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -77,6 +77,7 @@ import matplotlib.pyplot as plt
 from collections import namedtuple
 from scipy.optimize import root
 from autograd import jacobian
+# Thinly-wrapped numpy to enable automatic differentiation
 import autograd.numpy as np
 
 plt.rcParams["figure.figsize"] = (10, 5.7)

From d35e20e209411096556a777a9588676f6af61ecc Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 9 May 2023 23:47:55 +1000
Subject: [PATCH 03/23] update seealso at the beginning of the lecture

---
 lectures/newton_method.md | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 9926c91fe..3e357a96c 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -28,6 +28,11 @@ kernelspec:
 :depth: 2
 ```
 
+```{seealso}
+**GPU:** A version of this lecture which makes use of [jax](https://jax.readthedocs.io) to run the code
+on a `GPU` is [available here](https://jax.quantecon.org/newtons_method.html)
+```
+
 ## Overview
 
 Many economic problems involve finding [fixed
@@ -779,19 +784,17 @@ However, things will change when we move to higher dimensional problems.
 
 ### A High-Dimensional Problem
 
-Our next step is to investigate a large market with 5,000 goods.
-
-This section will be very slow to run. 
+Our next step is to investigate a large market with 3,000 goods.
 
-A JAX version of this lecture section using JAX accelerated linear algebra,
-automatic differentiation, and a GPU is available here TODO: link to JAX lecture
+A JAX version of this section using GPU accelerated linear algebra,
+automatic differentiation is available [here](https://jax.quantecon.org/newtons_method.html#application)
 
 (Spoiler: the JAX version can compute the result in just a few seconds)
 
-The excess demand function is essentially the same, but now the matrix $A$ is $5000 \times 5000$ and the parameter vectors $b$ and $c$ are $5000 \times 1$.
+The excess demand function is essentially the same, but now the matrix $A$ is $3000 \times 3000$ and the parameter vectors $b$ and $c$ are $3000 \times 1$.
 
 ```{code-cell} ipython3
-dim = 5000
+dim = 3000
 np.random.seed(123)
 
 # Create a random matrix A and normalize the rows to sum to one
@@ -820,7 +823,7 @@ p = newton(lambda p: e(p, A, b, c), init_p)
 np.max(np.abs(e(p, A, b, c)))
 ```
 
-With the same tolerance, SciPy's `root` function takes about the same time, but is less accurate
+With the same tolerance, SciPy's `root` function is slightly faster, but is less accurate
 
 ```{code-cell} ipython3
 %%time

From ac5063d541d9d29e50bc8e2badd6de3f6b691669 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Wed, 10 May 2023 00:19:26 +1000
Subject: [PATCH 04/23] add pip install

---
 lectures/newton_method.md | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 3e357a96c..fcfdb0118 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -73,7 +73,12 @@ Then we apply Newton's method to multi-dimensional settings to solve
 market for equilibria with multiple goods.
 
 At the end of the lecture we leverage the power of automatic
-differentiation in [`autograd`](https://github.com/HIPS/autograd) to solve a very high-dimensional equilibrium problem.
+differentiation in [`autograd`](https://github.com/HIPS/autograd) to solve a very high-dimensional equilibrium problem
+
+```{code-cell} ipython3
+!pip install autograd
+```
+
 
 We use the following imports in this lecture
 

From 2a5705f590b38db39983cdb84b811b56af50400a Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Wed, 10 May 2023 01:02:34 +1000
Subject: [PATCH 05/23] update to accept error

---
 lectures/newton_method.md | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index fcfdb0118..2f9d832fd 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -76,10 +76,11 @@ At the end of the lecture we leverage the power of automatic
 differentiation in [`autograd`](https://github.com/HIPS/autograd) to solve a very high-dimensional equilibrium problem
 
 ```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install autograd
 ```
 
-
 We use the following imports in this lecture
 
 ```{code-cell} ipython3
@@ -944,6 +945,8 @@ def multivariate_solow(k, A=A, s=s, α=α, δ=δ):
 Let's run through each starting value and see the output
 
 ```{code-cell} ipython3
+:tags: ["raises-exception"]
+
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
@@ -1066,6 +1069,8 @@ initLs = [np.repeat(5.0, 3),
 Let’s run through each initial guess and check the output
 
 ```{code-cell} ipython3
+:tags: ["raises-exception"]
+
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')

From 1e813e70229a879bd7e6bb602ce8f6cea94e4cd1 Mon Sep 17 00:00:00 2001
From: mmcky <mamckay@gmail.com>
Date: Wed, 10 May 2023 12:58:35 +1000
Subject: [PATCH 06/23] TMP: remove cache for a full build

---
 .github/workflows/ci.yml | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 96e909895..0703ea3b5 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -40,13 +40,13 @@ jobs:
       - name: Display Pip Versions
         shell: bash -l {0}
         run: pip list
-      - name: Download "build" folder (cache)
-        uses: dawidd6/action-download-artifact@v2
-        with:
-          workflow: cache.yml
-          branch: main
-          name: build-cache
-          path: _build
+      # - name: Download "build" folder (cache)
+      #   uses: dawidd6/action-download-artifact@v2
+      #   with:
+      #     workflow: cache.yml
+      #     branch: main
+      #     name: build-cache
+      #     path: _build
       # Build Assets (Download Notebooks and PDF via LaTeX)
       - name: Build Download Notebooks (sphinx-tojupyter)
         shell: bash -l {0}

From c6fc30c304d0e65927a2beebd1f5f998cce35a34 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Mon, 15 May 2023 13:49:06 +1000
Subject: [PATCH 07/23] fix bugs in career and mccall model

---
 lectures/career.md       | 39 +++++++++++++++++----------------
 lectures/mccall_model.md | 47 ++++++++++++++++++++--------------------
 2 files changed, 44 insertions(+), 42 deletions(-)

diff --git a/lectures/career.md b/lectures/career.md
index 54e0f8282..7a8b33c69 100644
--- a/lectures/career.md
+++ b/lectures/career.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -29,10 +31,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
 ```
 
@@ -46,7 +47,7 @@ This exposition draws on the presentation in {cite}`Ljungqvist2012`, section 6.5
 
 We begin with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 %matplotlib inline
 import matplotlib.pyplot as plt
 plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
@@ -165,7 +166,7 @@ Nice properties:
 
 Here's a figure showing the effect on the pmf of different shape parameters when $n=50$.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def gen_probs(n, a, b):
     probs = np.zeros(n+1)
     for k in range(n+1):
@@ -188,7 +189,7 @@ plt.show()
 We will first create a class `CareerWorkerProblem` which will hold the
 default parameterizations of the model and an initial guess for the value function.
 
-```{code-cell} python3
+```{code-cell} ipython3
 class CareerWorkerProblem:
 
     def __init__(self,
@@ -221,7 +222,7 @@ the corresponding Bellman operator $T$ and the greedy policy function.
 In this model, $T$ is defined by $Tv(\theta, \epsilon) = \max\{I, II, III\}$, where
 $I$, $II$ and $III$ are as given in {eq}`eyes`.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def operator_factory(cw, parallel_flag=True):
 
     """
@@ -277,7 +278,7 @@ def operator_factory(cw, parallel_flag=True):
 Lastly, `solve_model` will  take an instance of `CareerWorkerProblem` and
 iterate using the Bellman operator to find the fixed point of the Bellman equation.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def solve_model(cw,
                 use_parallel=True,
                 tol=1e-4,
@@ -311,7 +312,7 @@ def solve_model(cw,
 
 Here's the solution to the model -- an approximate value function
 
-```{code-cell} python3
+```{code-cell} ipython3
 cw = CareerWorkerProblem()
 T, get_greedy = operator_factory(cw)
 v_star = solve_model(cw, verbose=False)
@@ -333,7 +334,7 @@ plt.show()
 
 And here is the optimal policy
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(6, 6))
 tg, eg = np.meshgrid(cw.θ, cw.ϵ)
 lvls = (0.5, 1.5, 2.5, 3.5)
@@ -393,7 +394,7 @@ In reading the code, recall that `optimal_policy[i, j]` = policy at
 $(\theta_i, \epsilon_j)$ = either 1, 2 or 3; meaning 'stay put',
 'new job' and 'new life'.
 
-```{code-cell} python3
+```{code-cell} ipython3
 F = np.cumsum(cw.F_probs)
 G = np.cumsum(cw.G_probs)
 v_star = solve_model(cw, verbose=False)
@@ -465,7 +466,7 @@ Repeat the exercise with $\beta=0.99$ and interpret the change.
 
 The median for the original parameterization can be computed as follows
 
-```{code-cell} python3
+```{code-cell} ipython3
 cw = CareerWorkerProblem()
 F = np.cumsum(cw.F_probs)
 G = np.cumsum(cw.G_probs)
@@ -481,12 +482,12 @@ def passage_time(optimal_policy, F, G):
         if optimal_policy[i, j] == 1:    # Stay put
             return t
         elif optimal_policy[i, j] == 2:  # New job
-            j = qe.random.draw(G)
+            j = qe.random.draw(G, size=1)[0]
         else:                            # New life
-            i, j  = qe.random.draw(F), qe.random.draw(G)
+            i  = qe.random.draw(F, size=1)[0]
+            j  = qe.random.draw(G, size=1)[0]
         t += 1
-
-@njit(parallel=True)
+@njit
 def median_time(optimal_policy, F, G, M=25000):
     samples = np.empty(M)
     for i in prange(M):
@@ -521,7 +522,7 @@ figure -- interpret.
 
 Here is one solution
 
-```{code-cell} python3
+```{code-cell} ipython3
 cw = CareerWorkerProblem(G_a=100, G_b=100)
 T, get_greedy = operator_factory(cw)
 v_star = solve_model(cw, verbose=False)
diff --git a/lectures/mccall_model.md b/lectures/mccall_model.md
index b12d87b4f..0c286f52a 100644
--- a/lectures/mccall_model.md
+++ b/lectures/mccall_model.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.4
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -34,10 +36,9 @@ and the pros and cons as they themselves see them." -- Robert E. Lucas, Jr.
 
 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
 ```
 
@@ -59,7 +60,7 @@ As we'll see, McCall's model is not only interesting in its own right but also a
 
 Let's start with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 %matplotlib inline
 import matplotlib.pyplot as plt
 plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
@@ -343,21 +344,21 @@ $v$.
 Our default for $q$, the distribution of the state process, will be
 [Beta-binomial](https://en.wikipedia.org/wiki/Beta-binomial_distribution).
 
-```{code-cell} python3
+```{code-cell} ipython3
 n, a, b = 50, 200, 100                        # default parameters
 q_default = BetaBinomial(n, a, b).pdf()       # default choice of q
 ```
 
 Our default set of values for wages will be
 
-```{code-cell} python3
+```{code-cell} ipython3
 w_min, w_max = 10, 60
 w_default = np.linspace(w_min, w_max, n+1)
 ```
 
 Here's a plot of the probabilities of different wage outcomes:
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots()
 ax.plot(w_default, q_default, '-o', label='$q(w(i))$')
 ax.set_xlabel('wages')
@@ -372,7 +373,7 @@ We are going to use Numba to accelerate our code.
 
 The following helps Numba by providing some type
 
-```{code-cell} python3
+```{code-cell} ipython3
 mccall_data = [
     ('c', float64),      # unemployment compensation
     ('β', float64),      # discount factor
@@ -387,7 +388,7 @@ given the current state and an arbitrary feasible action.
 
 Default parameter values are embedded in the class.
 
-```{code-cell} python3
+```{code-cell} ipython3
 @jitclass(mccall_data)
 class McCallModel:
 
@@ -417,7 +418,7 @@ We will start from guess $v$ given by $v(i) = w(i) / (1 - β)$, which is the val
 
 Here's a function to implement this:
 
-```{code-cell} python3
+```{code-cell} ipython3
 def plot_value_function_seq(mcm, ax, num_plots=6):
     """
     Plot a sequence of value functions.
@@ -442,7 +443,7 @@ def plot_value_function_seq(mcm, ax, num_plots=6):
 
 Now let's create an instance of `McCallModel` and watch iterations  $T^k v$ converge from below:
 
-```{code-cell} python3
+```{code-cell} ipython3
 mcm = McCallModel()
 
 fig, ax = plt.subplots()
@@ -461,7 +462,7 @@ the reservation wage.
 
 We'll be using JIT compilation via Numba to turbocharge our loops.
 
-```{code-cell} python3
+```{code-cell} ipython3
 @jit(nopython=True)
 def compute_reservation_wage(mcm,
                              max_iter=500,
@@ -494,7 +495,7 @@ def compute_reservation_wage(mcm,
 
 The next line computes the reservation wage at  default parameters
 
-```{code-cell} python3
+```{code-cell} ipython3
 compute_reservation_wage(mcm)
 ```
 
@@ -506,7 +507,7 @@ parameters.
 In particular, let's look at what happens when we change $\beta$ and
 $c$.
 
-```{code-cell} python3
+```{code-cell} ipython3
 grid_size = 25
 R = np.empty((grid_size, grid_size))
 
@@ -519,7 +520,7 @@ for i, c in enumerate(c_vals):
         R[i, j] = compute_reservation_wage(mcm)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots()
 
 cs1 = ax.contourf(c_vals, β_vals, R.T, alpha=0.75)
@@ -612,7 +613,7 @@ The big difference here, however, is that we're iterating on a scalar $h$, rathe
 
 Here's an implementation:
 
-```{code-cell} python3
+```{code-cell} ipython3
 @jit(nopython=True)
 def compute_reservation_wage_two(mcm,
                                  max_iter=500,
@@ -730,7 +731,7 @@ Once your code is working, investigate how the reservation wage changes with $c$
 
 Here's one solution
 
-```{code-cell} python3
+```{code-cell} ipython3
 cdf = np.cumsum(q_default)
 
 @jit(nopython=True)
@@ -740,7 +741,7 @@ def compute_stopping_time(w_bar, seed=1234):
     t = 1
     while True:
         # Generate a wage draw
-        w = w_default[qe.random.draw(cdf)]
+        w = w_default[qe.random.draw(cdf, size=1)]
         # Stop when the draw is above the reservation wage
         if w >= w_bar:
             stopping_time = t
@@ -782,7 +783,7 @@ plt.show()
 
 Here is one solution:
 
-```{code-cell} python3
+```{code-cell} ipython3
 mccall_data_continuous = [
     ('c', float64),          # unemployment compensation
     ('β', float64),          # discount factor
@@ -832,7 +833,7 @@ $\beta$.
 
 We will do this using a contour plot.
 
-```{code-cell} python3
+```{code-cell} ipython3
 grid_size = 25
 R = np.empty((grid_size, grid_size))
 
@@ -845,7 +846,7 @@ for i, c in enumerate(c_vals):
         R[i, j] = compute_reservation_wage_continuous(mcmc)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots()
 
 cs1 = ax.contourf(c_vals, β_vals, R.T, alpha=0.75)

From d6ddfcfab698de00a0a3e97faca322c70c07b005 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Mon, 15 May 2023 18:03:32 +1000
Subject: [PATCH 08/23] revert changes in CI

---
 .github/workflows/ci.yml  | 14 +++++++-------
 lectures/newton_method.md |  4 ++--
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0703ea3b5..96e909895 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -40,13 +40,13 @@ jobs:
       - name: Display Pip Versions
         shell: bash -l {0}
         run: pip list
-      # - name: Download "build" folder (cache)
-      #   uses: dawidd6/action-download-artifact@v2
-      #   with:
-      #     workflow: cache.yml
-      #     branch: main
-      #     name: build-cache
-      #     path: _build
+      - name: Download "build" folder (cache)
+        uses: dawidd6/action-download-artifact@v2
+        with:
+          workflow: cache.yml
+          branch: main
+          name: build-cache
+          path: _build
       # Build Assets (Download Notebooks and PDF via LaTeX)
       - name: Build Download Notebooks (sphinx-tojupyter)
         shell: bash -l {0}
diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 2f9d832fd..2642c59b3 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -945,7 +945,7 @@ def multivariate_solow(k, A=A, s=s, α=α, δ=δ):
 Let's run through each starting value and see the output
 
 ```{code-cell} ipython3
-:tags: ["raises-exception"]
+:tags: [raises-exception]
 
 attempt = 1
 for init in initLs:
@@ -1069,7 +1069,7 @@ initLs = [np.repeat(5.0, 3),
 Let’s run through each initial guess and check the output
 
 ```{code-cell} ipython3
-:tags: ["raises-exception"]
+:tags: [raises-exception]
 
 attempt = 1
 for init in initLs:

From a829aeb1013a5b0aaabe481dad76ece57e327d1d Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Mon, 15 May 2023 20:51:38 +1000
Subject: [PATCH 09/23] remove the time magic

---
 lectures/newton_method.md | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 2642c59b3..2f7dbb169 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -945,7 +945,6 @@ def multivariate_solow(k, A=A, s=s, α=α, δ=δ):
 Let's run through each starting value and see the output
 
 ```{code-cell} ipython3
-:tags: [raises-exception]
 
 attempt = 1
 for init in initLs:
@@ -1069,15 +1068,13 @@ initLs = [np.repeat(5.0, 3),
 Let’s run through each initial guess and check the output
 
 ```{code-cell} ipython3
-:tags: [raises-exception]
-
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
 
     init = init.astype('float64')
 
-    %time p = newton(lambda p: e(p, A, b, c), \
+    p = newton(lambda p: e(p, A, b, c), \
                  init, \
                  tol=0.0)
     print('-'*64)

From 1aa6b9909871ed24ec1b55e53b03ee84ea7c2e22 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Mon, 15 May 2023 23:28:42 +1000
Subject: [PATCH 10/23] use try-catch instead relying on myst

---
 lectures/newton_method.md | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 2f7dbb169..4073e83b5 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -945,15 +945,16 @@ def multivariate_solow(k, A=A, s=s, α=α, δ=δ):
 Let's run through each starting value and see the output
 
 ```{code-cell} ipython3
-
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
-
-    %time k = newton(lambda k: multivariate_solow(k) - k, \
-                     init)
+    try:
+        %time k = newton(lambda k: multivariate_solow(k) - k, \
+                        init)
+    except Exception:
+        print('This iteration failed to converge')
     print('-'*64)
-    attempt +=1
+    attempt += 1
 ```
 
 We find that the results are invariant to the starting values given the well-defined property of this question.
@@ -1073,12 +1074,14 @@ for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
 
     init = init.astype('float64')
-
-    p = newton(lambda p: e(p, A, b, c), \
-                 init, \
-                 tol=0.0)
+    try:
+       %time p = newton(lambda p: e(p, A, b, c), \
+                    init, \
+                    tol=0.0)
+    except Exception:
+       print('This iteration failed to converge')
     print('-'*64)
-    attempt +=1
+    attempt += 1
 ```
 
 We can find that Newton's method may fail for some starting values.

From ef994e03929ed01598581be3c476faee29940816 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Mon, 15 May 2023 23:52:55 +1000
Subject: [PATCH 11/23] add raises-exception tag

---
 lectures/newton_method.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 4073e83b5..86da1e1e5 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1069,6 +1069,10 @@ initLs = [np.repeat(5.0, 3),
 Let’s run through each initial guess and check the output
 
 ```{code-cell} ipython3
+---
+tags: [raises-exception]
+---
+
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')

From 1cdb306af2d9789d665f81dd1ad17565b393b8d0 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 09:15:59 +1000
Subject: [PATCH 12/23] catch all errors

---
 lectures/newton_method.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 86da1e1e5..87e43ffcc 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -951,7 +951,7 @@ for init in initLs:
     try:
         %time k = newton(lambda k: multivariate_solow(k) - k, \
                         init)
-    except Exception:
+    except:
         print('This iteration failed to converge')
     print('-'*64)
     attempt += 1
@@ -1082,7 +1082,7 @@ for init in initLs:
        %time p = newton(lambda p: e(p, A, b, c), \
                     init, \
                     tol=0.0)
-    except Exception:
+    except:
        print('This iteration failed to converge')
     print('-'*64)
     attempt += 1

From e97d187461a87c6aca8617218e6227d1c934a97d Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 10:13:51 +1000
Subject: [PATCH 13/23] add a clear example for the function

---
 lectures/newton_method.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 87e43ffcc..0cd937a1b 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1095,6 +1095,10 @@ Sometimes it may take a few initial guesses to achieve convergence.
 Substitute the result back to the formula to check our result
 
 ```{code-cell} ipython3
+init = np.array([4.5, 0.1, 4.0])
+p = newton(lambda p: e(p, A, b, c), \
+                    init, \
+                    tol=0.0)
 e(p, A, b, c)
 ```
 

From 734410cc4dee6aab852a6af42bfc9de849183804 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 10:50:02 +1000
Subject: [PATCH 14/23] check output for the solution

---
 lectures/newton_method.md | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 0cd937a1b..d2b402a34 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1063,7 +1063,8 @@ c = np.array([1.0, 1.0, 1.0])
 
 initLs = [np.repeat(5.0, 3),
           np.ones(3),
-          np.array([4.5, 0.1, 4.0])] 
+          np.array([4.5, 0.1, 4.0]),
+          np.array([0, 1, 0])] 
 ```
 
 Let’s run through each initial guess and check the output
@@ -1094,14 +1095,6 @@ Sometimes it may take a few initial guesses to achieve convergence.
 
 Substitute the result back to the formula to check our result
 
-```{code-cell} ipython3
-init = np.array([4.5, 0.1, 4.0])
-p = newton(lambda p: e(p, A, b, c), \
-                    init, \
-                    tol=0.0)
-e(p, A, b, c)
-```
-
 We can see the result is very accurate.
 
 ```{solution-end}

From dc7d1e5597a6710541a44dd67bc88b098b04182d Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 11:44:06 +1000
Subject: [PATCH 15/23] add max iter

---
 lectures/newton_method.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index d2b402a34..cf5368a67 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1063,8 +1063,7 @@ c = np.array([1.0, 1.0, 1.0])
 
 initLs = [np.repeat(5.0, 3),
           np.ones(3),
-          np.array([4.5, 0.1, 4.0]),
-          np.array([0, 1, 0])] 
+          np.array([4.5, 0.1, 4.0])] 
 ```
 
 Let’s run through each initial guess and check the output
@@ -1080,9 +1079,10 @@ for init in initLs:
 
     init = init.astype('float64')
     try:
-       %time p = newton(lambda p: e(p, A, b, c), \
-                    init, \
-                    tol=0.0)
+       %time p = newton(lambda p: e(p, A, b, c), 
+                    init,
+                    tol=0.0,
+                    max_iter=15)
     except:
        print('This iteration failed to converge')
     print('-'*64)

From 1e988ac923d6d7d91b362d5bf0d9622b54961c31 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 14:52:33 +1000
Subject: [PATCH 16/23] update max_iter

---
 lectures/newton_method.md | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index cf5368a67..e35263173 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1079,9 +1079,9 @@ for init in initLs:
 
     init = init.astype('float64')
     try:
-       %time p = newton(lambda p: e(p, A, b, c), 
-                    init,
-                    tol=0.0,
+       %time p = newton(lambda p: e(p, A, b, c), \
+                    init, \
+                    tol=0.0, \
                     max_iter=15)
     except:
        print('This iteration failed to converge')
@@ -1095,6 +1095,10 @@ Sometimes it may take a few initial guesses to achieve convergence.
 
 Substitute the result back to the formula to check our result
 
+```{code-cell} ipython3
+e(p, A, b, c)
+```
+
 We can see the result is very accurate.
 
 ```{solution-end}

From 801ea6f693678c1c2ff575c3c1342098142b951a Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 15:09:17 +1000
Subject: [PATCH 17/23] check the output of the line before

---
 lectures/newton_method.md | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index e35263173..2b80aec14 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1095,10 +1095,6 @@ Sometimes it may take a few initial guesses to achieve convergence.
 
 Substitute the result back to the formula to check our result
 
-```{code-cell} ipython3
-e(p, A, b, c)
-```
-
 We can see the result is very accurate.
 
 ```{solution-end}

From 5e6dfb455cfaa6f570e5de8804a34c7bf3405371 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 15:25:40 +1000
Subject: [PATCH 18/23] fix convergence issue

---
 lectures/newton_method.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 2b80aec14..24d7c92f6 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -1081,7 +1081,7 @@ for init in initLs:
     try:
        %time p = newton(lambda p: e(p, A, b, c), \
                     init, \
-                    tol=0.0, \
+                    tol=1e-15, \
                     max_iter=15)
     except:
        print('This iteration failed to converge')

From 3d0456d626eed2934637bb377439a8f9d5cce07b Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 16 May 2023 15:49:34 +1000
Subject: [PATCH 19/23] finalize the code

---
 lectures/newton_method.md | 16 ++++++----------
 1 file changed, 6 insertions(+), 10 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 24d7c92f6..f763d1325 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -784,19 +784,14 @@ The result is very accurate.
 
 With the larger overhead, the speed is not better than the optimized `scipy` function.
 
-However, things will change when we move to higher dimensional problems.
-
-
 
 ### A High-Dimensional Problem
 
 Our next step is to investigate a large market with 3,000 goods.
 
-A JAX version of this section using GPU accelerated linear algebra,
+A JAX version of this section using GPU accelerated linear algebra and
 automatic differentiation is available [here](https://jax.quantecon.org/newtons_method.html#application)
 
-(Spoiler: the JAX version can compute the result in just a few seconds)
-
 The excess demand function is essentially the same, but now the matrix $A$ is $3000 \times 3000$ and the parameter vectors $b$ and $c$ are $3000 \times 1$.
 
 ```{code-cell} ipython3
@@ -829,7 +824,7 @@ p = newton(lambda p: e(p, A, b, c), init_p)
 np.max(np.abs(e(p, A, b, c)))
 ```
 
-With the same tolerance, SciPy's `root` function is slightly faster, but is less accurate
+With the same tolerance, SciPy's `root` function is slightly faster, but the result is less accurate
 
 ```{code-cell} ipython3
 %%time
@@ -845,9 +840,6 @@ p = solution.x
 np.max(np.abs(e(p, A, b, c)))
 ```
 
-The result is also less accurate.
-
-
 
 ## Exercises
 
@@ -1095,6 +1087,10 @@ Sometimes it may take a few initial guesses to achieve convergence.
 
 Substitute the result back to the formula to check our result
 
+```{code-cell} ipython3
+e(p, A, b, c)
+```
+
 We can see the result is very accurate.
 
 ```{solution-end}

From 5a14609bc1cc22939733aaaf653fa0d335e4acfc Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 23 May 2023 16:43:24 +1000
Subject: [PATCH 20/23] remove try-catch to raise errors

---
 lectures/newton_method.md | 20 ++++++--------------
 1 file changed, 6 insertions(+), 14 deletions(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index f763d1325..139950e30 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -940,11 +940,8 @@ Let's run through each starting value and see the output
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
-    try:
-        %time k = newton(lambda k: multivariate_solow(k) - k, \
-                        init)
-    except:
-        print('This iteration failed to converge')
+    %time k = newton(lambda k: multivariate_solow(k) - k, \
+                    init)
     print('-'*64)
     attempt += 1
 ```
@@ -1068,15 +1065,10 @@ tags: [raises-exception]
 attempt = 1
 for init in initLs:
     print(f'Attempt {attempt}: Starting value is {init} \n')
-
-    init = init.astype('float64')
-    try:
-       %time p = newton(lambda p: e(p, A, b, c), \
-                    init, \
-                    tol=1e-15, \
-                    max_iter=15)
-    except:
-       print('This iteration failed to converge')
+    %time p = newton(lambda p: e(p, A, b, c), \
+                init, \
+                tol=1e-15, \
+                max_iter=15)
     print('-'*64)
     attempt += 1
 ```

From a0b588a2360a45d80e8cbda68ca4af64ab587138 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <u6474961@anu.edu.au>
Date: Tue, 23 May 2023 17:10:04 +1000
Subject: [PATCH 21/23] update the runtime comparison

---
 lectures/newton_method.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lectures/newton_method.md b/lectures/newton_method.md
index 139950e30..b6651a38b 100644
--- a/lectures/newton_method.md
+++ b/lectures/newton_method.md
@@ -824,7 +824,7 @@ p = newton(lambda p: e(p, A, b, c), init_p)
 np.max(np.abs(e(p, A, b, c)))
 ```
 
-With the same tolerance, SciPy's `root` function is slightly faster, but the result is less accurate
+With the same tolerance, we compare the runtime and accuracy of Newton's method to SciPy's `root` function
 
 ```{code-cell} ipython3
 %%time

From 8a475c67dbc401edc2a79ef1b62f9bd0b8d66543 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <humphrey.yang@anu.edu.au>
Date: Tue, 30 May 2023 22:01:20 +1000
Subject: [PATCH 22/23] revert temprary bug fixes

---
 lectures/career.md       | 39 ++++++++++++++++-----------------
 lectures/mccall_model.md | 47 ++++++++++++++++++++--------------------
 2 files changed, 42 insertions(+), 44 deletions(-)

diff --git a/lectures/career.md b/lectures/career.md
index 7a8b33c69..54e0f8282 100644
--- a/lectures/career.md
+++ b/lectures/career.md
@@ -3,10 +3,8 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.14.4
 kernelspec:
-  display_name: Python 3 (ipykernel)
+  display_name: Python 3
   language: python
   name: python3
 ---
@@ -31,9 +29,10 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython3
-:tags: [hide-output]
-
+```{code-cell} ipython
+---
+tags: [hide-output]
+---
 !pip install quantecon
 ```
 
@@ -47,7 +46,7 @@ This exposition draws on the presentation in {cite}`Ljungqvist2012`, section 6.5
 
 We begin with some imports:
 
-```{code-cell} ipython3
+```{code-cell} ipython
 %matplotlib inline
 import matplotlib.pyplot as plt
 plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
@@ -166,7 +165,7 @@ Nice properties:
 
 Here's a figure showing the effect on the pmf of different shape parameters when $n=50$.
 
-```{code-cell} ipython3
+```{code-cell} python3
 def gen_probs(n, a, b):
     probs = np.zeros(n+1)
     for k in range(n+1):
@@ -189,7 +188,7 @@ plt.show()
 We will first create a class `CareerWorkerProblem` which will hold the
 default parameterizations of the model and an initial guess for the value function.
 
-```{code-cell} ipython3
+```{code-cell} python3
 class CareerWorkerProblem:
 
     def __init__(self,
@@ -222,7 +221,7 @@ the corresponding Bellman operator $T$ and the greedy policy function.
 In this model, $T$ is defined by $Tv(\theta, \epsilon) = \max\{I, II, III\}$, where
 $I$, $II$ and $III$ are as given in {eq}`eyes`.
 
-```{code-cell} ipython3
+```{code-cell} python3
 def operator_factory(cw, parallel_flag=True):
 
     """
@@ -278,7 +277,7 @@ def operator_factory(cw, parallel_flag=True):
 Lastly, `solve_model` will  take an instance of `CareerWorkerProblem` and
 iterate using the Bellman operator to find the fixed point of the Bellman equation.
 
-```{code-cell} ipython3
+```{code-cell} python3
 def solve_model(cw,
                 use_parallel=True,
                 tol=1e-4,
@@ -312,7 +311,7 @@ def solve_model(cw,
 
 Here's the solution to the model -- an approximate value function
 
-```{code-cell} ipython3
+```{code-cell} python3
 cw = CareerWorkerProblem()
 T, get_greedy = operator_factory(cw)
 v_star = solve_model(cw, verbose=False)
@@ -334,7 +333,7 @@ plt.show()
 
 And here is the optimal policy
 
-```{code-cell} ipython3
+```{code-cell} python3
 fig, ax = plt.subplots(figsize=(6, 6))
 tg, eg = np.meshgrid(cw.θ, cw.ϵ)
 lvls = (0.5, 1.5, 2.5, 3.5)
@@ -394,7 +393,7 @@ In reading the code, recall that `optimal_policy[i, j]` = policy at
 $(\theta_i, \epsilon_j)$ = either 1, 2 or 3; meaning 'stay put',
 'new job' and 'new life'.
 
-```{code-cell} ipython3
+```{code-cell} python3
 F = np.cumsum(cw.F_probs)
 G = np.cumsum(cw.G_probs)
 v_star = solve_model(cw, verbose=False)
@@ -466,7 +465,7 @@ Repeat the exercise with $\beta=0.99$ and interpret the change.
 
 The median for the original parameterization can be computed as follows
 
-```{code-cell} ipython3
+```{code-cell} python3
 cw = CareerWorkerProblem()
 F = np.cumsum(cw.F_probs)
 G = np.cumsum(cw.G_probs)
@@ -482,12 +481,12 @@ def passage_time(optimal_policy, F, G):
         if optimal_policy[i, j] == 1:    # Stay put
             return t
         elif optimal_policy[i, j] == 2:  # New job
-            j = qe.random.draw(G, size=1)[0]
+            j = qe.random.draw(G)
         else:                            # New life
-            i  = qe.random.draw(F, size=1)[0]
-            j  = qe.random.draw(G, size=1)[0]
+            i, j  = qe.random.draw(F), qe.random.draw(G)
         t += 1
-@njit
+
+@njit(parallel=True)
 def median_time(optimal_policy, F, G, M=25000):
     samples = np.empty(M)
     for i in prange(M):
@@ -522,7 +521,7 @@ figure -- interpret.
 
 Here is one solution
 
-```{code-cell} ipython3
+```{code-cell} python3
 cw = CareerWorkerProblem(G_a=100, G_b=100)
 T, get_greedy = operator_factory(cw)
 v_star = solve_model(cw, verbose=False)
diff --git a/lectures/mccall_model.md b/lectures/mccall_model.md
index c1a495f13..a3c02487e 100644
--- a/lectures/mccall_model.md
+++ b/lectures/mccall_model.md
@@ -3,10 +3,8 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.14.4
 kernelspec:
-  display_name: Python 3 (ipykernel)
+  display_name: Python 3
   language: python
   name: python3
 ---
@@ -36,9 +34,10 @@ and the pros and cons as they themselves see them." -- Robert E. Lucas, Jr.
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython3
-:tags: [hide-output]
-
+```{code-cell} ipython
+---
+tags: [hide-output]
+---
 !pip install quantecon
 ```
 
@@ -60,7 +59,7 @@ As we'll see, McCall's model is not only interesting in its own right but also a
 
 Let's start with some imports:
 
-```{code-cell} ipython3
+```{code-cell} ipython
 %matplotlib inline
 import matplotlib.pyplot as plt
 plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
@@ -344,21 +343,21 @@ $v$.
 Our default for $q$, the distribution of the state process, will be
 [Beta-binomial](https://en.wikipedia.org/wiki/Beta-binomial_distribution).
 
-```{code-cell} ipython3
+```{code-cell} python3
 n, a, b = 50, 200, 100                        # default parameters
 q_default = BetaBinomial(n, a, b).pdf()       # default choice of q
 ```
 
 Our default set of values for wages will be
 
-```{code-cell} ipython3
+```{code-cell} python3
 w_min, w_max = 10, 60
 w_default = np.linspace(w_min, w_max, n+1)
 ```
 
 Here's a plot of the probabilities of different wage outcomes:
 
-```{code-cell} ipython3
+```{code-cell} python3
 fig, ax = plt.subplots()
 ax.plot(w_default, q_default, '-o', label='$q(w(i))$')
 ax.set_xlabel('wages')
@@ -373,7 +372,7 @@ We are going to use Numba to accelerate our code.
 
 The following helps Numba by providing some type
 
-```{code-cell} ipython3
+```{code-cell} python3
 mccall_data = [
     ('c', float64),      # unemployment compensation
     ('β', float64),      # discount factor
@@ -388,7 +387,7 @@ given the current state and an arbitrary feasible action.
 
 Default parameter values are embedded in the class.
 
-```{code-cell} ipython3
+```{code-cell} python3
 @jitclass(mccall_data)
 class McCallModel:
 
@@ -418,7 +417,7 @@ We will start from guess $v$ given by $v(i) = w(i) / (1 - β)$, which is the val
 
 Here's a function to implement this:
 
-```{code-cell} ipython3
+```{code-cell} python3
 def plot_value_function_seq(mcm, ax, num_plots=6):
     """
     Plot a sequence of value functions.
@@ -443,7 +442,7 @@ def plot_value_function_seq(mcm, ax, num_plots=6):
 
 Now let's create an instance of `McCallModel` and watch iterations  $T^k v$ converge from below:
 
-```{code-cell} ipython3
+```{code-cell} python3
 mcm = McCallModel()
 
 fig, ax = plt.subplots()
@@ -462,7 +461,7 @@ the reservation wage.
 
 We'll be using JIT compilation via Numba to turbocharge our loops.
 
-```{code-cell} ipython3
+```{code-cell} python3
 @jit(nopython=True)
 def compute_reservation_wage(mcm,
                              max_iter=500,
@@ -495,7 +494,7 @@ def compute_reservation_wage(mcm,
 
 The next line computes the reservation wage at  default parameters
 
-```{code-cell} ipython3
+```{code-cell} python3
 compute_reservation_wage(mcm)
 ```
 
@@ -507,7 +506,7 @@ parameters.
 In particular, let's look at what happens when we change $\beta$ and
 $c$.
 
-```{code-cell} ipython3
+```{code-cell} python3
 grid_size = 25
 R = np.empty((grid_size, grid_size))
 
@@ -520,7 +519,7 @@ for i, c in enumerate(c_vals):
         R[i, j] = compute_reservation_wage(mcm)
 ```
 
-```{code-cell} ipython3
+```{code-cell} python3
 fig, ax = plt.subplots()
 
 cs1 = ax.contourf(c_vals, β_vals, R.T, alpha=0.75)
@@ -613,7 +612,7 @@ The big difference here, however, is that we're iterating on a scalar $h$, rathe
 
 Here's an implementation:
 
-```{code-cell} ipython3
+```{code-cell} python3
 @jit(nopython=True)
 def compute_reservation_wage_two(mcm,
                                  max_iter=500,
@@ -731,7 +730,7 @@ Once your code is working, investigate how the reservation wage changes with $c$
 
 Here's one solution
 
-```{code-cell} ipython3
+```{code-cell} python3
 cdf = np.cumsum(q_default)
 
 @jit(nopython=True)
@@ -741,7 +740,7 @@ def compute_stopping_time(w_bar, seed=1234):
     t = 1
     while True:
         # Generate a wage draw
-        w = w_default[qe.random.draw(cdf, size=1)]
+        w = w_default[qe.random.draw(cdf)]
         # Stop when the draw is above the reservation wage
         if w >= w_bar:
             stopping_time = t
@@ -783,7 +782,7 @@ plt.show()
 
 Here is one solution:
 
-```{code-cell} ipython3
+```{code-cell} python3
 mccall_data_continuous = [
     ('c', float64),          # unemployment compensation
     ('β', float64),          # discount factor
@@ -833,7 +832,7 @@ $\beta$.
 
 We will do this using a contour plot.
 
-```{code-cell} ipython3
+```{code-cell} python3
 grid_size = 25
 R = np.empty((grid_size, grid_size))
 
@@ -846,7 +845,7 @@ for i, c in enumerate(c_vals):
         R[i, j] = compute_reservation_wage_continuous(mcmc)
 ```
 
-```{code-cell} ipython3
+```{code-cell} python3
 fig, ax = plt.subplots()
 
 cs1 = ax.contourf(c_vals, β_vals, R.T, alpha=0.75)

From a56dd0d210256089b71fafad75f9d9217fab8538 Mon Sep 17 00:00:00 2001
From: Humphrey Yang <humphrey.yang@anu.edu.au>
Date: Wed, 31 May 2023 00:01:02 +1000
Subject: [PATCH 23/23] move exercises in mccall model

---
 lectures/mccall_model.md | 104 +++++++++++++++++++--------------------
 1 file changed, 50 insertions(+), 54 deletions(-)

diff --git a/lectures/mccall_model.md b/lectures/mccall_model.md
index a3c02487e..93f074215 100644
--- a/lectures/mccall_model.md
+++ b/lectures/mccall_model.md
@@ -661,6 +661,56 @@ Repeat a large number of times and take the average.
 Plot mean unemployment duration as a function of $c$ in `c_vals`.
 ```
 
+```{solution-start} mm_ex1
+:class: dropdown
+```
+
+Here's one solution
+
+```{code-cell} python3
+cdf = np.cumsum(q_default)
+
+@jit(nopython=True)
+def compute_stopping_time(w_bar, seed=1234):
+
+    np.random.seed(seed)
+    t = 1
+    while True:
+        # Generate a wage draw
+        w = w_default[qe.random.draw(cdf)]
+        # Stop when the draw is above the reservation wage
+        if w >= w_bar:
+            stopping_time = t
+            break
+        else:
+            t += 1
+    return stopping_time
+
+@jit(nopython=True)
+def compute_mean_stopping_time(w_bar, num_reps=100000):
+    obs = np.empty(num_reps)
+    for i in range(num_reps):
+        obs[i] = compute_stopping_time(w_bar, seed=i)
+    return obs.mean()
+
+c_vals = np.linspace(10, 40, 25)
+stop_times = np.empty_like(c_vals)
+for i, c in enumerate(c_vals):
+    mcm = McCallModel(c=c)
+    w_bar = compute_reservation_wage_two(mcm)
+    stop_times[i] = compute_mean_stopping_time(w_bar)
+
+fig, ax = plt.subplots()
+
+ax.plot(c_vals, stop_times, label="mean unemployment duration")
+ax.set(xlabel="unemployment compensation", ylabel="months")
+ax.legend()
+
+plt.show()
+```
+
+```{solution-end}
+```
 
 ```{exercise-start}
 :label: mm_ex2
@@ -722,60 +772,6 @@ Once your code is working, investigate how the reservation wage changes with $c$
 ```{exercise-end}
 ```
 
-## Solutions
-
-```{solution-start} mm_ex1
-:class: dropdown
-```
-
-Here's one solution
-
-```{code-cell} python3
-cdf = np.cumsum(q_default)
-
-@jit(nopython=True)
-def compute_stopping_time(w_bar, seed=1234):
-
-    np.random.seed(seed)
-    t = 1
-    while True:
-        # Generate a wage draw
-        w = w_default[qe.random.draw(cdf)]
-        # Stop when the draw is above the reservation wage
-        if w >= w_bar:
-            stopping_time = t
-            break
-        else:
-            t += 1
-    return stopping_time
-
-@jit(nopython=True)
-def compute_mean_stopping_time(w_bar, num_reps=100000):
-    obs = np.empty(num_reps)
-    for i in range(num_reps):
-        obs[i] = compute_stopping_time(w_bar, seed=i)
-    return obs.mean()
-
-c_vals = np.linspace(10, 40, 25)
-stop_times = np.empty_like(c_vals)
-for i, c in enumerate(c_vals):
-    mcm = McCallModel(c=c)
-    w_bar = compute_reservation_wage_two(mcm)
-    stop_times[i] = compute_mean_stopping_time(w_bar)
-
-fig, ax = plt.subplots()
-
-ax.plot(c_vals, stop_times, label="mean unemployment duration")
-ax.set(xlabel="unemployment compensation", ylabel="months")
-ax.legend()
-
-plt.show()
-```
-
-```{solution-end}
-```
-
-
 ```{solution-start} mm_ex2
 :class: dropdown
 ```