[WIP] Dynamic Firms

ogcore/SS.py

@@ -190,14 +190,31 @@ def inner_loop(outer_loop_vars, p, client):
 
     L = aggr.get_L(nssmat, p, 'SS')
     B = aggr.get_B(bssmat, p, 'SS', False)
-    K_demand_open = firm.get_K(L, p.world_int_rate[-1], p, 'SS')
-    K, K_d, K_f = aggr.get_K_splits(B, K_demand_open, D_d, p.zeta_K[-1])
+    z = firm.get_NPV_depr(r, p, 'SS')
+    V_d = B - D_d
+    if np.any(V_d < 0):
+        print('V_d has negative elements. Setting them ' +
+              'positive to prevent NAN.')
+        V_d = np.fmax(V_d, 0.05 * B)
+    V_f = p.zeta_K[-1] * V_d  # So that foreign equity holdings are some fraction of domestic holdings...
+    # this deviates from CBO method, but theirs doesn't work with dyn firms
+    # could also try to build in realistic repsonse to world and domestic rates - which CBO doesn't have
+    V = V_d + V_f
+    print("V values = ", V, V_d, V_f)
+    K, K_tau = firm.get_K_demand(None, V, None, z, p, 'SS')
+    print('K values = ', K, K_tau)
+
+    # K_demand_open = firm.get_K(L, p.world_int_rate[-1], p, 'SS')
+    # K, K_d, K_f = aggr.get_K_splits(B, K_demand_open, D_d, p.zeta_K[-1])
     Y = firm.get_Y(K, L, p, 'SS')
+    X = firm.get_X(z, K_tau)
     if p.zeta_K[-1] == 1.0:
         new_r = p.world_int_rate[-1]
     else:
-        new_r = firm.get_r(Y, K, p, 'SS')
-    new_w = firm.get_w_from_r(new_r, p, 'SS')
+        # new_r = firm.get_r(Y, K, p, 'SS')
+        new_r = firm.get_r(Y, K, K, V, V, X, X, p, 'SS')
+    # new_w = firm.get_w_from_r(new_r, p, 'SS')
+    new_w = firm.get_w(Y, L, p, 'SS')
 
     b_s = np.array(list(np.zeros(p.J).reshape(1, p.J)) +
                    list(bssmat[:-1, :]))
@@ -325,15 +342,26 @@ def SS_solver(bmat, nmat, r, BQ, TR, factor, Y, p, client,
     Bss = aggr.get_B(bssmat_splus1, p, 'SS', False)
     (Dss, D_d_ss, D_f_ss, new_borrowing, debt_service,
      new_borrowing_f) = fiscal.get_D_ss(r_gov_ss, Y, p)
-    K_demand_open_ss = firm.get_K(Lss, p.world_int_rate[-1], p, 'SS')
-    Kss, K_d_ss, K_f_ss = aggr.get_K_splits(
-        Bss, K_demand_open_ss, D_d_ss, p.zeta_K[-1])
+    # K_demand_open_ss = firm.get_K(Lss, p.world_int_rate[-1], p, 'SS')
+    # Kss, K_d_ss, K_f_ss = aggr.get_K_splits(
+    #     Bss, K_demand_open_ss, D_d_ss, p.zeta_K[-1])
+    zss = firm.get_NPV_depr(rss, p, 'SS')
+    V_d_ss = Bss - D_d_ss
+    if np.any(V_d_ss < 0):
+        print('V_d_ss has negative elements. Setting them ' +
+              'positive to prevent NAN.')
+        V_d_ss = np.fmax(V_d_ss, 0.05 * Bss)
+    V_f_ss = p.zeta_K[-1] * V_d_ss  # So that foreign equity holdings are some fraction of domestic holdings...
+    # this deviates from CBO method, but theirs doesn't work with dyn firms
+    # could also try to build in realistic repsonse to world and domestic rates - which CBO doesn't have
+    Vss = V_d_ss + V_f_ss
+    Kss, K_tau_ss = firm.get_K_demand(None, Vss, None, zss, p, 'SS')
     Yss = firm.get_Y(Kss, Lss, p, 'SS')
     r_hh_ss = aggr.get_r_hh(rss, r_gov_ss, Kss, Dss)
     # Note that implicity in this computation is that immigrants'
     # wealth is all in the form of private capital
-    I_d_ss = aggr.get_I(bssmat_splus1, K_d_ss, K_d_ss, p, 'SS')
-    Iss = aggr.get_I(bssmat_splus1, Kss, Kss, p, 'SS')
+    # I_d_ss = aggr.get_net_I(bssmat_splus1, K_d_ss, K_d_ss, p, 'SS')
+    Iss = aggr.get_net_I(bssmat_splus1, Kss, Kss, p, 'SS')
     wss = new_w
     BQss = new_BQ
     factor_ss = factor
@@ -379,10 +407,14 @@ def SS_solver(bmat, nmat, r, BQ, TR, factor, Y, p, client,
         new_borrowing, debt_service, p)
 
     # Compute total investment (not just domestic)
-    Iss_total = aggr.get_I(None, Kss, Kss, p, 'total_ss')
+    Iss_total = aggr.get_I(Kss, Kss, p.g_n_ss, p.g_y, p.delta)
 
     # solve resource constraint
     # net foreign borrowing
+    K_f_ss = 0  ## temporarily - will come back to open economy case
+    I_d_ss = Iss
+    print('Foreign debt holdings = ', D_f_ss)
+    print('Foreign capital holdings = ', K_f_ss)
     debt_service_f = fiscal.get_debt_service_f(r_hh_ss, D_f_ss)
     RC = aggr.resource_constraint(
         Yss, Css, Gss, I_d_ss, K_f_ss, new_borrowing_f, debt_service_f,
@@ -414,13 +446,15 @@ def SS_solver(bmat, nmat, r, BQ, TR, factor, Y, p, client,
             np.absolute(euler_savings).max())
 
     # Return dictionary of SS results
-    output = {'Kss': Kss, 'K_f_ss': K_f_ss, 'K_d_ss': K_d_ss,
+    output = {'Kss': Kss, 'K_tau_ss': K_tau_ss, 'Vss': Vss,
+              'V_f_ss': V_f_ss, 'V_d_ss': V_d_ss,
               'Bss': Bss, 'Lss': Lss, 'Css': Css, 'Iss': Iss,
-              'Iss_total': Iss_total, 'I_d_ss': I_d_ss, 'nssmat': nssmat,
+              'Iss_total': Iss_total, 'nssmat': nssmat,
               'Yss': Yss, 'Dss': Dss, 'D_f_ss': D_f_ss,
               'D_d_ss': D_d_ss, 'wss': wss, 'rss': rss,
-              'total_taxes_ss': taxss, 'ubissmat': ubissmat,
-              'r_gov_ss': r_gov_ss, 'r_hh_ss': r_hh_ss, 'theta': theta,
+              'r_gov_ss': r_gov_ss, 'r_hh_ss': r_hh_ss, 'zss': zss,
+              'theta': theta, 'total_taxes_ss': taxss,
+              'ubissmat': ubissmat,
               'BQss': BQss, 'factor_ss': factor_ss, 'bssmat_s': bssmat_s,
               'cssmat': cssmat, 'bssmat_splus1': bssmat_splus1,
               'yss_before_tax_mat': yss_before_tax_mat,

ogcore/TPI.py

@@ -80,6 +80,11 @@ def get_initial_SS_values(p):
     # placeholder for first-period wealth
     B0 = aggr.get_B(initial_b, p, 'SS', True)
 
+    # Get initial value of the capital stock and initial value of the
+    # depreciable basis for the capital stock
+    K0 = ss_baseline_vars['Kss']
+    K_tau0 = ss_baseline_vars['K_tau_ss']
+
     b_sinit = np.array(list(np.zeros(p.J).reshape(1, p.J)) +
                        list(initial_b[:-1]))
     b_splus1init = initial_b
@@ -93,7 +98,7 @@ def get_initial_SS_values(p):
     else:
         D0_baseline = None
 
-    initial_values = (B0, b_sinit, b_splus1init, factor, initial_b,
+    initial_values = (B0, K0, K_tau0, b_sinit, b_splus1init, factor, initial_b,
                       initial_n)
     baseline_values = (TRbaseline, Gbaseline, D0_baseline)
 
@@ -268,7 +273,8 @@ def inner_loop(guesses, outer_loop_vars, initial_values, ubi, j, ind, p):
             * n_mat (Numpy array): labor supply amounts, size = TxS
 
     '''
-    (K0, b_sinit, b_splus1init, factor, initial_b, initial_n) =\
+    # unpack variables and parameters pass to function
+    (B0, K0, K_tau0, b_sinit, b_splus1init, factor, initial_b, initial_n) =\
         initial_values
     guesses_b, guesses_n = guesses
     r, w, r_hh, BQ, TR, theta = outer_loop_vars
@@ -389,8 +395,8 @@ def run_TPI(p, client=None):
     '''
     # unpack tuples of parameters
     initial_values, ss_vars, theta, baseline_values = get_initial_SS_values(p)
-    (B0, b_sinit, b_splus1init, factor, initial_b, initial_n) =\
-        initial_values
+    (B0, K0, K_tau0, b_sinit, b_splus1init, factor, initial_b,
+     initial_n) = initial_values
     (TRbaseline, Gbaseline, D0_baseline) = baseline_values
 
     # Create time path of UBI household benefits and aggregate UBI outlays
@@ -416,17 +422,27 @@ def run_TPI(p, client=None):
     L_init[:p.T] = aggr.get_L(n_mat[:p.T], p, 'TPI')
     B_init[1:p.T] = aggr.get_B(b_mat[:p.T], p, 'TPI', False)[:p.T - 1]
     B_init[0] = B0
-    K_init = B_init * ss_vars['Kss'] / ss_vars['Bss']
-    K = K_init
-    K_d = K_init * ss_vars['K_d_ss'] / ss_vars['Kss']
-    K_f = K_init * ss_vars['K_f_ss'] / ss_vars['Kss']
+    V_init = B_init * ss_vars['Vss'] / ss_vars['Bss']
+    V = V_init
+    V_d = V_init * ss_vars['V_d_ss'] / ss_vars['Vss']
+    V_f = V_init * ss_vars['V_f_ss'] / ss_vars['Vss']
+    # compute z
+    z = np.ones(p.T + p.S) * ss_vars['zss']
+    K = np.ones(p.T + p.S) * ss_vars['Kss']
+    K_tau = np.ones(p.T + p.S) * ss_vars['K_tau_ss']
+    K[:p.T], K_tau[:p.T] = firm.get_K_demand(
+        K0, V[:p.T], K_tau0, z[:p.T], p, 'TPI')
+    X = firm.get_X(z, K_tau)
     L = L_init
     B = B_init
     Y = np.zeros_like(K)
     Y[:p.T] = firm.get_Y(K[:p.T], L[:p.T], p, 'TPI')
     Y[p.T:] = ss_vars['Yss']
     r = np.zeros_like(Y)
-    r[:p.T] = firm.get_r(Y[:p.T], K[:p.T], p, 'TPI')
+    r[:p.T] = firm.get_r(
+            Y[:p.T], K[:p.T], K[1:p.T+1], V[:p.T], V[1:p.T+1], X[:p.T],
+            X[1:p.T+1], p, 'TPI')
+    # r[:p.T] = firm.get_r(Y[:p.T], K[:p.T], p, 'TPI')
     r[p.T:] = ss_vars['rss']
     # For case where economy is small open econ
     r[p.zeta_K == 1] = p.world_int_rate[p.zeta_K == 1]
@@ -493,7 +509,7 @@ def run_TPI(p, client=None):
 
         r_hh[:p.T] = aggr.get_r_hh(r[:p.T], r_gov[:p.T], K[:p.T],
                                    D[:p.T])
-
+        print("INTREST RATES = ", r[:5], r_gov[:5], r_hh[:5])
         outer_loop_vars = (r, w, r_hh, BQ, TR, theta)
 
         euler_errors = np.zeros((p.T, 2 * p.S, p.J))
@@ -555,22 +571,44 @@ def run_TPI(p, client=None):
          debt_service, new_borrowing_f) =\
             fiscal.D_G_path(r_gov, dg_fixed_values, p)
         L[:p.T] = aggr.get_L(n_mat[:p.T], p, 'TPI')
+        print('Initial B, V, K = ', B[0], V[0], K[0])
         B[1:p.T] = aggr.get_B(bmat_splus1[:p.T], p, 'TPI',
                               False)[:p.T - 1]
-        K_demand_open = firm.get_K(
-            L[:p.T], p.world_int_rate[:p.T], p, 'TPI')
-        K[:p.T], K_d[:p.T], K_f[:p.T] = aggr.get_K_splits(
-            B[:p.T], K_demand_open, D_d[:p.T], p.zeta_K[:p.T])
+
+        z[:p.T] = firm.get_NPV_depr(r, p, 'TPI')[:p.T]
+        V_d[:p.T] = B[:p.T] - D_d[:p.T]
+        if np.any(V_d < 0):
+            print('V_d has negative elements. Setting them ' +
+                  'positive to prevent NAN.')
+            V_d = np.fmax(V_d, 0.05 * B)
+        V_f[:p.T] = p.zeta_K[:p.T] * V_d[:p.T]  # So that foreign equity holdings are some fraction of domestic holdings...
+        # this deviates from CBO method, but theirs doesn't work with dyn firms
+        # could also try to build in realistic repsonse to world and domestic rates - which CBO doesn't have
+        V[:p.T] = V_d[:p.T] + V_f[:p.T]
+        X = firm.get_X(z, K_tau)
+        K0 = V[0]  # cheat for now with zero adj costs -- otherwise issue because for some guesses of Y, K0 and V[0] not equal and yield negative interest rate
+        K[:p.T], K_tau[:p.T] = firm.get_K_demand(K0, V, K_tau0, z, p, 'TPI')
+
+        # K_demand_open = firm.get_K(
+        #     L[:p.T], p.world_int_rate[:p.T], p, 'TPI')
+        # K[:p.T], K_d[:p.T], K_f[:p.T] = aggr.get_K_splits(
+        #     B[:p.T], K_demand_open, D_d[:p.T], p.zeta_K[:p.T])
         Ynew = firm.get_Y(K[:p.T], L[:p.T], p, 'TPI')
         rnew = r.copy()
-        rnew[:p.T] = firm.get_r(Ynew[:p.T], K[:p.T], p, 'TPI')
+        # rnew[:p.T] = firm.get_r(Ynew[:p.T], K[:p.T], p, 'TPI')
+        rnew[:p.T] = firm.get_r(
+            Y[:p.T], K[:p.T], K[1:p.T+1], V[:p.T], V[1:p.T+1], X[:p.T],
+            X[1:p.T+1], p, 'TPI')
+        print('VARS ALONG THE WAY: ', r[:5], K[:5], Y[:5], V[:5], X[:5])
         # For case where economy is small open econ
         r[p.zeta_K == 1] = p.world_int_rate[p.zeta_K == 1]
         r_gov_new = fiscal.get_r_gov(rnew, p)
         r_hh_new = aggr.get_r_hh(rnew[:p.T], r_gov_new[:p.T], K[:p.T],
                                  Dnew[:p.T])
         # compute w
-        wnew = firm.get_w_from_r(rnew[:p.T], p, 'TPI')
+        # wnew = firm.get_w_from_r(rnew[:p.T], p, 'TPI')
+        wnew = w
+        wnew[:p.T] = firm.get_MPL(Ynew[:p.T], L[:p.T], p, 'TPI')
 
         b_mat_shift = np.append(np.reshape(initial_b, (1, p.S, p.J)),
                                 b_mat[:p.T - 1, :, :], axis=0)
@@ -591,7 +629,8 @@ def run_TPI(p, client=None):
             agg_pension_outlays[:p.T], UBI_outlays[:p.T], p, 'TPI')
 
         # update vars for next iteration
-        w[:p.T] = wnew[:p.T]
+        # w[:p.T] = wnew[:p.T]
+        w[:p.T] = utils.convex_combo(wnew[:p.T], w[:p.T], p.nu)
         r[:p.T] = utils.convex_combo(rnew[:p.T], r[:p.T], p.nu)
         BQ[:p.T] = utils.convex_combo(BQnew[:p.T], BQ[:p.T], p.nu)
         D[:p.T] = Dnew[:p.T]
@@ -665,18 +704,19 @@ def run_TPI(p, client=None):
     C = aggr.get_C(c_mat, p, 'TPI')
     # Note that implicity in this computation is that immigrants'
     # wealth is all in the form of private capital
-    I_d = aggr.get_I(bmat_splus1[:p.T], K_d[1:p.T + 1], K_d[:p.T], p,
-                     'TPI')
-    I = aggr.get_I(bmat_splus1[:p.T], K[1:p.T + 1], K[:p.T], p, 'TPI')
+    # I_d = aggr.get_net_I(
+    #     bmat_splus1[:p.T], K_d[1:p.T + 1], K_d[:p.T], p, 'TPI')
+    I = aggr.get_net_I(
+        bmat_splus1[:p.T], K[1:p.T + 1], K[:p.T], p, 'TPI')
     # solve resource constraint
     # foreign debt service costs
     debt_service_f = fiscal.get_debt_service_f(r_hh, D_f)
     RC_error = aggr.resource_constraint(
-        Y[:p.T - 1], C[:p.T - 1], G[:p.T - 1], I_d[:p.T - 1],
-        K_f[:p.T - 1], new_borrowing_f[:p.T - 1],
+        Y[:p.T - 1], C[:p.T - 1], G[:p.T - 1], I[:p.T - 1],
+        V_f[:p.T - 1], new_borrowing_f[:p.T - 1],
         debt_service_f[:p.T - 1], r_hh[:p.T - 1], p)
     # Compute total investment (not just domestic)
-    I_total = aggr.get_I(None, K[1:p.T + 1], K[:p.T], p, 'total_tpi')
+    I_total = aggr.get_net_I(None, K[1:p.T + 1], K[:p.T], p, 'total_tpi')
 
     # Compute resource constraint error
     rce_max = np.amax(np.abs(RC_error))
@@ -699,9 +739,9 @@ def run_TPI(p, client=None):
     ------------------------------------------------------------------------
     '''
 
-    output = {'Y': Y[:p.T], 'B': B, 'K': K, 'K_f': K_f, 'K_d': K_d,
+    output = {'Y': Y[:p.T], 'B': B, 'K': K, 'V_f': V_f, 'V_d': V_d,
               'L': L, 'C': C, 'I': I,
-              'I_total': I_total, 'I_d': I_d, 'BQ': BQ,
+              'I_total': I_total, 'BQ': BQ,
               'total_tax_revenue': total_tax_revenue,
               'business_tax_revenue': business_tax_revenue,
               'iit_payroll_tax_revenue': iit_payroll_tax_revenue,

ogcore/aggregates.py

@@ -43,9 +43,10 @@ def get_L(n, p, method):
     return L
 
 
-def get_I(b_splus1, K_p1, K, p, method):
+def get_net_I(b_splus1, K_p1, K, p, method):
     r'''
-    Calculate aggregate investment.
+    Calculate aggregate investment, accounting for flows from
+    immigrants.
 
     .. math::
         I_{t} = (1 + g_{n,t+1})e^{g_{y}}(K_{t+1} - \sum_{s=E}^{E+S}
@@ -93,6 +94,29 @@ def get_I(b_splus1, K_p1, K, p, method):
     return aggI
 
 
+def get_I(K_p1, K, g_np1, g_y, delta):
+    r'''
+    Calculate gross aggregate investment.
+
+    .. math::
+        I_{t} = (1 + g_{n,t+1})e^{g_{y}}(K_{t+1} - (1 - \delta)K_{t}
+
+    Args:
+        K_p1 (array_like): aggregate capital, one period ahead
+        K (array_like): aggregate capital
+        g_np1 (array_like): population growth rate one period ahead
+        g_y (scalar): growth rate in labor augmenting tech change
+        delta (scalar): rate of economic depreciation
+
+    Returns:
+        aggI (array_like): aggregate investment
+
+    '''
+    aggI = (1 + g_np1) * np.exp(g_y) * K_p1 - (1.0 - delta) * K
+
+    return aggI
+
+
 def get_B(b, p, method, preTP):
     r'''
     Calculate aggregate savings

ogcore/default_parameters.json

@@ -381,6 +381,25 @@
             }
         }
     },
+    "psi": {
+        "title": "Scale parameter in quadratic capital adjustment cost function",
+        "description": "Scale parameter in quadratic capital adjustment cost function",
+        "section_1": "Firm Parameters",
+        "section_2": "Investment costs",
+        "notes": "If psi = 0 then the dynamic firm problem collapses to the static firm problem.",
+        "type": "float",
+        "value": [
+            {
+                "value": 0.0
+            }
+        ],
+        "validators": {
+            "range": {
+                "min": 0.0,
+                "max": 10.0
+            }
+        }
+    },
     "Z": {
         "title": "Total factor productivity in firm production function",
         "description": "Total factor productivity in firm production function.  Set value for base year, click '+' to add value for next year.  All future years not specified are set to last value entered.",