Skip to content

Commit

Permalink
Replace brent_max with minimize_scalar.
Browse files Browse the repository at this point in the history
  • Loading branch information
@Jiarui-ZH
Jiarui-ZH committed Aug 7, 2024
1 parent 16413c8 commit cdef435
Showing 1 changed file with 74 additions and 61 deletions.
135 changes: 74 additions & 61 deletions lectures/ak2.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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]
Expand All @@ -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):
Expand Down Expand Up @@ -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']):
Expand All @@ -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]
Expand All @@ -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")
Expand All @@ -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
Expand Down

0 comments on commit cdef435

Please sign in to comment.