Skip to content

Commit

Permalink
Merge pull request #6 from jdebacker/add_debt
Browse files Browse the repository at this point in the history 
Add debt
  • Loading branch information
@jdebacker
jdebacker authored Oct 21, 2017
2 parents 8aecece + 3023f2f commit 8eae279
Show file tree
Hide file tree
Showing 8 changed files with 1,207 additions and 319 deletions.
70 changes: 44 additions & 26 deletions Code/SS.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,10 @@


@numba.jit
def find_SD(PF, Pi, sizez, sizek, Gamma_initial):
def find_SD(PF_k, PF_b, Pi, sizez, sizek, sizeb, Gamma_initial):
'''
------------------------------------------------------------------------
Compute the stationary distribution of firms over (k, z)
Compute the stationary distribution of firms over (z, k)
------------------------------------------------------------------------
SDtol = tolerance required for convergence of SD
SDdist = distance between last two distributions
Expand All @@ -30,17 +30,19 @@ def find_SD(PF, Pi, sizez, sizek, Gamma_initial):
------------------------------------------------------------------------
'''
Gamma = Gamma_initial
SDtol = 1e-12
SDtol = 1e-8#1e-12
SDdist = 7
SDiter = 0
SDmaxiter = 1000
SDmaxiter = 2000
while SDdist > SDtol and SDmaxiter > SDiter:
HGamma = np.zeros((sizez, sizek))
HGamma = np.zeros((sizez, sizek, sizeb))
for i in range(sizez): # z
for j in range(sizek): # k
for m in range(sizez): # z'
HGamma[m, PF[i, j]] = \
HGamma[m, PF[i, j]] + Pi[i, m] * Gamma[i, j]
for m in range(sizeb): # b
for ii in range(sizez): # z'
HGamma[ii, PF_k[i, j, m], PF_b[i, j, m]] = \
(HGamma[ii, PF_k[i, j, m], PF_b[i, j, m]] +
Pi[i, ii] * Gamma[i, j, m])
SDdist = (np.absolute(HGamma - Gamma)).max()
Gamma = HGamma
SDiter += 1
Expand All @@ -52,28 +54,41 @@ def find_SD(PF, Pi, sizez, sizek, Gamma_initial):
print('Stationary distribution did not converge')

# Check if state space is binding
if Gamma.sum(axis=0)[-1] > 0.002:
if (Gamma.sum(axis=0)).sum(axis=1)[-1] > 0.002:
print('Stationary distribution is binding on k-grid. Consider ' +
'increasing the upper bound.')
if (Gamma.sum(axis=0)).sum(axis=0)[-1] > 0.01:
print('Stationary distribution is binding on b-grid. Consider ' +
'increasing the upper bound.')
if (Gamma.sum(axis=0)).sum(axis=0)[0] > 0.01:
print('Stationary distribution is binding on b-grid. Consider ' +
'decreasing the lower bound.')

return Gamma


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
(r, alpha_k, alpha_l, delta, psi, betafirm, kgrid, zgrid, bgrid, Pi,
eta0, eta1, eta2, s, sizek, sizez, sizeb, 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()
op, e, l_d, y, eta, collateral_constraint = VFI.get_firmobjects(r, w, zgrid, kgrid, bgrid,
alpha_k, alpha_l, delta,
psi, eta0, eta1, eta2, s,
sizez, sizek, sizeb, tax_params)
VF, PF_k, PF_b, optK, optI, optB = VFI.VFI(e, eta, collateral_constraint, betafirm, delta,
kgrid, bgrid, Pi, sizez, sizek,
sizeb, tax_params, VF_initial)
Gamma = find_SD(PF_k, PF_b, Pi, sizez, sizek, sizeb, Gamma_initial)
L_d = (Gamma.sum(axis=2) * l_d).sum()
Y = (Gamma.sum(axis=2) * y).sum()
I = (Gamma * optI).sum()
Psi = (Gamma * VFI.adj_costs(optK, K, delta, psi)).sum()
k3grid = np.tile(np.reshape(kgrid, (1, sizek, 1)), (sizez, 1, sizeb))
Psi = (Gamma * VFI.adj_costs(optK, k3grid, delta, psi)).sum()
C = Y - I - Psi
# note that financial frictions not here- they aren't real costs,
# rather they are costs paid by firms and recieved for financial
# interemediaries and flow to households as income
L_s = get_L_s(w, C, h)
# print('Labor demand and supply = ', L_d, L_s)
MCdist = np.absolute(L_d - L_s)
Expand All @@ -91,8 +106,9 @@ 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
(r, alpha_k, alpha_l, delta, psi, betafirm, kgrid, zgrid, bgrid, Pi,
eta0, eta1, eta2, s, sizek, sizez, sizeb, 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
Expand Down Expand Up @@ -128,8 +144,9 @@ def golden_ratio_eqm(lb, ub, args, tolerance=1e-4):
# 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)
args = (r, alpha_k, alpha_l, delta, psi, betafirm, kgrid, zgrid,
bgrid, Pi, eta0, eta1, eta2, s, sizek, sizez, sizeb, h,
tax_params, VF1, Gamma1)
f1, VF1, Gamma1 = Market_Clearing(x1, args)
else:
# print('f2 < f1')
Expand All @@ -151,8 +168,9 @@ def golden_ratio_eqm(lb, ub, args, tolerance=1e-4):
# 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)
args = (r, alpha_k, alpha_l, delta, psi, betafirm, kgrid, zgrid,
bgrid, Pi, eta0, eta1, eta2, s, sizek, sizez, sizeb, 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
Expand Down
168 changes: 128 additions & 40 deletions Code/VFI.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,37 @@
# imports
import numba
import numpy as np
import time


@numba.jit
def create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, Vmat, tax_params):
def create_EV(Pi, V, sizez, sizek, sizeb):
'''
Compute expectation of continuation value function.
Args:
Pi: 2D array, transition probabilities for exogenous state var
V: 3D array, value function
Returns:
EV: 3D array, expected value function (conditional on z)
'''
EV = np.zeros_like(V)
# start = time.time()
for i in range(sizez): # loop over z
for jj in range(sizek): # loop over k'
for mm in range(sizeb): # loop over b'
for ii in range(sizez): # loop over z'
EV[i, jj, mm] = EV[i, jj, mm] + Pi[i, ii] * V[ii, jj, mm]
# end = time.time()
# print('EV loop takes ', end-start, ' seconds to complete.')

return EV



@numba.jit
def create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, sizeb, tax_params):
'''
------------------------------------------------------------------------
This function loops over the state and control variables, operating on the
Expand Down Expand Up @@ -46,14 +73,21 @@ def create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, Vmat, tax_params):
RETURNS: Vmat
------------------------------------------------------------------------
'''
tau_i, tau_d, tau_g, tau_c = tax_params
tau_i, tau_d, tau_g, tau_c, f_e, f_b = tax_params
# initialize Vmat array
# start = time.time()
Vmat = np.empty((sizez, sizek, sizek, sizeb, sizeb))
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]) * (e[i, j, m] >= 0) +
((e[i, j, m] + eta[i, j, m]) *
(e[i, j, m] < 0)) + betafirm * EV[i, m])
for jj in range(sizek): # loop over k'
for m in range(sizeb): # loop over b
for mm in range(sizeb): # loop over b'
Vmat[i, j, jj, m, mm] = ((((1 - tau_d) / (1 - tau_g)) *
e[i, j, jj, m, mm]) * (e[i, j, jj, m, mm] >= 0) +
((e[i, j, jj, m, mm] + eta[i, j, jj, m, mm]) *
(e[i, j, jj, m, mm] < 0)) + betafirm * EV[i, jj, mm])
# end = time.time()
# print('Vmat loop takes ', end-start, ' seconds to complete.')

return Vmat

Expand All @@ -75,8 +109,8 @@ def adj_costs(kprime, k, delta, psi):


@numba.jit
def get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, eta0, eta1,
sizez, sizek, tax_params):
def get_firmobjects(r, w, zgrid, kgrid, bgrid, alpha_k, alpha_l, delta, psi, eta0, eta1,
eta2, s, sizez, sizek, sizeb, tax_params):
'''
-------------------------------------------------------------------------
Generating possible cash flow levels
Expand All @@ -92,35 +126,59 @@ def get_firmobjects(w, z, K, alpha_k, alpha_l, delta, psi, eta0, eta1,
today (state), and choice of capital stock tomorrow (control)
-------------------------------------------------------------------------
'''
tau_i, tau_d, tau_g, tau_c = tax_params
tau_i, tau_d, tau_g, tau_c, f_e, f_b = tax_params
# Initialize arrays
op = np.empty((sizez, sizek))
l_d = np.empty((sizez, sizek))
y = np.empty((sizez, sizek))
e = np.empty((sizez, sizek, sizek))
for i in range(sizez):
for j in range(sizek):
e = np.empty((sizez, sizek, sizek, sizeb, sizeb))
collateral_constraint = np.empty((sizez, sizek, sizek, sizeb, sizeb))
# start = time.time()
for i in range(sizez): # loop over z
for j in range(sizek): # loop over k
op[i, j] = ((1 - alpha_l) * ((alpha_l / w) **
(alpha_l / (1 - alpha_l))) *
((z[i] * (K[j] ** alpha_k)) **
((zgrid[i] * (kgrid[j] ** alpha_k)) **
(1 / (1 - alpha_l))))
l_d[i, j] = (((alpha_l / w) ** (1 / (1 - alpha_l))) *
(z[i] ** (1 / (1 - alpha_l))) *
(K[j] ** (alpha_k / (1 - alpha_l))))
y[i, j] = z[i] * (K[j] ** alpha_k) * (l_d[i, j] ** alpha_l)
for m in range(sizek):
e[i, j, m] = ((1 - tau_c) * op[i, j] + (delta * tau_c
* K[j]) - K[m]
+ ((1 - delta)
* K[j]) -
adj_costs(K[m], K[j], delta, psi))
(zgrid[i] ** (1 / (1 - alpha_l))) *
(kgrid[j] ** (alpha_k / (1 - alpha_l))))
y[i, j] = zgrid[i] * (kgrid[j] ** alpha_k) * (l_d[i, j] ** alpha_l)
for m in range(sizeb): # loop over b
for jj in range(sizek): # loop over k'
for mm in range(sizeb): # loop over b'
e[i, j, jj, m, mm] =\
(((1 - tau_c) * op[i, j]) +
(delta * (1 - f_e) * tau_c * kgrid[j]) +
(f_e * tau_c * (kgrid[jj] > ((1 - delta) *
kgrid[j])) *
(kgrid[jj] - ((1 - delta) * kgrid[j]))) -
(kgrid[jj] - ((1 - delta) * kgrid[j])) -
adj_costs(kgrid[jj], kgrid[j], delta, psi) +
bgrid[mm] - ((1 + r) * bgrid[m]) +
(r * tau_c * bgrid[m]) -
(r * tau_c * (1 - f_b) * bgrid[m] *
(bgrid[m] > 0)))
# collateral_constraint[i, j, jj, m, mm] =\
# (((1-tau_c) * op[0, jj] + tau_c * delta *
# kgrid[jj] + s * kgrid[jj]) <=
# (((1 + r) * bgrid[mm]) - (r * tau_c * bgrid[mm])))
collateral_constraint[i, j, jj, m, mm] =\
(((1 + r) * bgrid[mm] - (tau_c * r * bgrid[mm])
+ (r * tau_c * (1 - f_b) * bgrid[mm] *
(bgrid[mm] > 0))) >
(((1 - tau_c) * op[0, jj]) +
((1 - f_e) * tau_c * delta * kgrid[jj]) +
s * kgrid[jj]))
eta = (-1 * eta0 + eta1 * e - eta2 * (e ** 2)) * (e < 0)
# end = time.time()
# print('Firm objects loop takes ', end-start, ' seconds to complete.')

eta = (eta0 + eta1 * e) * (e < 0)
return op, e, l_d, y, eta, collateral_constraint

return op, e, l_d, y, eta


def VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, tax_params, VF_initial):
def VFI(e, eta, collateral_constraint, betafirm, delta, kgrid, bgrid, Pi, sizez, sizek, sizeb,
tax_params, VF_initial):
'''
------------------------------------------------------------------------
Value Function Iteration
Expand All @@ -143,31 +201,59 @@ def VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, tax_params, VF_initial):
possible value of the state variables (k)
------------------------------------------------------------------------
'''
tau_i, tau_d, tau_g, tau_c = tax_params
tau_i, tau_d, tau_g, tau_c, f_e, f_b = tax_params
VFtol = 1e-6
VFdist = 7.0
VFmaxiter = 3000
V = VF_initial
Vmat = np.empty((sizez, sizek, sizek)) # initialize Vmat matrix
Vstore = np.empty((sizez, sizek, VFmaxiter)) # initialize Vstore array
#Vstore = np.empty((sizez, sizek, sizeb, VFmaxiter)) # initialize Vstore array
VFiter = 1
# while VFdist > VFtol and VFiter < VFmaxiter:
# TV = V
# EV = create_EV(Pi, V, sizez, sizek, sizeb) # expected VF (expectation over z')
# Vmat = create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, sizeb,
# tax_params) + (collateral_constraint * -1000000000)
#
# #Vstore[:, :, :, VFiter] = V.reshape(sizez, sizek, sizeb) # store value function
# # at each iteration for graphing later
# # apply max operator to Vmat (to get V(z,k,b))
# V = (Vmat.max(axis=4)).max(axis=2)
# PF_k = np.argmax(Vmat.max(axis=4), axis=2)
# PF_b = np.argmax(Vmat.max(axis=2), axis=3)
# VFdist = (np.absolute(V - TV)).max() # check distance between value
# # function for this iteration and value function from past iteration
# # print('VF iteration: ', VFiter)
# VFiter += 1

VFflag = 0
PF_k_old = np.zeros_like(V)
PF_b_old = np.zeros_like(V)
while VFdist > VFtol and VFiter < VFmaxiter:
TV = V
EV = np.dot(Pi, V) # expected VF (expectation over z')
Vmat = create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, Vmat, tax_params)
EV = create_EV(Pi, V, sizez, sizek, sizeb) # expected VF (expectation over z')
Vmat = create_Vmat(EV, e, eta, betafirm, Pi, sizez, sizek, sizeb,
tax_params) + (collateral_constraint * -1000000000)

if VFiter%10 == 0:

Vstore[:, :, VFiter] = V.reshape(sizez, sizek) # store value function
# at each iteration for graphing later
V = Vmat.max(axis=2) # apply max operator to Vmat (to get V(k))
PF = np.argmax(Vmat, axis=2)
V = (Vmat.max(axis=4)).max(axis=2)
PF_k = np.argmax(Vmat.max(axis=4), axis=2)
PF_b = np.argmax(Vmat.max(axis=2), axis=3)
if (np.absolute(PF_k-PF_k_old).max() == 0) & (np.absolute(PF_b-PF_b_old).max() == 0):
VFiter = VFmaxiter
VFflag = 1
else:
PF_k_old = PF_k
PF_b_old = PF_b
V = (Vmat.max(axis=4)).max(axis=2)
VFdist = (np.absolute(V - TV)).max() # check distance between value
# function for this iteration and value function from past iteration
# print('VF iteration: ', VFiter)
VFiter += 1

if VFiter < VFmaxiter:
print('Value function converged after this many iterations:', VFiter)
else:
elif VFflag == 0:
print('Value function did not converge')

VF = V # solution to the functional equation
Expand All @@ -180,7 +266,9 @@ def VFI(e, eta, betafirm, delta, K, Pi, sizez, sizek, tax_params, VF_initial):
optI = (sizez, sizek) vector, optimal choice of investment for each (z, k)
------------------------------------------------------------------------
'''
optK = K[PF]
optI = optK - (1 - delta) * K
optK = kgrid[PF_k]
k3grid = np.tile(np.reshape(kgrid, (1, sizek, 1)), (sizez, 1, sizeb))
optI = optK - (1 - delta) * k3grid
optB = bgrid[PF_b]

return VF, PF, optK, optI
return VF, PF_k, PF_b, optK, optI, optB
Loading

0 comments on commit 8eae279

Please sign in to comment.