From ae2383950292f774b6e63b482bba4b38a332b824 Mon Sep 17 00:00:00 2001 From: John Ryan Date: Thu, 26 Dec 2024 12:02:37 -0500 Subject: [PATCH 1/2] Add HSV estimates of mtrs --- iot/inverse_optimal_tax.py | 49 +++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 16 deletions(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 7288897..1dcb5d3 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -4,6 +4,7 @@ import scipy from statsmodels.nonparametric.kernel_regression import KernelReg from scipy.interpolate import UnivariateSpline +from scipy.linalg import lstsq class IOT: @@ -126,22 +127,23 @@ def compute_mtr_dist( * mtr_prime (array_like): rate of change in marginal tax rates for each income bin """ - bins = 1000 # number of equal-width bins - data.loc[:, ["z_bin"]] = pd.cut( - data[income_measure], bins, include_lowest=True - ) - binned_data = pd.DataFrame( - data[["mtr", income_measure, "z_bin", weight_var]] - .groupby(["z_bin"], observed=False) - .apply(lambda x: wm(x[["mtr", income_measure]], x[weight_var])) - ) - # make column 0 into two columns - binned_data[["mtr", income_measure]] = pd.DataFrame( - binned_data[0].tolist(), index=binned_data.index - ) - binned_data.drop(columns=0, inplace=True) - binned_data.reset_index(inplace=True) + if mtr_smoother == "kreg": + bins = 1000 # number of equal-width bins + data.loc[:, ["z_bin"]] = pd.cut( + data[income_measure], bins, include_lowest=True + ) + binned_data = pd.DataFrame( + data[["mtr", income_measure, "z_bin", weight_var]] + .groupby(["z_bin"], observed=False) + .apply(lambda x: wm(x[["mtr", income_measure]], x[weight_var])) + ) + # make column 0 into two columns + binned_data[["mtr", income_measure]] = pd.DataFrame( + binned_data[0].tolist(), index=binned_data.index + ) + binned_data.drop(columns=0, inplace=True) + binned_data.reset_index(inplace=True) mtr_function = KernelReg( binned_data["mtr"].dropna(), binned_data[income_measure].dropna(), @@ -149,10 +151,25 @@ def compute_mtr_dist( reg_type="ll", ) mtr, _ = mtr_function.fit(self.z) + mtr_prime = np.gradient(mtr, edge_order=2) + elif mtr_smoother == "HSV": + # estimate the HSV function on mtrs via weighted least squares + X = data[income_measure].values + X = np.column_stack((np.ones(len(X)), X)) + w = data[weight_var].values + w_sqrt = np.sqrt(w) + y = data["mtr"].values + X_weighted = X * w_sqrt[:, np.newaxis] + y_weighted = y * w_sqrt + coef, _, _, _ = lstsq(X_weighted, y_weighted) + tau = -coef[1] + lambda_param = np.exp(coef[0]) / (1-tau) + mtr = 1 - lambda_param * (1-tau) * self.z ** (-tau) + mtr_prime = lambda_param * tau * (1-tau) * self.z ** (-tau - 1) else: print("Please enter a value mtr_smoother method") assert False - mtr_prime = np.gradient(mtr, edge_order=2) + return mtr, mtr_prime From 90a7fd042e128225a6d08ba621af157520aa53da Mon Sep 17 00:00:00 2001 From: John Ryan Date: Mon, 30 Dec 2024 10:50:45 -0500 Subject: [PATCH 2/2] Fix HSV --- iot/inverse_optimal_tax.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 1dcb5d3..9ce432a 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -154,11 +154,11 @@ def compute_mtr_dist( mtr_prime = np.gradient(mtr, edge_order=2) elif mtr_smoother == "HSV": # estimate the HSV function on mtrs via weighted least squares - X = data[income_measure].values + X = np.log(data[income_measure].values) X = np.column_stack((np.ones(len(X)), X)) - w = data[weight_var].values + w = np.array(data[weight_var].values) w_sqrt = np.sqrt(w) - y = data["mtr"].values + y = np.log(1-data["mtr"].values) X_weighted = X * w_sqrt[:, np.newaxis] y_weighted = y * w_sqrt coef, _, _, _ = lstsq(X_weighted, y_weighted) @@ -352,7 +352,7 @@ def sw_weights(self): + ((self.theta_z * self.eti * self.mtr) / (1 - self.mtr)) + ((self.eti * self.z * self.mtr_prime) / (1 - self.mtr) ** 2) ) - integral = np.trapz(g_z, self.z) + integral = np.trapz(g_z * self.f, self.z) g_z = g_z / integral # use Lockwood and Weinzierl formula, which should be equivalent but using numerical differentiation bracket_term = ( @@ -360,11 +360,11 @@ def sw_weights(self): - self.F - (self.mtr / (1 - self.mtr)) * self.eti * self.z * self.f ) - # d_dz_bracket = np.gradient(bracket_term, edge_order=2) - d_dz_bracket = np.diff(bracket_term) / np.diff(self.z) - d_dz_bracket = np.append(d_dz_bracket, d_dz_bracket[-1]) + d_dz_bracket = np.gradient(bracket_term, edge_order=2) + # d_dz_bracket = np.diff(bracket_term) / np.diff(self.z) + # d_dz_bracket = np.append(d_dz_bracket, d_dz_bracket[-1]) g_z_numerical = -(1 / self.f) * d_dz_bracket - integral = np.trapz(g_z_numerical, self.z) + integral = np.trapz(g_z_numerical * self.f, self.z) g_z_numerical = g_z_numerical / integral return g_z, g_z_numerical