Merge pull request #4 from jdebacker/equity_costs

Equity costs

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
@@ -45,11 +45,11 @@ def find_SD(PF, Pi, sizez, sizek):
         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:
@@ -59,28 +59,108 @@ def find_SD(PF, Pi, sizez, sizek):
     return Gamma
 
 
-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)
+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 = L_d - L_s
+    # print('Labor demand and supply = ', L_d, L_s)
+    MCdist = np.absolute(L_d - L_s)
 
-    return MCdist
+    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

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
@@ -160,10 +164,10 @@ def VFI(e, betafirm, delta, K, Pi, sizez, sizek, tax_params):
         # 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
 

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
@@ -132,27 +138,29 @@
 Solve for general equilibrium
 ------------------------------------------------------------------------
 '''
-start_time = time.clock()
-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),
-                     xtol=1e-4, full_output=True)
-print('SS results: ', results)
-w = results[0]
-GE_time = time.clock() - start_time
-print('Solving the GE model took ', GE_time, ' seconds to solve')
+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))
 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('Solving the GE model took ', end_time - start_time, ' seconds to solve')
+print('SS wage rate: ', w)
+
 
 '''
 ------------------------------------------------------------------------
 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())
 
 '''