From 73a9849bc5857ab8eec218d0610885b5889482b7 Mon Sep 17 00:00:00 2001 From: Evan Magnusson Date: Thu, 16 Jul 2015 14:45:30 -0700 Subject: [PATCH] documented misc_funcs and tax_funcs --- Python/CURRENT_DEBUGGED_VERSION/misc_funcs.py | 26 +++- .../run_simulation.py | 6 +- Python/CURRENT_DEBUGGED_VERSION/tax_funcs.py | 111 ++++++++++++++++-- 3 files changed, 132 insertions(+), 11 deletions(-) diff --git a/Python/CURRENT_DEBUGGED_VERSION/misc_funcs.py b/Python/CURRENT_DEBUGGED_VERSION/misc_funcs.py index 46498da35..e2725b1b1 100644 --- a/Python/CURRENT_DEBUGGED_VERSION/misc_funcs.py +++ b/Python/CURRENT_DEBUGGED_VERSION/misc_funcs.py @@ -1,9 +1,12 @@ ''' ------------------------------------------------------------------------ -Last updated 6/19/2015 +Last updated 7/16/2015 Miscellaneous functions for SS and TPI. +This python files calls: + OUTPUT/Saved_moments/wealth_data_moments.pkl + ------------------------------------------------------------------------ ''' @@ -24,6 +27,11 @@ def perc_dif_func(simul, data): ''' Used to calculate the absolute percent difference between the data moments and model moments + Inputs: + simul = model moments (any shape) + data = data moments (same shape as simul) + Output: + output = absolute percent difference between data and model moments (same shape as simul) ''' frac = (simul - data)/data output = np.abs(frac) @@ -34,6 +42,12 @@ def convex_combo(var1, var2, params): ''' Takes the convex combination of two variables, where nu is the value between 0 and 1 in params. + Inputs: + var1 = (any shape) + var2 = (same shape as var1) + params = parameters list from model (list) (only nu is needed...perhaps it should just take that as an input) + Outputs: + combo = convex combination of var1 and var2 (same shape as var1) ''' 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 combo = nu * var1 + (1-nu)*var2 @@ -45,13 +59,23 @@ def check_wealth_calibration(wealth_model, factor_model, params): Creates a vector of the percent differences between the model and data wealth moments for the two age groups for each J group. + Inputs: + wealth_model = model wealth levels (SxJ array) + factor_model = factor to convert wealth levels to dollars (scalar) + params = parameters list from model (list) + Outputs: + wealth_fits = Fits for how well the model wealth levels match the data wealth levels ((2*J)x1 array) ''' + # Import the wealth data moments wealth_dict = pickle.load(open("OUTPUT/Saved_moments/wealth_data_moments.pkl", "r")) # Set lowest ability group's wealth to be a positive, not negative, number for the calibration wealth_dict['wealth_data_array'][2:26, 0] = 500.0 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 + # Convert wealth levels from model to dollar terms wealth_model_dollars = wealth_model * factor_model wealth_fits = np.zeros(2*J) + # Look at the percent difference between the fits for the first age group (20-44) and second age group (45-65) + # The wealth data moment indices are shifted because they start at age 18 wealth_fits[0::2] = perc_dif_func(np.mean(wealth_model_dollars[:24], axis=0), np.mean(wealth_dict['wealth_data_array'][2:26], axis=0)) wealth_fits[1::2] = perc_dif_func(np.mean(wealth_model_dollars[24:45], axis=0), np.mean(wealth_dict['wealth_data_array'][26:47], axis=0)) return wealth_fits diff --git a/Python/CURRENT_DEBUGGED_VERSION/run_simulation.py b/Python/CURRENT_DEBUGGED_VERSION/run_simulation.py index 81fa12c1c..fcd45da5e 100644 --- a/Python/CURRENT_DEBUGGED_VERSION/run_simulation.py +++ b/Python/CURRENT_DEBUGGED_VERSION/run_simulation.py @@ -148,7 +148,8 @@ d_tax_income = .219 # Wealth tax params # These are non-calibrated values, h and m just need -# need to be nonzero to avoid errors +# need to be nonzero to avoid errors. When p_wealth +# is zero, there is no wealth tax. h_wealth = 0.1 m_wealth = 1.0 p_wealth = 0.0 @@ -161,11 +162,12 @@ maxiter = 250 mindist_SS = 1e-9 mindist_TPI = 1e-6 -nu = .40 +nu = .4 flag_graphs = False get_baseline = True # Calibration parameters calibrate_model = False +# These guesses are close to the calibrated values chi_b_guess = np.array([2, 10, 90, 350, 1700, 22000, 120000]) chi_n_guess = np.array([47.12000874 , 22.22762421 , 14.34842241 , 10.67954008 , 8.41097278 , 7.15059004 , 6.46771332 , 5.85495452 , 5.46242013 , 5.00364263 diff --git a/Python/CURRENT_DEBUGGED_VERSION/tax_funcs.py b/Python/CURRENT_DEBUGGED_VERSION/tax_funcs.py index e28dab894..918030bc9 100644 --- a/Python/CURRENT_DEBUGGED_VERSION/tax_funcs.py +++ b/Python/CURRENT_DEBUGGED_VERSION/tax_funcs.py @@ -1,6 +1,6 @@ ''' ------------------------------------------------------------------------ -Last updated 6/19/2015 +Last updated 7/16/2015 Functions for taxes in SS and TPI. @@ -15,8 +15,7 @@ ------------------------------------------------------------------------ Tax functions ------------------------------------------------------------------------ - The first function gets the replacement rate values for the payroll - tax. The next 4 functions are the wealth and income tax functions, + The next 4 functions are the wealth and income tax functions, with their derivative functions. The remaining functions are used to get the total amount of taxes and lump sum taxes. ------------------------------------------------------------------------ @@ -24,7 +23,20 @@ def replacement_rate_vals(nssmat, wss, factor_ss, e, J, omega_SS, lambdas): - # Import data need to compute replacement rates, outputed from SS.py + ''' + Calculates replacement rate values for the payroll tax. + Inputs: + nssmat = labor participation rate values (SxJ array or Sx1 array) + wss = wage rate (scalar) + factor_ss = factor that converts income to dollars (scalar) + e = ability levels (SxJ array or Sx1 array) + J = number of ability types (scalar) + omega_SS = population weights by age (Sx1 array) + lambdas = ability weights (Jx1 array or scalar) + Outputs: + theta = replacement rates for each ability type (Jx1 array) + ''' + # Do a try/except, depending on whether the arrays are 1 or 2 dimensional try: AIME = ((wss * factor_ss * e * nssmat)*omega_SS).sum(0) * lambdas / 12.0 PIA = np.zeros(J) @@ -58,6 +70,14 @@ def replacement_rate_vals(nssmat, wss, factor_ss, e, J, omega_SS, lambdas): def tau_wealth(b, params): + ''' + Calculates tau_wealth based on the wealth level for an individual + Inputs: + b = wealth holdings of an individual (various length arrays or scalar) + params = parameter list of model + Outputs: + tau_w = tau_wealth (various length arrays or scalar) + ''' 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 h = h_wealth m = m_wealth @@ -67,6 +87,14 @@ def tau_wealth(b, params): def tau_w_prime(b, params): + ''' + Calculates derivative of tau_wealth based on the wealth level for an individual + Inputs: + b = wealth holdings of an individual (various length arrays or scalar) + params = parameter list of model (list) + Outputs: + tau_w_prime = derivative of tau_wealth (various length arrays or scalar) + ''' 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 h = h_wealth m = m_wealth @@ -77,8 +105,17 @@ def tau_w_prime(b, params): def tau_income(r, b, w, e, n, factor, params): ''' - Gives income tax value at a - certain income level + Gives income tax value for a certain income level + Inputs: + r = interest rate (various length list or scalar) + b = wealth holdings (various length array or scalar) + w = wage (various length list or scalar) + e = ability level (various size array or scalar) + n = labor participation rate (various length array or scalar) + factor = scaling factor (scalar) + params = parameter list of model (list) + Output: + tau = tau_income (various length array or scalar) ''' 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 a = a_tax_income @@ -95,8 +132,17 @@ def tau_income(r, b, w, e, n, factor, params): def tau_income_deriv(r, b, w, e, n, factor, params): ''' - Gives derivative of income tax value at a - certain income level + Gives derivative of income tax value at a certain income level + Inputs: + r = interest rate (various length list or scalar) + b = wealth holdings (various length array or scalar) + w = wage (various length list or scalar) + e = ability level (various size array or scalar) + n = labor participation rate (various length array or scalar) + factor = scaling factor (scalar) + params = parameter list of model (list) + Output: + tau = derivative of tau_income (various length array or scalar) ''' 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 a = a_tax_income @@ -112,6 +158,25 @@ def tau_income_deriv(r, b, w, e, n, factor, params): def get_lump_sum(r, b, w, e, n, BQ, lambdas, factor, weights, method, params, theta, tau_bq): + ''' + Gives lump sum tax value. + Inputs: + r = interest rate (various length list or scalar) + b = wealth holdings (various length array or scalar) + w = wage (various length list or scalar) + e = ability level (various size array or scalar) + n = labor participation rate (various length array or scalar) + BQ = Bequest values (various length array or scalar) + lambdas = ability levels (Jx1 array or scalar) + factor = scaling factor (scalar) + weights = population weights (various length array or scalar) + method = 'SS' or 'TPI', depending on the shape of arrays + params = parameter list of model (list) + theta = replacement rate values (Jx1 array or scalar) + tau_bq = bequest tax values (Jx1 array or scalar) + Output: + T_H = lump sum tax (Tx1 array or scalar) + ''' 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 I = r * b + w * e * n T_I = tau_income(r, b, w, e, n, factor, params) * I @@ -129,12 +194,37 @@ def get_lump_sum(r, b, w, e, n, BQ, lambdas, factor, weights, method, params, th def total_taxes(r, b, w, e, n, BQ, lambdas, factor, T_H, j, method, shift, params, theta, tau_bq): + ''' + Gives net taxes values. + Inputs: + r = interest rate (various length list or scalar) + b = wealth holdings (various length array or scalar) + w = wage (various length list or scalar) + e = ability level (various size array or scalar) + n = labor participation rate (various length array or scalar) + BQ = Bequest values (various length array or scalar) + lambdas = ability levels (Jx1 array or scalar) + factor = scaling factor (scalar) + T_H = net taxes (Tx1 array or scalar) + j = Which ability level is being computed, if doing one ability level at a time (scalar) + method = 'SS' or 'TPI' or 'TPI_scalar', depending on the shape of arrays + shift = Computing for periods 0--s or 1--(s+1) (bool) (True for 1--(s+1)) + params = parameter list of model (list) + theta = replacement rate values (Jx1 array or scalar) + tau_bq = bequest tax values (Jx1 array or scalar) + Output: + total_taxes = net taxes (various length array or scalar) + ''' 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 I = r * b + w * e * n T_I = tau_income(r, b, w, e, n, factor, params) * I T_P = tau_payroll * w * e * n T_W = tau_wealth(b, params) * b if method == 'SS': + # Depending on if we are looking at b_s or b_s+1, the + # entry for retirement will change (it shifts back one). + # The shift boolean makes sure we start replacement rates + # at the correct age. if shift is False: T_P[retire:] -= theta * w else: @@ -142,6 +232,9 @@ def total_taxes(r, b, w, e, n, BQ, lambdas, factor, T_H, j, method, shift, param T_BQ = tau_bq * BQ / lambdas elif method == 'TPI': if shift is False: + # retireTPI is different from retire, because in TPI we are counting backwards + # with different length lists. This will always be the correct location + # of retirement, depending on the shape of the lists. retireTPI = (retire - S) else: retireTPI = (retire-1 - S) @@ -152,6 +245,8 @@ def total_taxes(r, b, w, e, n, BQ, lambdas, factor, T_H, j, method, shift, param T_P[:, retire:, :] -= theta.reshape(1, 1, J) * w T_BQ = tau_bq.reshape(1, 1, J) * BQ / lambdas elif method == 'TPI_scalar': + # The above methods won't work if scalars are used. This option is only called by the + # SS_TPI_firstdoughnutring function in TPI. T_P -= theta[j] * w T_BQ = tau_bq[j] * BQ / lambdas total_taxes = T_I + T_P + T_BQ + T_W - T_H