Merge pull request PSLmodels#206 from jdebacker/testing

Representative agent model

Python/execute.py

@@ -32,7 +32,8 @@ def runner(output_base, baseline_dir, baseline=False, analytical_mtrs=True, age_
             pass
 
     if run_micro:
-        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, reform=reform, guid=guid)
+        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, age_specific=age_specific, 
+                                     start_year=user_params['start_year'], reform=reform, guid=guid)
     print ("in runner, baseline is ", baseline)
     run_params = ogusa.parameters.get_parameters(baseline=baseline, guid=guid)
     run_params['analytical_mtrs'] = analytical_mtrs
@@ -167,7 +168,8 @@ def runner_SS(output_base, baseline_dir, baseline=False, analytical_mtrs=True, a
             pass
 
     if run_micro:
-        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, reform=reform, guid=guid)
+        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, age_specific=age_specific, 
+                                     start_year=user_params['start_year'], reform=reform, guid=guid)
     print ("in runner, baseline is ", baseline)
     run_params = ogusa.parameters.get_parameters(baseline=baseline, guid=guid)
     run_params['analytical_mtrs'] = analytical_mtrs

Python/execute_large.py

@@ -33,7 +33,8 @@ def runner(output_base, baseline_dir, baseline=False, analytical_mtrs=True, age_
             pass
 
     if run_micro:
-        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, age_specific=age_specific, reform=reform, guid=guid)
+        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, age_specific=age_specific, 
+                                     start_year=user_params['start_year'], reform=reform, guid=guid)
     print ("in runner, baseline is ", baseline)
     run_params = ogusa.parameters.get_parameters(baseline=baseline, guid=guid)
     run_params['analytical_mtrs'] = analytical_mtrs
@@ -168,7 +169,8 @@ def runner_SS(output_base, baseline_dir, baseline=False, analytical_mtrs=True, a
             pass
 
     if run_micro:
-        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, age_specific=age_specific, reform=reform, guid=guid)
+        txfunc.get_tax_func_estimate(baseline=baseline, analytical_mtrs=analytical_mtrs, age_specific=age_specific, 
+                                     start_year=user_params['start_year'], reform=reform, guid=guid)
     print ("in runner, baseline is ", baseline)
     run_params = ogusa.parameters.get_parameters(baseline=baseline, guid=guid)
     run_params['analytical_mtrs'] = analytical_mtrs

Python/ogusa/SS.py

@@ -395,9 +395,14 @@ def SS_fsolve(guesses, b_guess_init, n_guess_init, chi_n, chi_b, tax_params, par
     error2 = new_r - r
     error3 = new_T_H - T_H
     error4 = new_factor/1000000 - factor/1000000
+
+    print 'mean income in model and data: ', average_income_model, mean_income_data
+    print 'model income with factor: ', average_income_model*factor
+
     print 'errors: ', error1, error2, error3, error4
     print 'T_H: ', new_T_H
     print 'factor: ', new_factor
+    print 'interest rate: ', new_r
 
     # Check and punish violations
     if r <= 0:
@@ -505,6 +510,7 @@ def SS_fsolve_reform(guesses, b_guess_init, n_guess_init, factor, chi_n, chi_b,
     print 'errors: ', error1, error2, error3
     print 'T_H: ', new_T_H
 
+
     # Check and punish violations
     if r <= 0:
         error1 += 1e9
@@ -798,8 +804,8 @@ def run_steady_state(income_tax_parameters, ss_parameters, iterative_params, bas
 
     mtrx_ss = tax.MTR_labor(rss, bssmat_s, wss, e, nssmat, factor_ss, analytical_mtrs, etr_params_3D, mtrx_params_3D)
 
-    #np.savetxt("mtr_ss_capital.csv", mtry_ss, delimiter=",")
-    #np.savetxt("mtr_ss_labor.csv", mtrx_ss, delimiter=",")
+    # np.savetxt("mtr_ss_capital.csv", mtry_ss, delimiter=",")
+    # np.savetxt("mtr_ss_labor.csv", mtrx_ss, delimiter=",")
 
     taxss_params = (J,S, retire, np.tile(np.reshape(etr_params,(S,1,etr_params.shape[1])),(1,J,1)),
                     h_wealth, p_wealth, m_wealth, tau_payroll)

Python/ogusa/get_micro_data.py

@@ -17,7 +17,7 @@
 import pickle
 
 
-def get_data(baseline=False, reform={}):
+def get_data(baseline=False, start_year=2016, reform={}):
     '''
     --------------------------------------------------------------------
     This function creates dataframes of micro data from the 
@@ -38,9 +38,9 @@ def get_data(baseline=False, reform={}):
     calc1 = Calculator(records=records1, policy=policy1)
 
     # this increment_year function extrapolates all PUF variables to the next year
-    # so this step takes the calculator to 2015
-    calc1.increment_year()
-    calc1.increment_year()
+    # so this step takes the calculator to the start_year
+    for i in range(start_year-2013):
+        calc1.increment_year()
 
     # running all the functions and calculates taxes
     calc1.calc_all()
@@ -61,10 +61,14 @@ def get_data(baseline=False, reform={}):
 
     # note that use total pension income (e01500) since don't have both the 
     # taxable (e01700) and non-taxable pension income separately
-    # don't appear to have variable for non-taxable IRS distributions
+    # don't appear to have variable for non-taxable IRA distributions
+    # capital_income_sources = ('e00300', 'e00400', 'e00600',
+    #                             'e00650', 'e01400',
+    #                             'e01500', 'e02000',
+    #                             'p22250','p23250')
     capital_income_sources = ('e00300', 'e00400', 'e00600',
                                 'e00650', 'e01400',
-                                'e01500', 'e02000',
+                                'e01700', 'e02000',
                                 'p22250','p23250')
 
     # calculating MTRs separately - can skip items with zero tax
@@ -111,7 +115,7 @@ def get_data(baseline=False, reform={}):
     # dictionary of data frames to return
     micro_data_dict = {}
 
-    micro_data_dict['2015'] = DataFrame(data = temp,
+    micro_data_dict[str(start_year)] = DataFrame(data = temp,
                       columns = ['MTR wage', 'MTR self-employed Wage', 'MTR capital income','Age',
                                  'Wage and Salaries', 'Self-Employed Income','Wage + Self-Employed Income',
                                  'Adjusted Total income', 'Total Tax Liability', 'Year', 'Weights'])
@@ -141,7 +145,7 @@ def get_data(baseline=False, reform={}):
                        columns = ['MTR wage', 'MTR self-employed Wage','MTR capital income','Age',
                                   'Wage and Salaries', 'Self-Employed Income','Wage + Self-Employed Income',
                                   'Adjusted Total income','Total Tax Liability','Year', 'Weights'])
-        print 'year: ', i
+        print 'year: ', str(calc1.current_year) 
 
     if reform:
         pkl_path = "micro_data_policy.pkl"

Python/ogusa/income.py

@@ -54,7 +54,7 @@
                      1.11e+00, -0.93939272, 1.60e+00, 1.89e+00])
 ages = np.linspace(21, 80, 60)
 ages = np.tile(ages.reshape(60, 1), (1, 7))
-income_profiles = constant + one * ages + two * ages ** 2 + three * ages ** 3
+income_profiles = constant + (one * ages) + (two * (ages ** 2)) + (three * (ages ** 3))
 income_profiles = np.exp(income_profiles)
 
 

Python/ogusa/parameters.py

@@ -24,6 +24,7 @@
 import pickle
 import txfunc
 import elliptical_u_est
+import matplotlib.pyplot as plt
 
 
 DATASET = 'REAL'
@@ -317,9 +318,24 @@ def get_full_parameters(baseline, guid, user_modifiable, metadata):
 
     mean_income_data = dict_params['tfunc_avginc'][0]
 
-    etr_params = dict_params['tfunc_etr_params_S'][:S,:BW,:]
-    mtrx_params = dict_params['tfunc_mtrx_params_S'][:S,:BW,:]
-    mtry_params = dict_params['tfunc_mtry_params_S'][:S,:BW,:]
+    # etr_params = dict_params['tfunc_etr_params_S'][:S,:BW,:]
+    # mtrx_params = dict_params['tfunc_mtrx_params_S'][:S,:BW,:]
+    # mtry_params = dict_params['tfunc_mtry_params_S'][:S,:BW,:]
+
+    # set etrs and mtrs to constant rates over income/age
+    etr_params = np.zeros((S,BW,10))
+    mtrx_params = np.zeros((S,BW,10))
+    mtry_params = np.zeros((S,BW,10))
+    etr_params[:,:,7] = dict_params['tfunc_avg_etr']
+    mtrx_params[:,:,7] = dict_params['tfunc_avg_mtrx']
+    mtry_params[:,:,7] = dict_params['tfunc_avg_mtry']
+    etr_params[:,:,9] = dict_params['tfunc_avg_etr']
+    mtrx_params[:,:,9] = dict_params['tfunc_avg_mtrx']
+    mtry_params[:,:,9] = dict_params['tfunc_avg_mtry']
+    etr_params[:,:,5] = 1.0
+    mtrx_params[:,:,5] = 1.0
+    mtry_params[:,:,5] = 1.0
+
 
     # To zero out income taxes, uncomment the following 3 lines:
     # etr_params[:,:,6:] = 0.0
@@ -350,9 +366,11 @@ def get_full_parameters(baseline, guid, user_modifiable, metadata):
     flag_graphs = False
     #   Calibration parameters
     # These guesses are close to the calibrated values
+    chi_b_guess = np.ones((J,)) * 80.0
     #chi_b_guess = np.array([0.7, 0.7, 1.0, 1.2, 1.2, 1.2, 1.4])
-    #chi_b_guess = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
-    chi_b_guess = np.array([5, 10, 90, 250, 250, 250, 250])
+    #chi_b_guess = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 4.0, 10.0])
+    #chi_b_guess = np.array([5, 10, 90, 250, 250, 250, 250])
+    #chi_b_guess = np.array([2, 10, 90, 350, 1700, 22000, 120000])
     chi_n_guess = np.array([38.12000874, 33.22762421, 25.34842241, 26.67954008, 24.41097278, 
                             23.15059004, 22.46771332, 21.85495452, 21.46242013, 22.00364263, 
                             21.57322063, 21.53371545, 21.29828515, 21.10144524, 20.8617942, 
@@ -395,7 +413,10 @@ def get_full_parameters(baseline, guid, user_modifiable, metadata):
     g_n_vector = np.tile(g_n_ss,(T+S,))
 
 
-    e = get_e(S, J, starting_age, ending_age, lambdas, omega_SS, flag_graphs)
+    e_hetero = get_e(S, J, starting_age, ending_age, lambdas, omega_SS, flag_graphs)
+    e = np.tile(((e_hetero*lambdas).sum(axis=1)).reshape(S,1),(1,J))
+    e /= (e * omega_SS.reshape(S, 1)* lambdas.reshape(1, J)).sum()
+
 
     allvars = dict(locals())
 

Python/ogusa/tax.py

@@ -293,6 +293,7 @@ def MTR_capital(r, b, w, e, n, factor, analytical_mtrs, etr_params, mtry_params)
         num = (A*(x**2)) + (B*(y**2)) + (C*x*y) + (D*x) + (E*y)
         denom = (A*(x**2)) + (B*(y**2)) + (C*x*y) + (D*x) + (E*y) + F
         tau =  (Phi*(num/denom)) + K
+
     return tau
 
 

Python/ogusa/txfunc.py

@@ -427,16 +427,16 @@ def find_outliers(sse_mat, age_vec, se_mult, varstr, graph=False):
         # Plot sum of squared errors of tax functions over age for each
         # year of budget window
         fig, ax = plt.subplots()
-        plt.plot(age_vec, sse_mat[:,0], label='2015')
-        plt.plot(age_vec, sse_mat[:,1], label='2016')
-        plt.plot(age_vec, sse_mat[:,2], label='2017')
-        plt.plot(age_vec, sse_mat[:,3], label='2018')
-        plt.plot(age_vec, sse_mat[:,4], label='2019')
-        plt.plot(age_vec, sse_mat[:,5], label='2020')
-        plt.plot(age_vec, sse_mat[:,6], label='2021')
-        plt.plot(age_vec, sse_mat[:,7], label='2022')
-        plt.plot(age_vec, sse_mat[:,8], label='2023')
-        plt.plot(age_vec, sse_mat[:,9], label='2024')
+        plt.plot(age_vec, sse_mat[:,0], label=str(start_year))
+        plt.plot(age_vec, sse_mat[:,1], label=str(start_year+1))
+        plt.plot(age_vec, sse_mat[:,2], label=str(start_year+2))
+        plt.plot(age_vec, sse_mat[:,3], label=str(start_year+3))
+        plt.plot(age_vec, sse_mat[:,4], label=str(start_year+4))
+        plt.plot(age_vec, sse_mat[:,5], label=str(start_year+5))
+        plt.plot(age_vec, sse_mat[:,6], label=str(start_year+6))
+        plt.plot(age_vec, sse_mat[:,7], label=str(start_year+7))
+        plt.plot(age_vec, sse_mat[:,8], label=str(start_year+8))
+        plt.plot(age_vec, sse_mat[:,9], label=str(start_year+9))
         # for the minor ticks, use no labels; default NullFormatter
         minorLocator   = MultipleLocator(1)
         ax.xaxis.set_minor_locator(minorLocator)
@@ -471,16 +471,16 @@ def find_outliers(sse_mat, age_vec, se_mult, varstr, graph=False):
             # Plot sum of squared errors of tax functions over age for
             # each year of budget window
             fig, ax = plt.subplots()
-            plt.plot(age_vec, sse_mat_new[:,0], label='2015')
-            plt.plot(age_vec, sse_mat_new[:,1], label='2016')
-            plt.plot(age_vec, sse_mat_new[:,2], label='2017')
-            plt.plot(age_vec, sse_mat_new[:,3], label='2018')
-            plt.plot(age_vec, sse_mat_new[:,4], label='2019')
-            plt.plot(age_vec, sse_mat_new[:,5], label='2020')
-            plt.plot(age_vec, sse_mat_new[:,6], label='2021')
-            plt.plot(age_vec, sse_mat_new[:,7], label='2022')
-            plt.plot(age_vec, sse_mat_new[:,8], label='2023')
-            plt.plot(age_vec, sse_mat_new[:,9], label='2024')
+            plt.plot(age_vec, sse_mat_new[:,0], label=str(start_year))
+            plt.plot(age_vec, sse_mat_new[:,1], label=str(start_year+1))
+            plt.plot(age_vec, sse_mat_new[:,2], label=str(start_year+2))
+            plt.plot(age_vec, sse_mat_new[:,3], label=str(start_year+3))
+            plt.plot(age_vec, sse_mat_new[:,4], label=str(start_year+4))
+            plt.plot(age_vec, sse_mat_new[:,5], label=str(start_year+5))
+            plt.plot(age_vec, sse_mat_new[:,6], label=str(start_year+6))
+            plt.plot(age_vec, sse_mat_new[:,7], label=str(start_year+7))
+            plt.plot(age_vec, sse_mat_new[:,8], label=str(start_year+8))
+            plt.plot(age_vec, sse_mat_new[:,9], label=str(start_year+9))
             # for the minor ticks, use no labels; default NullFormatter
             minorLocator   = MultipleLocator(1)
             ax.xaxis.set_minor_locator(minorLocator)
@@ -503,16 +503,16 @@ def find_outliers(sse_mat, age_vec, se_mult, varstr, graph=False):
                 # Plot sum of squared errors of tax functions over age
                 # for each year of budget window
                 fig, ax = plt.subplots()
-                plt.plot(age_vec, sse_mat_new2[:,0], label='2015')
-                plt.plot(age_vec, sse_mat_new2[:,1], label='2016')
-                plt.plot(age_vec, sse_mat_new2[:,2], label='2017')
-                plt.plot(age_vec, sse_mat_new2[:,3], label='2018')
-                plt.plot(age_vec, sse_mat_new2[:,4], label='2019')
-                plt.plot(age_vec, sse_mat_new2[:,5], label='2020')
-                plt.plot(age_vec, sse_mat_new2[:,6], label='2021')
-                plt.plot(age_vec, sse_mat_new2[:,7], label='2022')
-                plt.plot(age_vec, sse_mat_new2[:,8], label='2023')
-                plt.plot(age_vec, sse_mat_new2[:,9], label='2024')
+                plt.plot(age_vec, sse_mat_new2[:,0], label=str(start_year))
+                plt.plot(age_vec, sse_mat_new2[:,1], label=str(start_year+1))
+                plt.plot(age_vec, sse_mat_new2[:,2], label=str(start_year+2))
+                plt.plot(age_vec, sse_mat_new2[:,3], label=str(start_year+3))
+                plt.plot(age_vec, sse_mat_new2[:,4], label=str(start_year+4))
+                plt.plot(age_vec, sse_mat_new2[:,5], label=str(start_year+5))
+                plt.plot(age_vec, sse_mat_new2[:,6], label=str(start_year+6))
+                plt.plot(age_vec, sse_mat_new2[:,7], label=str(start_year+7))
+                plt.plot(age_vec, sse_mat_new2[:,8], label=str(start_year+8))
+                plt.plot(age_vec, sse_mat_new2[:,9], label=str(start_year+9))
                 # for the minor ticks, use no labels; default NullFormatter
                 minorLocator   = MultipleLocator(1)
                 ax.xaxis.set_minor_locator(minorLocator)
@@ -586,7 +586,7 @@ def replace_outliers(param_arr, sse_big_mat,):
 
 
 
-def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,reform={}):
+def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False, start_year=2016,reform={}):
     '''
     --------------------------------------------------------------------
     This function estimates functions for the ETR, MTR on Labor Income,
@@ -649,6 +649,9 @@ def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,r
     mtry_sumsq_arr = np.zeros((s_max - s_min + 1, tpers))
     mtry_obs_arr = np.zeros((s_max - s_min + 1, tpers))
     AvgInc = np.zeros(tpers)
+    AvgETR = np.zeros(tpers)
+    AvgMTRx = np.zeros(tpers)
+    AvgMTRy = np.zeros(tpers)
     TotPop_yr = np.zeros(tpers)
     PopPct_age = np.zeros((s_max-s_min+1, tpers))
     desc_data = False
@@ -659,7 +662,7 @@ def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,r
     # cmap1 = matplotlib.cm.get_cmap('jet')
     # cmap1 = matplotlib.cm.get_cmap('coolwarm')
 
-    beg_yr = int(2015)
+    beg_yr = int(start_year)
     end_yr = int(beg_yr+BW-1)
     years_list = np.arange(beg_yr, end_yr + 1)
     start_time = time.clock()
@@ -671,7 +674,7 @@ def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,r
     '''
 
     # call tax caculator and get microdata
-    micro_data = get_micro_data.get_data(baseline=baseline, reform=reform)
+    micro_data = get_micro_data.get_data(baseline=baseline, start_year=beg_yr, reform=reform)
     """if reform:
         micro_data = pickle.load( open( "micro_data_policy.pkl", "rb" ) )
     else:
@@ -724,6 +727,21 @@ def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,r
             (((data['Adjusted Total income'] * data['Weights']).sum())
             / data['Weights'].sum())
 
+        # Calculate average ETR and MTRs (weight by population weights
+        #    and income) for each year
+        AvgETR[t-beg_yr] = \
+            (((data['Effective Tax Rate']*data['Adjusted Total income'] * data['Weights']).sum())
+            / (data['Adjusted Total income']*data['Weights']).sum())
+
+        AvgMTRx[t-beg_yr] = \
+            (((data['MTR Labor']*data['Adjusted Total income'] * data['Weights']).sum())
+            / (data['Adjusted Total income']*data['Weights']).sum())
+
+        AvgMTRy[t-beg_yr] = \
+            (((data['MTR capital income']*data['Adjusted Total income'] * data['Weights']).sum())
+            / (data['Adjusted Total income']*data['Weights']).sum())
+
+
         # Calculate total population in each year
         TotPop_yr[t-beg_yr] = data['Weights'].sum()
 
@@ -1696,6 +1714,7 @@ def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,r
     # Save tax function parameters array and computation time in pickle
     dict_params = dict([('tfunc_etr_params_S', etrparam_arr_S),
         ('tfunc_mtrx_params_S', mtrxparam_arr_S), ('tfunc_mtry_params_S', mtryparam_arr_S), ('tfunc_avginc', AvgInc),
+        ('tfunc_avg_etr', AvgETR), ('tfunc_avg_mtrx', AvgMTRx), ('tfunc_avg_mtry', AvgMTRy),
         ('tfunc_etr_sumsq', etr_sumsq_arr), ('tfunc_mtrx_sumsq', mtrx_sumsq_arr), ('tfunc_mtry_sumsq', mtry_sumsq_arr),
         ('tfunc_etr_obs', etr_obs_arr), ('tfunc_mtrx_obs', mtrx_obs_arr), ('tfunc_mtry_obs', mtry_obs_arr),
         ('tfunc_time', elapsed_time)])
@@ -1705,9 +1724,9 @@ def tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False,r
     return dict_params
 
 
-def get_tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False, reform={}, guid=''):
+def get_tax_func_estimate(baseline=False, analytical_mtrs=True, age_specific=False, start_year=2016, reform={}, guid=''):
     # Code to run manually from here:
-    dict_params = tax_func_estimate(baseline, analytical_mtrs, age_specific, reform)
+    dict_params = tax_func_estimate(baseline, analytical_mtrs, age_specific, start_year, reform)
     if baseline:
         baseline_pckl = "TxFuncEst_baseline{}.pkl".format(guid)
         pkl_path = os.path.join(TAX_ESTIMATE_PATH, baseline_pckl)