using scipy.optimize to replace brent_max test

Author: Jiarui-ZH

lectures/ak2.md

@@ -15,7 +15,7 @@ kernelspec:
 
 In addition to what’s in Anaconda, this lecture will need the following libraries:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-output]
 
 !pip install --upgrade quantecon
@@ -407,15 +407,16 @@ $$
 
 ### Implementation
 
-```{code-cell}
+```{code-cell} ipython3
 import numpy as np
 import matplotlib.pyplot as plt
 from quantecon.optimize import brent_max
+from scipy.optimize import minimize_scalar
 ```
 
 For parameters $\alpha = 0.3$ and $\beta = 0.5$, let's compute  $\hat{K}$:
 
-```{code-cell}
+```{code-cell} ipython3
 # parameters
 α = 0.3
 β = 0.5
@@ -433,8 +434,7 @@ Knowing $\hat K$, we can calculate other equilibrium objects.
 
 Let's first define  some Python helper functions.
 
-```{code-cell}
-
+```{code-cell} ipython3
 def K_to_Y(K, α):
 
     return K ** α
@@ -465,28 +465,28 @@ def K_to_C(K, D, τ, r, α, β):
 
 We can use these helper functions to obtain steady state values $\hat{Y}$, $\hat{r}$, and $\hat{W}$ associated with  steady state values $\hat{K}$ and $\hat{r}$.
 
-```{code-cell}
+```{code-cell} ipython3
 Y_hat, r_hat, W_hat = K_to_Y(K_hat, α), K_to_r(K_hat, α), K_to_W(K_hat, α)
 Y_hat, r_hat, W_hat
 ```
 
 Since  steady state government debt $\hat{D}$ is  $0$, all taxes are  used to pay for government expenditures
 
-```{code-cell}
+```{code-cell} ipython3
 G_hat = τ_hat * Y_hat
 G_hat
 ```
 
 We use the optimal consumption plans to find  steady state consumptions for  young and  old
 
-```{code-cell}
+```{code-cell} ipython3
 Cy_hat, Co_hat = K_to_C(K_hat, D_hat, τ_hat, r_hat, α, β)
 Cy_hat, Co_hat
 ```
 
 Let's store the steady state quantities and prices using an array called `init_ss`
 
-```{code-cell}
+```{code-cell} ipython3
 init_ss = np.array([K_hat, Y_hat, Cy_hat, Co_hat,     # quantities
                     W_hat, r_hat,                     # prices
                     τ_hat, D_hat, G_hat               # policies
@@ -536,7 +536,7 @@ We'll then compute the single remaining undetermined policy variable from the go
 
 When we simulate  transition paths, it is useful to distinguish  **state variables** at time $t$  such as $K_t, Y_t, D_t, W_t, r_t$ from  **control variables** that include $C_{yt}, C_{ot}, \tau_{t}, G_t$.
 
-```{code-cell}
+```{code-cell} ipython3
 class ClosedFormTrans:
     """
     This class simulates length T transitional path of a economy
@@ -681,7 +681,7 @@ class ClosedFormTrans:
 
 We can create an instance `closed` for model parameters $\{\alpha, \beta\}$ and use it for various fiscal policy experiments.
 
-```{code-cell}
+```{code-cell} ipython3
 closed = ClosedFormTrans(α, β)
 ```
 
@@ -711,7 +711,7 @@ The first step is to prepare sequences of policy variables that describe  fiscal
 
 We must define  sequences of government expenditure $\{G_t\}_{t=0}^{T}$ and debt level $\{D_t\}_{t=0}^{T+1}$ in advance, then pass them  to the solver.
 
-```{code-cell}
+```{code-cell} ipython3
 T = 20
 
 # tax cut
@@ -730,7 +730,7 @@ Let's use the `simulate` method of `closed` to compute dynamic transitions.
 
 Note that we leave `τ_pol` as `None`, since the tax rates need to be determined to satisfy the government budget constraint.
 
-```{code-cell}
+```{code-cell} ipython3
 quant_seq1, price_seq1, policy_seq1 = closed.simulate(T, init_ss,
                                                       D_pol=D_seq,
                                                       G_pol=G_seq)
@@ -739,7 +739,7 @@ closed.plot()
 
 We can also  experiment with a lower tax cut rate, such as $0.2$.
 
-```{code-cell}
+```{code-cell} ipython3
 # lower tax cut rate
 τ0 = 0.15 * (1 - 0.2)
 
@@ -753,7 +753,7 @@ quant_seq2, price_seq2, policy_seq2 = closed.simulate(T, init_ss,
                                                       G_pol=G_seq)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig, axs = plt.subplots(3, 3, figsize=(14, 10))
 
 # quantities
@@ -797,7 +797,7 @@ The government targets  the same tax rate $\tau_t=\hat{\tau}$ and to accumulate
 
 To conduct  this experiment, we pass `τ_seq` and `G_seq` as inputs  and let `D_pol`  be determined along the path by satisfying the government budget constraint.
 
-```{code-cell}
+```{code-cell} ipython3
 # government expenditure cut by a half
 G_seq = τ_hat * 0.5 * Y_hat * np.ones(T+1)
 
@@ -812,7 +812,7 @@ As the government accumulates the asset and uses it in production, the  rental r
 
 As a result,  the ratio  $-\frac{D_t}{K_t}$ of the  government asset to  physical capital used in production will increase over time
 
-```{code-cell}
+```{code-cell} ipython3
 plt.plot(range(T+1), -closed.policy_seq[:-1, 1] / closed.quant_seq[:, 0])
 plt.xlabel('t')
 plt.title('-D/K');
@@ -839,7 +839,7 @@ But now let  the government cut its  expenditures only  at $t=0$.
 
 From $t \geq 1$, the government expeditures  return to  $\hat{G}$  and  $\tau_t$ adjusts to maintain the   asset level $-D_t = -D_1$.
 
-```{code-cell}
+```{code-cell} ipython3
 # sequence of government purchase
 G_seq = τ_hat * Y_hat * np.ones(T+1)
 G_seq[0] = 0
@@ -911,7 +911,7 @@ Let's implement this "guess and verify" approach
 
 We start by defining the Cobb-Douglas utility function
 
-```{code-cell}
+```{code-cell} ipython3
 def U(Cy, Co, β):
 
     return (Cy ** β) * (Co ** (1-β))
@@ -921,7 +921,7 @@ We use `Cy_val` to compute the lifetime value of an arbitrary consumption plan,
 
 Note that it requires knowing future prices $r_{t+1}$ and tax rate $\tau_{t+1}$.
 
-```{code-cell}
+```{code-cell} ipython3
 def Cy_val(Cy, W, r_next, τ, τ_next, δy, δo_next, β):
 
     # Co given by the budget constraint
@@ -934,14 +934,16 @@ An optimal consumption plan $C_y^*$ can be found by maximizing `Cy_val`.
 
 Here is an example that computes optimal consumption $C_y^*=\hat{C}_y$ in the steady state  with $\delta_{yt}=\delta_{ot}=0,$ like one that we studied earlier
 
-```{code-cell}
+```{code-cell} ipython3
 W, r_next, τ, τ_next = W_hat, r_hat, τ_hat, τ_hat
 δy, δo_next = 0, 0
 
-Cy_opt, U_opt, _ = brent_max(Cy_val,            # maximand
-                             1e-6,              # lower bound
-                             W*(1-τ)-δy-1e-6,   # upper bound
-                             args=(W, r_next, τ, τ_next, δy, δo_next, β))
+result = minimize_scalar(lambda Cy: -Cy_val(Cy, W, r_next, τ, τ_next, δy, δo_next, β), 
+                         bounds=(1e-6, W*(1-τ)-δy-1e-6), 
+                         method='bounded')
+
+Cy_opt = result.x
+U_opt = -result.fun
 
 Cy_opt, U_opt
 ```
@@ -950,7 +952,7 @@ Let's define a Python class `AK2` that  computes the transition paths  with the
 
 It can handle   nonzero lump sum taxes
 
-```{code-cell}
+```{code-cell} ipython3
 class AK2():
     """
     This class simulates length T transitional path of a economy
@@ -1121,13 +1123,13 @@ class AK2():
 
 We can initialize an instance of class `AK2` with model parameters $\{\alpha, \beta\}$ and then use it to conduct fiscal policy experiments.
 
-```{code-cell}
+```{code-cell} ipython3
 ak2 = AK2(α, β)
 ```
 
 We first examine that the "guess and verify" method leads to the same numerical results as we obtain with the closed form solution when lump sum taxes are muted
 
-```{code-cell}
+```{code-cell} ipython3
 δy_seq = np.ones(T+2) * 0.
 δo_seq = np.ones(T+2) * 0.
 
@@ -1141,22 +1143,22 @@ D_pol[0] = D_hat
 D_pol[1:] = D1
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 quant_seq3, price_seq3, policy_seq3 = ak2.simulate(T, init_ss,
                                                    δy_seq, δo_seq,
                                                    D_pol=D_pol, G_pol=G_pol,
                                                    verbose=True)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 ak2.plot()
 ```
 
 Next, we  activate  lump sum taxes. 
 
 Let's alter  our  {ref}`exp-tax-cut`  fiscal policy experiment by assuming that  the government also increases  lump sum taxes for both  young and old  people $\delta_{yt}=\delta_{ot}=0.005, t\geq0$.
 
-```{code-cell}
+```{code-cell} ipython3
 δy_seq = np.ones(T+2) * 0.005
 δo_seq = np.ones(T+2) * 0.005
 
@@ -1170,7 +1172,7 @@ quant_seq4, price_seq4, policy_seq4 = ak2.simulate(T, init_ss,
 
 Note how   "crowding out"  has been  mitigated.
 
-```{code-cell}
+```{code-cell} ipython3
 fig, axs = plt.subplots(3, 3, figsize=(14, 10))
 
 # quantities
@@ -1224,7 +1226,7 @@ We can  use our code to compute the transition ignited by  launching this system
 
 Let's compare the results to the {ref}`exp-tax-cut`.
 
-```{code-cell}
+```{code-cell} ipython3
 δy_seq = np.ones(T+2) * Cy_hat * 0.1
 δo_seq = np.ones(T+2) * -Cy_hat * 0.1
 
@@ -1235,7 +1237,7 @@ quant_seq5, price_seq5, policy_seq5 = ak2.simulate(T, init_ss,
                                                    D_pol=D_pol, G_pol=G_pol)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig, axs = plt.subplots(3, 3, figsize=(14, 10))
 
 # quantities