Replace brent_max with minimize_scalar.

Author: Jiarui-ZH

lectures/ak2.md

@@ -551,34 +551,35 @@ class ClosedFormTrans:
         self.α, self.β = α, β
 
     def simulate(self,
-                T,           # length of transitional path to simulate
-                init_ss,     # initial steady state
-                τ_pol=None,  # sequence of tax rates
-                D_pol=None,  # sequence of government debt levels
-                G_pol=None): # sequence of government purchases
-
+                 T,           # length of transitional path to simulate
+                 init_ss,     # initial steady state
+                 τ_pol=None,  # sequence of tax rates
+                 D_pol=None,  # sequence of government debt levels
+                 G_pol=None): # sequence of government purchases
+    
         α, β = self.α, self.β
-
+    
         # unpack the steady state variables
         K_hat, Y_hat, Cy_hat, Co_hat = init_ss[:4]
         W_hat, r_hat = init_ss[4:6]
         τ_hat, D_hat, G_hat = init_ss[6:9]
-
+    
         # initialize array containers
         # K, Y, Cy, Co
         quant_seq = np.empty((T+1, 4))
-
+    
         # W, r
         price_seq = np.empty((T+1, 2))
-
+    
         # τ, D, G
         policy_seq = np.empty((T+2, 3))
-
+    
         # t=0, starting from steady state
         K0, Y0 = K_hat, Y_hat
         W0, r0 = W_hat, r_hat
         D0 = D_hat
-
+        δo = 0  # Define δo (set to 0 for this example, adjust as needed)
+    
         # fiscal policy
         if τ_pol is None:
             D1 = D_pol[1]
@@ -592,27 +593,27 @@ class ClosedFormTrans:
             D1 = D_pol[1]
             τ0 = τ_pol[0]
             G0 = τ0 * (Y0 + r0 * D0) + D1 - (1 + r0) * D0
-
+    
         # optimal consumption plans
         Cy0, Co0 = K_to_C(K0, D0, τ0, r0, α, β)
-
+    
         # t=0 economy
         quant_seq[0, :] = K0, Y0, Cy0, Co0
         price_seq[0, :] = W0, r0
         policy_seq[0, :] = τ0, D0, G0
         policy_seq[1, 1] = D1
-
+    
         # starting from t=1 to T
         for t in range(1, T+1):
-
+    
             # transition of K
             K_old, τ_old = quant_seq[t-1, 0], policy_seq[t-1, 0]
             D = policy_seq[t, 1]
             K = K_old ** α * (1 - τ_old) * (1 - α) * (1 - β) - D
-
+    
             # output, capital return, wage
             Y, r, W = K_to_Y(K, α), K_to_r(K, α), K_to_W(K, α)
-
+    
             # to satisfy the government budget constraint
             if τ_pol is None:
                 D = D_pol[t]
@@ -629,22 +630,32 @@ class ClosedFormTrans:
                 D_next = D_pol[t+1]
                 τ = τ_pol[t]
                 G = τ * (Y + r * D) + D_next - (1 + r) * D
-
+    
             # optimal consumption plans
-            Cy, Co = K_to_C(K, D, τ, r, α, β)
-
+            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 = Cy_opt
+            Co = (1 + r * (1 - τ)) * (K + D) - δo
+    
             # store time t economy aggregates
             quant_seq[t, :] = K, Y, Cy, Co
             price_seq[t, :] = W, r
             policy_seq[t, 0] = τ
             policy_seq[t+1, 1] = D_next
             policy_seq[t, 2] = G
-
+    
         self.quant_seq = quant_seq
         self.price_seq = price_seq
         self.policy_seq = policy_seq
-
+    
         return quant_seq, price_seq, policy_seq
+    
+    
 
     def plot(self):
 
@@ -967,57 +978,57 @@ class AK2():
         self.α, self.β = α, β
 
     def simulate(self,
-                T,           # length of transitional path to simulate
-                init_ss,     # initial steady state
-                δy_seq,      # sequence of lump sum tax for the young
-                δo_seq,      # sequence of lump sum tax for the old
-                τ_pol=None,  # sequence of tax rates
-                D_pol=None,  # sequence of government debt levels
-                G_pol=None,  # sequence of government purchases
-                verbose=False,
-                max_iter=500,
-                tol=1e-5):
-
+                 T,           # length of transitional path to simulate
+                 init_ss,     # initial steady state
+                 δy_seq,      # sequence of lump sum tax for the young
+                 δo_seq,      # sequence of lump sum tax for the old
+                 τ_pol=None,  # sequence of tax rates
+                 D_pol=None,  # sequence of government debt levels
+                 G_pol=None,  # sequence of government purchases
+                 verbose=False,
+                 max_iter=500,
+                 tol=1e-5):
+    
         α, β = self.α, self.β
-
+    
         # unpack the steady state variables
         K_hat, Y_hat, Cy_hat, Co_hat = init_ss[:4]
         W_hat, r_hat = init_ss[4:6]
         τ_hat, D_hat, G_hat = init_ss[6:9]
-
+    
         # K, Y, Cy, Co
         quant_seq = np.empty((T+2, 4))
-
+    
         # W, r
         price_seq = np.empty((T+2, 2))
-
+    
         # τ, D, G
         policy_seq = np.empty((T+2, 3))
         policy_seq[:, 1] = D_pol
         policy_seq[:, 2] = G_pol
-
+    
         # initial guesses of prices
         price_seq[:, 0] = np.ones(T+2) * W_hat
         price_seq[:, 1] = np.ones(T+2) * r_hat
-
+    
         # initial guesses of policies
         policy_seq[:, 0] = np.ones(T+2) * τ_hat
-
+    
         # t=0, starting from steady state
         quant_seq[0, :2] = K_hat, Y_hat
-
+    
         if verbose:
             # prepare to plot iterations until convergence
             fig, axs = plt.subplots(1, 3, figsize=(14, 4))
-
+    
         # containers for checking convergence
         price_seq_old = np.empty_like(price_seq)
         policy_seq_old = np.empty_like(policy_seq)
-
+    
         # start iteration
         i_iter = 0
         while True:
-
+    
             if verbose:
                 # plot current prices at ith iteration
                 for i, name in enumerate(['W', 'r']):
@@ -1029,11 +1040,11 @@ class AK2():
                 axs[2].legend(bbox_to_anchor=(1.05, 1), loc='upper left')
                 axs[2].set_title('τ')
                 axs[2].set_xlabel('t')
-
+    
             # store old prices from last iteration
             price_seq_old[:] = price_seq
             policy_seq_old[:] = policy_seq
-
+    
             # start updating quantities and prices
             for t in range(T+1):
                 K, Y = quant_seq[t, :2]
@@ -1043,40 +1054,41 @@ class AK2():
                 τ_next, D_next, G_next = policy_seq[t+1, :]
                 δy, δo = δy_seq[t], δo_seq[t]
                 δy_next, δo_next = δy_seq[t+1], δo_seq[t+1]
-
+    
                 # consumption for the old
                 Co = (1 + r * (1 - τ)) * (K + D) - δo
-
+    
                 # optimal consumption for the young
-                out = brent_max(Cy_val, 1e-6, W*(1-τ)-δy-1e-6,
-                                args=(W, r_next, τ, τ_next,
-                                      δy, δo_next, β))
-                Cy = out[0]
-
+                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 = result.x
+    
                 quant_seq[t, 2:] = Cy, Co
                 τ_num = ((1 + r) * D + G - D_next - δy - δo)
                 τ_denom = (Y + r * D)
                 policy_seq[t, 0] = τ_num / τ_denom
-
+    
                 # saving of the young
                 A_next = W * (1 - τ) - δy - Cy
-
+    
                 # transition of K
                 K_next = A_next - D_next
                 Y_next = K_to_Y(K_next, α)
                 W_next, r_next = K_to_W(K_next, α), K_to_r(K_next, α)
-
+    
                 quant_seq[t+1, :2] = K_next, Y_next
                 price_seq[t+1, :] = W_next, r_next
-
+    
             i_iter += 1
-
+    
             if (np.max(np.abs(price_seq_old - price_seq)) < tol) & \
                (np.max(np.abs(policy_seq_old - policy_seq)) < tol):
                 if verbose:
                     print(f"Converge using {i_iter} iterations")
                 break
-
+    
             if i_iter > max_iter:
                 if verbose:
                     print(f"Fail to converge using {i_iter} iterations")
@@ -1085,9 +1097,10 @@ class AK2():
         self.quant_seq = quant_seq
         self.price_seq = price_seq
         self.policy_seq = policy_seq
-
+    
         return quant_seq, price_seq, policy_seq
 
+
     def plot(self):
 
         quant_seq = self.quant_seq