Merge pull request PSLmodels#88 from evan-magnusson/firmprob

[Development_firmproblem] [WIP] Multitude of bug fixes, mostly involving resource constraint violations

Python/Development_firmproblem/SS.py

@@ -1,6 +1,6 @@
 '''
 ------------------------------------------------------------------------
-Last updated: 6/4/2015
+Last updated: 7/13/2015
 
 Calculates steady state of OLG model with S age cohorts.
 
@@ -42,7 +42,7 @@
 S            = number of periods an individual lives
 J            = number of different ability groups
 T            = number of time periods until steady state is reached
-lambdas_init  = percent of each age cohort in each ability group
+lambdas  = percent of each age cohort in each ability group
 starting_age = age of first members of cohort
 ending age   = age of the last members of cohort
 E            = number of cohorts before S=1
@@ -91,7 +91,6 @@
     for key in variables:
         globals()[key] = variables[key]
 
-
 '''
 ------------------------------------------------------------------------
     Define Functions
@@ -102,22 +101,22 @@
 income_tax_params = [a_tax_income, b_tax_income, c_tax_income, d_tax_income]
 wealth_tax_params = [h_wealth, p_wealth, m_wealth]
 ellipse_params = [b_ellipse, upsilon]
-parameters = [J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, tau_payroll, retire, mean_income_data] + income_tax_params + wealth_tax_params + ellipse_params
+parameters = [J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, g_n_ss, tau_payroll, retire, mean_income_data] + income_tax_params + wealth_tax_params + ellipse_params
 iterative_params = [maxiter, mindist_SS]
 
 # Functions
 
 
 def Euler_equation_solver(guesses, r, w, T_H, factor, j, params, chi_b, chi_n, tau_bq, rho, lambdas, weights, e):
-    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
+    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, g_n_ss, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
     b_guess = np.array(guesses[:S])
     n_guess = np.array(guesses[S:])
     b_s = np.array([0] + list(b_guess[:-1]))
     b_splus1 = b_guess
     b_splus2 = np.array(list(b_guess[1:]) + [0])
 
-    BQ = (1.0/(1+g_n_ss))*house.get_BQ(r, b_splus1, weights[:, j], rho)
-    theta = tax.replacement_rate_vals(n_guess, w, factor, e[:,j], J, weights[:, j])
+    BQ = house.get_BQ(r, b_splus1, weights, lambdas[j], rho, g_n_ss)
+    theta = tax.replacement_rate_vals(n_guess, w, factor, e[:,j], J, weights, lambdas[j])
 
     error1 = house.euler_savings_func(w, r, e[:, j], n_guess, b_s, b_splus1, b_splus2, BQ, factor, T_H, chi_b[j], params, theta, tau_bq[j], rho, lambdas[j])
     error2 = house.euler_labor_leisure_func(w, r, e[:, j], n_guess, b_s, b_splus1, BQ, factor, T_H, chi_n, params, theta, tau_bq[j], lambdas[j])
@@ -138,7 +137,7 @@ def Euler_equation_solver(guesses, r, w, T_H, factor, j, params, chi_b, chi_n, t
 
 
 def SS_solver(b_guess_init, n_guess_init, wguess, rguess, T_Hguess, factorguess, chi_n, chi_b, params, iterative_params, tau_bq, rho, lambdas, weights, e):
-    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
+    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, g_n_ss, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
     maxiter, mindist_SS = iterative_params
     w = wguess
     r = rguess
@@ -160,17 +159,17 @@ def SS_solver(b_guess_init, n_guess_init, wguess, rguess, T_Hguess, factorguess,
             nssmat[:,j] = solutions[S:]
             # print np.array(Euler_equation_solver(np.append(bssmat[:, j], nssmat[:, j]), r, w, T_H, factor, j, params, chi_b, chi_n, theta, tau_bq, rho, lambdas, e)).max()
 
-        K = (1.0/(1+g_n_ss)) * house.get_K(bssmat, weights)
-        L = firm.get_L(e, nssmat, weights)
+        K = house.get_K(bssmat, weights.reshape(S, 1), lambdas.reshape(1, J), g_n_ss)
+        L = firm.get_L(e, nssmat, weights.reshape(S, 1), lambdas.reshape(1, J))
         Y = firm.get_Y(K, L, params)
         new_r = firm.get_r(Y, K, params)
         new_w = firm.get_w(Y, L, params)
         b_s = np.array(list(np.zeros(J).reshape(1, J)) + list(bssmat[:-1, :]))
-        average_income_model = ((new_r * b_s + new_w * e * nssmat) * weights).sum()
+        average_income_model = ((new_r * b_s + new_w * e * nssmat) * weights.reshape(S, 1) * lambdas.reshape(1, J)).sum()
         new_factor = mean_income_data / average_income_model 
-        new_BQ = (1.0/(1+g_n_ss))*house.get_BQ(new_r, bssmat, weights, rho.reshape(S, 1))
-        theta = tax.replacement_rate_vals(nssmat, new_w, new_factor, e, J, weights)
-        new_T_H = tax.get_lump_sum(new_r, b_s, new_w, e, nssmat, new_BQ, lambdas, factor, weights, 'SS', params, theta, tau_bq)
+        new_BQ = house.get_BQ(new_r, bssmat, weights.reshape(S, 1), lambdas.reshape(1, J), rho.reshape(S, 1), g_n_ss)
+        theta = tax.replacement_rate_vals(nssmat, new_w, new_factor, e, J, weights.reshape(S, 1), lambdas)
+        new_T_H = tax.get_lump_sum(new_r, b_s, new_w, e, nssmat, new_BQ, lambdas.reshape(1, J), factor, weights.reshape(S, 1), 'SS', params, theta, tau_bq)
 
         r = misc_funcs.convex_combo(new_r, r, params)
         w = misc_funcs.convex_combo(new_w, w, params)
@@ -212,7 +211,7 @@ def function_to_minimize(chi_params_scalars, chi_params_init, params, weights_SS
         The max absolute deviation between the actual and simulated
             wealth moments
     '''
-    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
+    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, g_n_ss, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
     chi_params_init *= chi_params_scalars
     # print 'Print Chi_b: ', chi_params_init[:J]
     # print 'Scaling vals:', chi_params_scalars[:J]
@@ -337,24 +336,23 @@ def function_to_minimize(chi_params_scalars, chi_params_init, params, weights_SS
 nssmat = solutions[S * J:2*S*J].reshape(S, J)
 wss, rss, factor_ss, T_Hss = solutions[2*S*J:]
 
-Kss = (1.0/(1+g_n_ss)) * house.get_K(bssmat_splus1, omega_SS)
-Lss = firm.get_L(e, nssmat, omega_SS)
+Kss = house.get_K(bssmat_splus1, omega_SS.reshape(S, 1), lambdas, g_n_ss)
+Lss = firm.get_L(e, nssmat, omega_SS.reshape(S, 1), lambdas)
 Yss = firm.get_Y(Kss, Lss, parameters)
 
-Iss = delta*Kss
+Iss = firm.get_I(Kss, Kss, delta, g_y, g_n_ss)
 
-theta = tax.replacement_rate_vals(nssmat, wss, factor_ss, e, J, omega_SS)
-BQss = (1.0/(1+g_n_ss))*house.get_BQ(rss, bssmat_s, omega_SS, rho.reshape(S, 1))
+theta = tax.replacement_rate_vals(nssmat, wss, factor_ss, e, J, omega_SS.reshape(S, 1), lambdas)
+BQss = house.get_BQ(rss, bssmat_splus1, omega_SS.reshape(S, 1), lambdas, rho.reshape(S, 1), g_n_ss)
 b_s = np.array(list(np.zeros(J).reshape((1, J))) + list(bssmat))
 taxss = tax.total_taxes(rss, b_s, wss, e, nssmat, BQss, lambdas, factor_ss, T_Hss, None, 'SS', False, parameters, theta, tau_bq)
 cssmat = house.get_cons(rss, b_s, wss, e, nssmat, BQss.reshape(1, J), lambdas.reshape(1, J), bssmat_splus1, parameters, taxss)
 
-Css = (cssmat * omega_SS).sum()
+Css = house.get_C(cssmat, omega_SS.reshape(S, 1), lambdas)
 
-diffe = Yss - (Css + Iss)
+resource_constraint = Yss - (Css + Iss)
 
-print diffe
-print diffe / Yss
+print 'Resource Constraint Difference:', resource_constraint
 
 house.constraint_checker_SS(bssmat, nssmat, cssmat, parameters)
 

Python/Development_firmproblem/SS_graphs.py

@@ -26,14 +26,18 @@
 from mpl_toolkits.mplot3d import Axes3D
 import cPickle as pickle
 
+import firm_funcs as firm
+import household_funcs as house
+
 '''
 ------------------------------------------------------------------------
     Functions
 ------------------------------------------------------------------------
 '''
 
 
-def the_inequalizer(dist, weights):
+def the_inequalizer(dist, pop_weights, ability_weights, S, J):
+    weights = np.tile(pop_weights.reshape(S, 1), (1, J)) * ability_weights.reshape(1, J)
     flattened_dist = dist.flatten()
     flattened_weights = weights.flatten()
     idx = np.argsort(flattened_dist)
@@ -84,16 +88,16 @@ def the_inequalizer(dist, weights):
 Kss_init = Kss
 Lss_init = Lss
 
-Css = (cssmat * omega_SS).sum()
-Css_init = Css
-income_init = cssmat + delta * bssmat_splus1
+Css_init = house.get_C(cssmat, omega_SS.reshape(S, 1), lambdas)
+iss_init = firm.get_I(bssmat_splus1, bssmat_splus1, delta, g_y, g_n_ss)
+income_init = cssmat + iss_init
 # print (income_init*omega_SS).sum()
 # print Css + delta * Kss
 # print Kss
 # print Lss
 # print Css_init
 # print (utility_init * omega_SS).sum()
-the_inequalizer(income_init, omega_SS)
+the_inequalizer(income_init, omega_SS, lambdas, S, J)
 
 
 
@@ -285,15 +289,16 @@ def the_inequalizer(dist, weights):
 
 
 
-Css = (cssmat * omega_SS).sum()
-income = cssmat + delta * bssmat_splus1
+Css = house.get_C(cssmat, omega_SS.reshape(S, 1), lambdas)
+iss = firm.get_I(bssmat_splus1, bssmat_splus1, delta, g_y, g_n_ss)
+income = cssmat + iss
 # print (income*omega_SS).sum()
 # print Css + delta * Kss
 # print Kss
 # print Lss
 # print Css
 # print (utility * omega_SS).sum()
-# the_inequalizer(yss, omega_SS)
+# the_inequalizer(yss, omega_SS, lambdas, S, J)
 
 print (Lss - Lss_init)/Lss_init
 

Python/Development_firmproblem/TPI.py

@@ -119,10 +119,10 @@
 income_tax_params = [a_tax_income, b_tax_income, c_tax_income, d_tax_income]
 wealth_tax_params = [h_wealth, p_wealth, m_wealth]
 ellipse_params = [b_ellipse, upsilon]
-parameters = [J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, tau_payroll, retire, mean_income_data] + income_tax_params + wealth_tax_params + ellipse_params
+parameters = [J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, g_n_ss, tau_payroll, retire, mean_income_data] + income_tax_params + wealth_tax_params + ellipse_params
 
-N_tilde = omega.sum(1).sum(1)
-omega_stationary = omega / N_tilde.reshape(T+S, 1, 1)
+N_tilde = omega.sum(1)
+omega_stationary = omega / N_tilde.reshape(T+S, 1)
 
 if get_baseline:
     initial_b = bssmat_splus1
@@ -131,16 +131,16 @@
     initial_b = bssmat_init
     initial_n = nssmat_init
 # the following needs a g_n term
-K0 = house.get_K(initial_b, omega_stationary[0])
+K0 = house.get_K(initial_b, omega_stationary[0].reshape(S, 1), lambdas, g_n_vector[0])
 b_sinit = np.array(list(np.zeros(J).reshape(1, J)) + list(initial_b[:-1]))
 b_splus1init = initial_b
-L0 = firm.get_L(e, initial_n, omega_stationary[0])
+L0 = firm.get_L(e, initial_n, omega_stationary[0].reshape(S, 1), lambdas)
 Y0 = firm.get_Y(K0, L0, parameters)
 w0 = firm.get_w(Y0, L0, parameters)
 r0 = firm.get_r(Y0, K0, parameters)
 # the following needs a g_n term
-BQ0 = house.get_BQ(r0, initial_b, omega_stationary[0], rho.reshape(S, 1))
-T_H_0 = tax.get_lump_sum(r0, b_sinit, w0, e, initial_n, BQ0, lambdas, factor_ss, omega_stationary[0], 'SS', parameters, theta, tau_bq)
+BQ0 = house.get_BQ(r0, initial_b, omega_stationary[0].reshape(S, 1), lambdas, rho.reshape(S, 1), g_n_vector[0])
+T_H_0 = tax.get_lump_sum(r0, b_sinit, w0, e, initial_n, BQ0, lambdas, factor_ss, omega_stationary[0].reshape(S, 1), 'SS', parameters, theta, tau_bq)
 tax0 = tax.total_taxes(r0, b_sinit, w0, e, initial_n, BQ0, lambdas, factor_ss, T_H_0, None, 'SS', False, parameters, theta, tau_bq)
 c0 = house.get_cons(r0, b_sinit, w0, e, initial_n, BQ0.reshape(1, J), lambdas.reshape(1, J), b_splus1init, parameters, tax0)
 
@@ -188,7 +188,7 @@ def Steady_state_TPI_solver(guesses, winit, rinit, BQinit, T_H_init, factor, j,
         Value of Euler error. (as an 2*S*J x 1 list)
     '''
 
-    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
+    J, S, T, beta, sigma, alpha, Z, delta, ltilde, nu, g_y, g_n_ss, tau_payroll, retire, mean_income_data, a_tax_income, b_tax_income, c_tax_income, d_tax_income, h_wealth, p_wealth, m_wealth, b_ellipse, upsilon = params
     length = len(guesses)/2
     b_guess = np.array(guesses[:length])
     n_guess = np.array(guesses[length:])
@@ -325,20 +325,20 @@ def Steady_state_TPI_solver(guesses, winit, rinit, BQinit, T_H_init, factor, j,
     #     print 't-loop:', euler_errors.max()
 
     b_mat[0, :, :] = initial_b
-    Kinit = (omega_stationary[:T, :, :] * b_mat[:T, :, :]).sum(2).sum(1)
-    Linit = (omega_stationary[:T, :, :] * e.reshape(
+    Kinit = (omega_stationary[:T, :].reshape(T, S, 1) * b_mat[:T, :, :] * lambdas.reshape(1, 1, J)).sum(2).sum(1) / (1.0 + g_n_vector[:T])
+    Linit = (omega_stationary[:T, :].reshape(T, S, 1) * lambdas.reshape(1, 1, J) * e.reshape(
         1, S, J) * n_mat[:T, :, :]).sum(2).sum(1)
 
     Ynew = firm.get_Y(Kinit, Linit, parameters)
     wnew = firm.get_w(Ynew, Linit, parameters)
     rnew = firm.get_r(Ynew, Kinit, parameters)
     # the following needs a g_n term
-    BQnew = (1+rnew.reshape(T, 1))*(b_mat[:T] * omega_stationary[:T] * rho.reshape(1, S, 1)).sum(1)
+    BQnew = (1+rnew.reshape(T, 1))*(b_mat[:T] * omega_stationary[:T].reshape(T, S, 1) * lambdas.reshape(1, 1, J) * rho.reshape(1, S, 1)).sum(1) / (1.0 + g_n_vector[:T].reshape(T, 1))
     bmat_s = np.zeros((T, S, J))
     bmat_s[:, 1:, :] = b_mat[:T, :-1, :]
     T_H_new = np.array(list(tax.get_lump_sum(rnew.reshape(T, 1, 1), bmat_s, wnew.reshape(
         T, 1, 1), e.reshape(1, S, J), n_mat[:T], BQnew.reshape(T, 1, J), lambdas.reshape(
-        1, 1, J), factor_ss, omega_stationary[:T], 'TPI', parameters, theta, tau_bq)) + [T_Hss]*S)
+        1, 1, J), factor_ss, omega_stationary[:T].reshape(T, S, 1), 'TPI', parameters, theta, tau_bq)) + [T_Hss]*S)
 
     winit[:T] = misc_funcs.convex_combo(wnew, winit[:T], parameters)
     rinit[:T] = misc_funcs.convex_combo(rnew, rinit[:T], parameters)
@@ -396,8 +396,8 @@ def Steady_state_TPI_solver(guesses, winit, rinit, BQinit, T_H_init, factor, j,
 
 b_mat[0, :, :] = initial_b
 
-Kpath_TPI = list(Kinit) + list(np.ones(10)*Kss)
-Lpath_TPI = list(Linit) + list(np.ones(10)*Lss)
+Kpath_TPI = np.array(list(Kinit) + list(np.ones(10)*Kss))
+Lpath_TPI = np.array(list(Linit) + list(np.ones(10)*Lss))
 BQpath_TPI = np.array(list(BQinit) + list(np.ones((10, J))*BQss))
 
 
@@ -408,12 +408,17 @@ def Steady_state_TPI_solver(guesses, winit, rinit, BQinit, T_H_init, factor, j,
 
 tax_path = tax.total_taxes(rinit[:T].reshape(T, 1, 1), b_s, winit[:T].reshape(T, 1, 1), e.reshape(
     1, S, J), n_mat[:T], BQinit[:T, :].reshape(T, 1, J), lambdas, factor_ss, T_H_init[:T].reshape(T, 1, 1), None, 'TPI', False, parameters, theta, tau_bq)
-cinit = house.get_cons(rinit[:T].reshape(T, 1, 1), b_s, winit[:T].reshape(T, 1, 1), e.reshape(1, S, J), n_mat[:T], BQinit[:T].reshape(T, 1, J), lambdas.reshape(1, 1, J), b_splus1, parameters, tax_path)
+c_path = house.get_cons(rinit[:T].reshape(T, 1, 1), b_s, winit[:T].reshape(T, 1, 1), e.reshape(1, S, J), n_mat[:T], BQinit[:T].reshape(T, 1, J), lambdas.reshape(1, 1, J), b_splus1, parameters, tax_path)
+
+Y_path = firm.get_Y(Kpath_TPI[:T], Lpath_TPI[:T], parameters)
+C_path = (c_path * omega_stationary[:T].reshape(T, S, 1) * lambdas).sum(1).sum(1)
+I_path = firm.get_I(Kpath_TPI[1:T+1], Kpath_TPI[:T], delta, g_y, g_n_vector[:T])
+print 'Resource Constraint Difference:', Y_path - C_path - I_path
 
 print'Checking time path for violations of constaints.'
 for t in xrange(T):
     house.constraint_checker_TPI(b_mat[t, :-1, :], n_mat[
-        t], cinit[t], t, parameters, N_tilde)
+        t], c_path[t], t, parameters, N_tilde)
 
 '''
 ------------------------------------------------------------------------
@@ -429,7 +434,7 @@ def Steady_state_TPI_solver(guesses, winit, rinit, BQinit, T_H_init, factor, j,
 ------------------------------------------------------------------------
 '''
 
-var_names = ['Kpath_TPI', 'b_mat', 'cinit',
+var_names = ['Kpath_TPI', 'b_mat', 'c_path',
              'eul_savings', 'eul_laborleisure', 'Lpath_TPI', 'BQpath_TPI',
              'n_mat', 'rinit', 'winit', 'Yinit', 'T_H_init', 'tax_path']
 dictionary = {}