From c3d4ac6251f81dbf97e0ecbe63bbbd091c3fbeaf Mon Sep 17 00:00:00 2001 From: jdebacker Date: Fri, 6 Oct 2017 07:53:16 -0400 Subject: [PATCH 1/3] add fin frictions --- Code/SS.py | 26 ++++++++++++++++---------- Code/VFI.py | 20 ++++++++++++-------- Code/execute.py | 25 ++++++++++++++++++------- 3 files changed, 46 insertions(+), 25 deletions(-) diff --git a/Code/SS.py b/Code/SS.py index 6c4f9ba..3f2f942 100644 --- a/Code/SS.py +++ b/Code/SS.py @@ -16,7 +16,7 @@ @numba.jit -def find_SD(PF, Pi, sizez, sizek): +def find_SD(PF, Pi, sizez, sizek, Gamma_initial): ''' ------------------------------------------------------------------------ Compute the stationary distribution of firms over (k, z) @@ -29,7 +29,7 @@ def find_SD(PF, Pi, sizez, sizek): HGamma = operated on stationary distribution ------------------------------------------------------------------------ ''' - Gamma = np.ones((sizez, sizek)) * (1 / (sizek * sizez)) + Gamma = Gamma_initial SDtol = 1e-12 SDdist = 7 SDiter = 0 @@ -60,14 +60,16 @@ def find_SD(PF, Pi, sizez, sizek): def GE_loop(w, *args): - (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, sizek, sizez, h, - tax_params) = args - op, e, l_d, y = VFI.get_firmobjects(w, z, K, alpha_k, alpha_l, - delta, psi, sizez, sizek, - tax_params) - VF, PF, optK, optI = VFI.VFI(e, betafirm, delta, K, Pi, sizez, sizek, - tax_params) - Gamma = find_SD(PF, Pi, sizez, sizek) + global VF_initial, Gamma_initial + (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, sizek, + sizez, h, tax_params, VF_initial, Gamma_initial) = args + # print('VF_initial = ', VF_initial[:4, 50:60]) + op, e, l_d, y, eta = VFI.get_firmobjects(w, z, K, alpha_k, alpha_l, + delta, psi, eta0, eta1, + sizez, sizek, tax_params) + VF, PF, optK, optI = VFI.VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, + tax_params, VF_initial) + Gamma = find_SD(PF, Pi, sizez, sizek, Gamma_initial) L_d = (Gamma * l_d).sum() Y = (Gamma * y).sum() I = (Gamma * optI).sum() @@ -76,6 +78,10 @@ def GE_loop(w, *args): L_s = get_L_s(w, C, h) print('Labor demand and supply = ', L_d, L_s) MCdist = L_d - L_s + VF_initial = VF + Gamma_initial = Gamma + # print('VF_initial = ', VF_initial[:4, 50:60]) + return MCdist diff --git a/Code/VFI.py b/Code/VFI.py index d48aad5..1d10b37 100644 --- a/Code/VFI.py +++ b/Code/VFI.py @@ -16,7 +16,7 @@ @numba.jit -def create_Vmat(EV, e, betafirm, Pi, sizez, sizek, Vmat, tax_params): +def create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, Vmat, tax_params): ''' ------------------------------------------------------------------------ This function loops over the state and control variables, operating on the @@ -50,8 +50,9 @@ def create_Vmat(EV, e, betafirm, Pi, sizez, sizek, Vmat, tax_params): for i in range(sizez): # loop over z for j in range(sizek): # loop over k for m in range(sizek): # loop over k' - Vmat[i, j, m] = (((1 - tau_d) / (1 - tau_g)) * e[i, j, m] - + betafirm * EV[i, m]) + Vmat[i, j, m] = (((1 - tau_d) / (1 - tau_g)) * + (e[i, j, m] - eta[i, j, m]) + + betafirm * EV[i, m]) return Vmat @@ -73,7 +74,8 @@ def adj_costs(kprime, k, delta, psi): @numba.jit -def get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, sizez, sizek, tax_params): +def get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, eta0, eta1, + sizez, sizek, tax_params): ''' ------------------------------------------------------------------------- Generating possible cash flow levels @@ -112,10 +114,12 @@ def get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, sizez, sizek, tax_par * K[j]) - adj_costs(K[m], K[j], delta, psi)) - return op, e, l_d, y + eta = (eta0 + eta1 * e) * (e < 0) + return op, e, l_d, y, eta -def VFI(e, betafirm, delta, K, Pi, sizez, sizek, tax_params): + +def VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, tax_params, VF_initial): ''' ------------------------------------------------------------------------ Value Function Iteration @@ -142,14 +146,14 @@ def VFI(e, betafirm, delta, K, Pi, sizez, sizek, tax_params): VFtol = 1e-6 VFdist = 7.0 VFmaxiter = 3000 - V = np.zeros((sizez, sizek)) # initial guess at value function + V = VF_initial Vmat = np.empty((sizez, sizek, sizek)) # initialize Vmat matrix Vstore = np.empty((sizez, sizek, VFmaxiter)) # initialize Vstore array VFiter = 1 while VFdist > VFtol and VFiter < VFmaxiter: TV = V EV = np.dot(Pi, V) # expected VF (expectation over z') - Vmat = create_Vmat(EV, e, betafirm, Pi, sizez, sizek, Vmat, tax_params) + Vmat = create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, Vmat, tax_params) Vstore[:, :, VFiter] = V.reshape(sizez, sizek) # store value function # at each iteration for graphing later diff --git a/Code/execute.py b/Code/execute.py index 5be318c..e85fe69 100644 --- a/Code/execute.py +++ b/Code/execute.py @@ -1,6 +1,7 @@ # TODO: (0) check equations. Also, can I replicate GM (2010)??(1) add costly equity, (2) add debt, (3) add more tax params - tax deprec rate, interest deduct, etc (as outline in OG-USA guide), # (4) think more about tables and figures, (5) update figures scripot (6) add tables script, # (7) write moments script to generate moments (maybe do before tables), (8) estimation script +#(8) work on passing past solution of VFI so have good starting values ''' ------------------------------------------------------------------------ @@ -23,6 +24,7 @@ import scipy.optimize as opt import os import time +import numpy as np import grids import VFI @@ -71,6 +73,10 @@ rho = 0.7605 sigma_eps = 0.213 +# financial frictions +eta0 = 0.04 # fixed cost to issuing equity +eta1 = 0.02 # linear cost to issuing equity + # taxes tau_i = 0.25 tau_d = 0.25 @@ -133,9 +139,14 @@ ------------------------------------------------------------------------ ''' start_time = time.clock() +VF_initial = np.zeros((sizez, sizek)) # initial guess at Value Function +# initial guess at stationary distribution +Gamma_initial = np.ones((sizez, sizek)) * (1 / (sizek * sizez)) results = opt.bisect(SS.GE_loop, 0.1, 2, args=(alpha_k, alpha_l, delta, psi, - betafirm, K, z, Pi, sizek, - sizez, h, tax_params), + betafirm, K, z, Pi, eta0, + eta1, sizek, sizez, h, + tax_params, VF_initial, + Gamma_initial), xtol=1e-4, full_output=True) print('SS results: ', results) w = results[0] @@ -148,11 +159,11 @@ Find model outputs given eq'm wage rate ------------------------------------------------------------------------ ''' -op, e, l_d, y = VFI.get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, - sizez, sizek, tax_params) -VF, PF, optK, optI = VFI.VFI(e, betafirm, delta, K, Pi, sizez, sizek, - tax_params) -Gamma = SS.find_SD(PF, Pi, sizez, sizek) +op, e, l_d, y, eta = VFI.get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, + eta0, eta1, sizez, sizek, tax_params) +VF, PF, optK, optI = VFI.VFI(e, eta, betafirm, delta, K, Pi, + sizez, sizek, tax_params, VF_initial) +Gamma = SS.find_SD(PF, Pi, sizez, sizek, Gamma_initial) print('Sum of Gamma = ', Gamma.sum()) ''' From c3804f4690f3764329cb18e4c42cedff3f1ce99b Mon Sep 17 00:00:00 2001 From: jdebacker Date: Fri, 6 Oct 2017 15:52:27 -0400 Subject: [PATCH 2/3] add custom golden rule search for finding GE --- Code/SS.py | 113 +++++++++++++++++++++++++++++++++++++++++++++--- Code/VFI.py | 8 ++-- Code/execute.py | 19 ++++++-- 3 files changed, 127 insertions(+), 13 deletions(-) diff --git a/Code/SS.py b/Code/SS.py index 3f2f942..65ab189 100644 --- a/Code/SS.py +++ b/Code/SS.py @@ -45,11 +45,11 @@ def find_SD(PF, Pi, sizez, sizek, Gamma_initial): Gamma = HGamma SDiter += 1 - if SDiter < SDmaxiter: - print('Stationary distribution converged after this many iterations: ', - SDiter) - else: - print('Stationary distribution did not converge') + # if SDiter < SDmaxiter: + # # print('Stationary distribution converged after this many iterations: ', + # SDiter) + # else: + # print('Stationary distribution did not converge') # Check if state space is binding if Gamma.sum(axis=0)[-1] > 0.002: @@ -76,7 +76,7 @@ def GE_loop(w, *args): Psi = (Gamma * VFI.adj_costs(optK, K, delta, psi)).sum() C = Y - I - Psi L_s = get_L_s(w, C, h) - print('Labor demand and supply = ', L_d, L_s) + # print('Labor demand and supply = ', L_d, L_s) MCdist = L_d - L_s VF_initial = VF Gamma_initial = Gamma @@ -86,7 +86,108 @@ def GE_loop(w, *args): return MCdist +def Market_Clearing(w, args): + (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, sizek, + sizez, h, tax_params, VF_initial, Gamma_initial) = args + # print('VF_initial = ', VF_initial[:4, 50:60]) + op, e, l_d, y, eta = VFI.get_firmobjects(w, z, K, alpha_k, alpha_l, + delta, psi, eta0, eta1, + sizez, sizek, tax_params) + VF, PF, optK, optI = VFI.VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, + tax_params, VF_initial) + Gamma = find_SD(PF, Pi, sizez, sizek, Gamma_initial) + L_d = (Gamma * l_d).sum() + Y = (Gamma * y).sum() + I = (Gamma * optI).sum() + Psi = (Gamma * VFI.adj_costs(optK, K, delta, psi)).sum() + C = Y - I - Psi + L_s = get_L_s(w, C, h) + # print('Labor demand and supply = ', L_d, L_s) + MCdist = np.absolute(L_d - L_s) + + return MCdist, VF, Gamma + + def get_L_s(w, C, h): L_s = w / (h * C) return L_s + + +def golden_ratio_eqm(lb, ub, args, tolerance=1e-4): + ''' + Use the golden section search method to find the GE + ''' + (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, sizek, + sizez, h, tax_params, VF_initial, Gamma_initial) = args + golden_ratio = 2 / (np.sqrt(5) + 1) + + # Use the golden ratio to set the initial test points + x1 = ub - golden_ratio * (ub - lb) + x2 = lb + golden_ratio * (ub - lb) + + # Evaluate the function at the test points + f1, VF1, Gamma1 = Market_Clearing(x1, args) + f2, VF2, Gamma2 = Market_Clearing(x2, args) + + iteration = 0 + while (np.absolute(ub - lb) > tolerance): + iteration += 1 + # print('Iteration #', iteration) + # print('f1 =', f1) + # print('f2 =', f2) + + if (f2 > f1): + # then the minimum is to the left of x2 + # let x2 be the new upper bound + # let x1 be the new upper test point + # print('f2 > f1') + # Set the new upper bound + ub = x2 + # print('New Upper Bound =', ub) + # print('New Lower Bound =', lb) + # Set the new upper test point + # Use the special result of the golden ratio + x2 = x1 + # print('New Upper Test Point = ', x2) + f2 = f1 + + # Set the new lower test point + x1 = ub - golden_ratio * (ub - lb) + # print('New Lower Test Point = ', x1) + args = (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, + eta0, eta1, sizek, sizez, h, tax_params, VF1, Gamma1) + f1, VF1, Gamma1 = Market_Clearing(x1, args) + else: + # print('f2 < f1') + # the minimum is to the right of x1 + # let x1 be the new lower bound + # let x2 be the new lower test point + + # Set the new lower bound + lb = x1 + # print('New Upper Bound =', ub) + # print('New Lower Bound =', lb) + + # Set the new lower test point + x1 = x2 + # print('New Lower Test Point = ', x1) + + f1 = f2 + + # Set the new upper test point + x2 = lb + golden_ratio * (ub - lb) + # print('New Upper Test Point = ', x2) + args = (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, + eta0, eta1, sizek, sizez, h, tax_params, VF2, Gamma2) + f2, VF2, Gamma2 = Market_Clearing(x2, args) + + # Use the mid-point of the final interval as the estimate of the optimzer + # print('', '\n') + # print('Final Lower Bound =', lb, '\n') + # print('Final Upper Bound =', ub, '\n') + est_min = (lb + ub)/2 + # print('Estimated Minimizer =', est_min, '\n') + # print('MC_dist =', f2, '\n') + + return est_min diff --git a/Code/VFI.py b/Code/VFI.py index 1d10b37..9c96e13 100644 --- a/Code/VFI.py +++ b/Code/VFI.py @@ -164,10 +164,10 @@ def VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, tax_params, VF_initial): # print('VF iteration: ', VFiter) VFiter += 1 - if VFiter < VFmaxiter: - print('Value function converged after this many iterations:', VFiter) - else: - print('Value function did not converge') + # if VFiter < VFmaxiter: + # print('Value function converged after this many iterations:', VFiter) + # else: + # print('Value function did not converge') VF = V # solution to the functional equation diff --git a/Code/execute.py b/Code/execute.py index e85fe69..5d0fced 100644 --- a/Code/execute.py +++ b/Code/execute.py @@ -138,21 +138,34 @@ Solve for general equilibrium ------------------------------------------------------------------------ ''' -start_time = time.clock() VF_initial = np.zeros((sizez, sizek)) # initial guess at Value Function # initial guess at stationary distribution Gamma_initial = np.ones((sizez, sizek)) * (1 / (sizek * sizez)) +start_time = time.time() results = opt.bisect(SS.GE_loop, 0.1, 2, args=(alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, sizek, sizez, h, tax_params, VF_initial, Gamma_initial), xtol=1e-4, full_output=True) +end_time = time.time() +print('Bisection time = ', end_time - start_time) print('SS results: ', results) w = results[0] -GE_time = time.clock() - start_time -print('Solving the GE model took ', GE_time, ' seconds to solve') +# GE_time = time.clock() - start_time +# print('Solving the GE model took ', GE_time, ' seconds to solve') print('SS wage rate = ', w) +gr_args = (alpha_k, alpha_l, delta, psi, + betafirm, K, z, Pi, eta0, eta1, + sizek, sizez, h, tax_params, + VF_initial, Gamma_initial) +start_time = time.time() +w = SS.golden_ratio_eqm(0.1, 2, gr_args, tolerance=1e-4) +end_time = time.time() +print('Golden Rule time = ', end_time - start_time) +print('SS wage rate: ', w) +quit() + ''' ------------------------------------------------------------------------ From a0a1cb8c2d24f2a9ce635ae50af7d288030f3d62 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Fri, 6 Oct 2017 15:55:32 -0400 Subject: [PATCH 3/3] clean up old root finder --- Code/SS.py | 27 --------------------------- Code/execute.py | 22 +++------------------- 2 files changed, 3 insertions(+), 46 deletions(-) diff --git a/Code/SS.py b/Code/SS.py index 65ab189..bdab6b4 100644 --- a/Code/SS.py +++ b/Code/SS.py @@ -59,33 +59,6 @@ def find_SD(PF, Pi, sizez, sizek, Gamma_initial): return Gamma -def GE_loop(w, *args): - global VF_initial, Gamma_initial - (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, sizek, - sizez, h, tax_params, VF_initial, Gamma_initial) = args - # print('VF_initial = ', VF_initial[:4, 50:60]) - op, e, l_d, y, eta = VFI.get_firmobjects(w, z, K, alpha_k, alpha_l, - delta, psi, eta0, eta1, - sizez, sizek, tax_params) - VF, PF, optK, optI = VFI.VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, - tax_params, VF_initial) - Gamma = find_SD(PF, Pi, sizez, sizek, Gamma_initial) - L_d = (Gamma * l_d).sum() - Y = (Gamma * y).sum() - I = (Gamma * optI).sum() - Psi = (Gamma * VFI.adj_costs(optK, K, delta, psi)).sum() - C = Y - I - Psi - L_s = get_L_s(w, C, h) - # print('Labor demand and supply = ', L_d, L_s) - MCdist = L_d - L_s - VF_initial = VF - Gamma_initial = Gamma - # print('VF_initial = ', VF_initial[:4, 50:60]) - - - return MCdist - - def Market_Clearing(w, args): (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, sizek, sizez, h, tax_params, VF_initial, Gamma_initial) = args diff --git a/Code/execute.py b/Code/execute.py index 5d0fced..10fbd34 100644 --- a/Code/execute.py +++ b/Code/execute.py @@ -141,30 +141,14 @@ VF_initial = np.zeros((sizez, sizek)) # initial guess at Value Function # initial guess at stationary distribution Gamma_initial = np.ones((sizez, sizek)) * (1 / (sizek * sizez)) -start_time = time.time() -results = opt.bisect(SS.GE_loop, 0.1, 2, args=(alpha_k, alpha_l, delta, psi, - betafirm, K, z, Pi, eta0, - eta1, sizek, sizez, h, - tax_params, VF_initial, - Gamma_initial), - xtol=1e-4, full_output=True) -end_time = time.time() -print('Bisection time = ', end_time - start_time) -print('SS results: ', results) -w = results[0] -# GE_time = time.clock() - start_time -# print('Solving the GE model took ', GE_time, ' seconds to solve') print('SS wage rate = ', w) -gr_args = (alpha_k, alpha_l, delta, psi, - betafirm, K, z, Pi, eta0, eta1, - sizek, sizez, h, tax_params, - VF_initial, Gamma_initial) +gr_args = (alpha_k, alpha_l, delta, psi, betafirm, K, z, Pi, eta0, eta1, + sizek, sizez, h, tax_params, VF_initial, Gamma_initial) start_time = time.time() w = SS.golden_ratio_eqm(0.1, 2, gr_args, tolerance=1e-4) end_time = time.time() -print('Golden Rule time = ', end_time - start_time) +print('Solving the GE model took ', end_time - start_time, ' seconds to solve') print('SS wage rate: ', w) -quit() '''