From 003d214d731ac399262edccb398f6263b015ed1c Mon Sep 17 00:00:00 2001 From: jdebacker Date: Thu, 26 Oct 2023 16:51:28 -0400 Subject: [PATCH 1/9] change ic_elast to eti --- examples/ComparingCandidatePlatforms.ipynb | 2 +- examples/example.py | 2 +- iot/inverse_optimal_tax.py | 12 ++++++------ iot/iot_user.py | 14 +++++++------- iot/tests/test_inverse_optimal_tax.py | 4 ++-- 5 files changed, 17 insertions(+), 17 deletions(-) diff --git a/examples/ComparingCandidatePlatforms.ipynb b/examples/ComparingCandidatePlatforms.ipynb index 195334a..01bd9af 100644 --- a/examples/ComparingCandidatePlatforms.ipynb +++ b/examples/ComparingCandidatePlatforms.ipynb @@ -110,7 +110,7 @@ " mtr_wrt=\"e00200p\",\n", " income_measure=income_measure,\n", " weight_var=\"s006\",\n", - " inc_elast=0.25,\n", + " eti=0.25,\n", ")" ] }, diff --git a/examples/example.py b/examples/example.py index 061bc15..40ba9ec 100644 --- a/examples/example.py +++ b/examples/example.py @@ -21,7 +21,7 @@ baseline_policies=[None, None], labels=["2017 Law", "Biden 2020"], years=[2017, 2020], - inc_elast=0.2, + eti=0.2, ) # %% diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 20849ec..c74bd5f 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -19,7 +19,7 @@ class IOT: weight_var, mtr income_measure (str): name of income measure from data to use weight_var (str): name of weight measure from data to use - inc_elast (scalar): compensated elasticity of taxable income + eti (scalar): compensated elasticity of taxable income w.r.t. the marginal tax rate bandwidth (scalar): size of income bins in units of income lower_bound (scalar): minimum income to consider @@ -38,7 +38,7 @@ def __init__( data, income_measure="e00200", weight_var="s006", - inc_elast=0.25, + eti=0.25, bandwidth=1000, lower_bound=0, upper_bound=500000, @@ -54,7 +54,7 @@ def __init__( # (data[income_measure] >= lower_bound) # & (data[income_measure] <= upper_bound) # ] - self.inc_elast = inc_elast + self.eti = eti self.z, self.F, self.f, self.f_prime = self.compute_income_dist( data, income_measure, weight_var, dist_type, kde_bw ) @@ -243,9 +243,9 @@ def sw_weights(self): """ g_z = ( 1 - + ((self.theta_z * self.inc_elast * self.mtr) / (1 - self.mtr)) + + ((self.theta_z * self.eti * self.mtr) / (1 - self.mtr)) + ( - (self.inc_elast * self.z * self.mtr_prime) + (self.eti * self.z * self.mtr_prime) / (1 - self.mtr) ** 2 ) ) @@ -253,7 +253,7 @@ def sw_weights(self): bracket_term = ( 1 - self.F - - (self.mtr / (1 - self.mtr)) * self.inc_elast * self.z * 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) diff --git a/iot/iot_user.py b/iot/iot_user.py index 046f950..c5a12d6 100644 --- a/iot/iot_user.py +++ b/iot/iot_user.py @@ -30,7 +30,7 @@ class iot_comparison: mtr_wtr (str): name of income source to compute MTR on income_measure (str): name of income measure from data to use weight_var (str): name of weight measure from data to use - inc_elast (scalar): compensated elasticity of taxable income + eti (scalar): compensated elasticity of taxable income w.r.t. the marginal tax rate bandwidth (scalar): size of income bins in units of income lower_bound (scalar): minimum income to consider @@ -55,7 +55,7 @@ def __init__( mtr_wrt="e00200p", income_measure="e00200", weight_var="s006", - inc_elast=0.25, + eti=0.25, bandwidth=1000, lower_bound=0, upper_bound=500000, @@ -103,7 +103,7 @@ def __init__( j, income_measure=income_measure, weight_var=weight_var, - inc_elast=inc_elast, + eti=eti, bandwidth=bandwidth, lower_bound=lower_bound, upper_bound=upper_bound, @@ -208,15 +208,15 @@ def JJZFig4(self, policy="Current Law"): # g1 with mtr_prime = 0 g1 = ( 0 - + ((df.theta_z * self.iot[k].inc_elast * df.mtr) / (1 - df.mtr)) - + ((self.iot[k].inc_elast * df.z * 0) / (1 - df.mtr) ** 2) + + ((df.theta_z * self.iot[k].eti * df.mtr) / (1 - df.mtr)) + + ((self.iot[k].eti * df.z * 0) / (1 - df.mtr) ** 2) ) # g2 with theta_z = 0 g2 = ( 0 - + ((0 * self.iot[k].inc_elast * df.mtr) / (1 - df.mtr)) + + ((0 * self.iot[k].eti * df.mtr) / (1 - df.mtr)) + ( - (self.iot[k].inc_elast * df.z * df.mtr_prime) + (self.iot[k].eti * df.z * df.mtr_prime) / (1 - df.mtr) ** 2 ) ) diff --git a/iot/tests/test_inverse_optimal_tax.py b/iot/tests/test_inverse_optimal_tax.py index df8d823..53d71e1 100644 --- a/iot/tests/test_inverse_optimal_tax.py +++ b/iot/tests/test_inverse_optimal_tax.py @@ -157,7 +157,7 @@ def test_IOT_df(): # ) # iot2 = copy.deepcopy(iot1) # iot2.theta_z = np.array([1.7, 2.4, 99.0, 1.5, 1.5, 1.5]) -# iot2.inc_elast = np.array([0.3, 0.1, 0.0, 0.4, 0.4, 0.4]) +# iot2.eti = np.array([0.3, 0.1, 0.0, 0.4, 0.4, 0.4]) # iot2.mtr = np.array([0.25, 0.2, 0.25, 0.25, 0.25, 0.0]) # iot2.z = np.array([5000.0, 5000.0, 5000.0, 5000.0, 300.0, 300.0]) # iot2.mtr_prime = np.array([0.03, 0.03, 0.03, 0.0, 0.0, 0.0]) @@ -285,7 +285,7 @@ def test_IOT_df(): # weight_var=weight_var, # dist_type=dist_type, # mtr_smoother=mtr_smoother, -# inc_elast=elasticity, +# eti=elasticity, # ) # if sim_dist_type == "exponential": # g_z_test = iot_test.g_z From 8e0deae2601c7d0bc15355687dda56ff5e7580e4 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Thu, 26 Oct 2023 22:35:40 -0400 Subject: [PATCH 2/9] add ability to enter varying eti --- iot/inverse_optimal_tax.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index c74bd5f..28439c9 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -3,7 +3,7 @@ import scipy.stats as st import scipy from statsmodels.nonparametric.kernel_regression import KernelReg - +from scipy.interpolate import UnivariateSpline class IOT: """ @@ -54,14 +54,29 @@ def __init__( # (data[income_measure] >= lower_bound) # & (data[income_measure] <= upper_bound) # ] - self.eti = eti + # Get income distribution self.z, self.F, self.f, self.f_prime = self.compute_income_dist( data, income_measure, weight_var, dist_type, kde_bw ) + # see if eti is a scalar + if isinstance(eti, float): + self.eti = eti + else: # if not, then it should be a dict with keys containing lists as values + # check that same number of ETI values as knot points + assert len(eti["knot_points"]) == len(eti["eti_values"]) + # want to interpolate across income distribution with knot points + # assume that eti can't go beyond 1 (or the max of the eti_values provided) + eti_spl = UnivariateSpline( + eti["knot_points"], eti["eti_values"], k=1, s=0 + ) + self.eti = eti_spl(self.z) + # compute marginal tax rate schedule self.mtr, self.mtr_prime = self.compute_mtr_dist( data, weight_var, income_measure, mtr_smoother, mtr_smooth_param ) + # compute theta_z, the elasticity of the tax base self.theta_z = 1 + ((self.z * self.f_prime) / self.f) + # compute the social welfare weights self.g_z, self.g_z_numerical = self.sw_weights() def df(self): From fba48593a94d6a2326645ab46f6652a8de5e42e6 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Thu, 26 Oct 2023 22:36:05 -0400 Subject: [PATCH 3/9] black format --- iot/inverse_optimal_tax.py | 8 +++----- iot/iot_user.py | 5 +---- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 28439c9..1c7b387 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -5,6 +5,7 @@ from statsmodels.nonparametric.kernel_regression import KernelReg from scipy.interpolate import UnivariateSpline + class IOT: """ Constructor for the IOT class. @@ -68,7 +69,7 @@ def __init__( # assume that eti can't go beyond 1 (or the max of the eti_values provided) eti_spl = UnivariateSpline( eti["knot_points"], eti["eti_values"], k=1, s=0 - ) + ) self.eti = eti_spl(self.z) # compute marginal tax rate schedule self.mtr, self.mtr_prime = self.compute_mtr_dist( @@ -259,10 +260,7 @@ def sw_weights(self): g_z = ( 1 + ((self.theta_z * self.eti * self.mtr) / (1 - self.mtr)) - + ( - (self.eti * self.z * self.mtr_prime) - / (1 - self.mtr) ** 2 - ) + + ((self.eti * self.z * self.mtr_prime) / (1 - self.mtr) ** 2) ) # use Lockwood and Weinzierl formula, which should be equivalent but using numerical differentiation bracket_term = ( diff --git a/iot/iot_user.py b/iot/iot_user.py index c5a12d6..cc652f6 100644 --- a/iot/iot_user.py +++ b/iot/iot_user.py @@ -215,10 +215,7 @@ def JJZFig4(self, policy="Current Law"): g2 = ( 0 + ((0 * self.iot[k].eti * df.mtr) / (1 - df.mtr)) - + ( - (self.iot[k].eti * df.z * df.mtr_prime) - / (1 - df.mtr) ** 2 - ) + + ((self.iot[k].eti * df.z * df.mtr_prime) / (1 - df.mtr) ** 2) ) plot_df = pd.DataFrame( { From 9ade6e9dfcfdb72fbea3fe7c1ffb7cd4adbfa758 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Fri, 27 Oct 2023 14:24:33 -0400 Subject: [PATCH 4/9] add function to back out implies ETI --- iot/inverse_optimal_tax.py | 42 +++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 1c7b387..98d403e 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -244,7 +244,8 @@ def sw_weights(self): r""" Returns the social welfare weights for a given tax policy. - See Jacobs, Jongen, and Zoutman (2017) + See Jacobs, Jongen, and Zoutman (2017) and + Lockwood and Weinzierl (2016) for details. .. math:: g_{z} = 1 + \theta_z \varepsilon^{c}\frac{T'(z)}{(1-T'(z))} + @@ -275,6 +276,45 @@ def sw_weights(self): return g_z, g_z_numerical +def find_eti(iot1, iot2, g_z_type="g_z_numerical"): + """ + This function solves for the ETI that would result in the + policy represented via MTRs in iot2 be consistent with the + social welfare function inferred from the policies of iot1. + + .. math:: + \varepilon_{z} = \frac{(1-T'(z))}{T'(z)}\frac{(1-F(z))}{zf(z)}\int_{z}^{\infty}\frac{1-g_{\tilde{z}}{1-F(y)}dF(\tilde{z}) + + Args: + iot1 (IOT): IOT class instance representing baseline policy + iot2 (IOT): IOT class instance representing reform policy + g_z_type (str): type of social welfare function to use + Options are: + * 'g_z' for the analytical formula + * 'g_z_numerical' for the numerical approximation + + Returns: + eti_beliefs (array-like): vector of ETI beliefs over z + """ + if g_z_type == "g_z": + g_z = iot1.g_z + else: + g_z = iot1.g_z_numerical + # The equation below is a simplication of the above to make the integration easier + eti_beliefs_lw = ( + ((1 - iot2.mtr) / iot2.mtr) + * ((1 - iot2.F) / (iot2.z * iot2.f)) + * (1 - iot2.F - (g_z.sum() - np.cumsum(g_z))) + ) + # derivation from JJZ analytical solution that doesn't involved integration + eti_beliefs_jjz = (g_z - 1) / ( + (iot2.theta_z * (iot2.mtr / (1 - iot2.mtr))) + + (iot2.z * (iot2.mtr_prime / (1 - iot2.mtr) ** 2)) + ) + + return eti_beliefs_lw, eti_beliefs_jjz + + def wm(value, weight): """ Weighted mean function that allows for zero division From 44547a0ed7342208ddcc04cd9ddd5866cb15ca51 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Sat, 28 Oct 2023 16:44:34 -0400 Subject: [PATCH 5/9] use cubic spline to interpolate ETI(z) if enough data --- iot/inverse_optimal_tax.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 98d403e..58f3c3b 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -67,8 +67,12 @@ def __init__( assert len(eti["knot_points"]) == len(eti["eti_values"]) # want to interpolate across income distribution with knot points # assume that eti can't go beyond 1 (or the max of the eti_values provided) + if len(eti["knot_points"]) > 3: + spline_order = 3 + else: + spline_order = 1 eti_spl = UnivariateSpline( - eti["knot_points"], eti["eti_values"], k=1, s=0 + eti["knot_points"], eti["eti_values"], k=spline_order, s=0 ) self.eti = eti_spl(self.z) # compute marginal tax rate schedule @@ -276,7 +280,7 @@ def sw_weights(self): return g_z, g_z_numerical -def find_eti(iot1, iot2, g_z_type="g_z_numerical"): +def find_eti(iot1, iot2, g_z_type="g_z"): """ This function solves for the ETI that would result in the policy represented via MTRs in iot2 be consistent with the @@ -302,9 +306,8 @@ def find_eti(iot1, iot2, g_z_type="g_z_numerical"): g_z = iot1.g_z_numerical # The equation below is a simplication of the above to make the integration easier eti_beliefs_lw = ( - ((1 - iot2.mtr) / iot2.mtr) - * ((1 - iot2.F) / (iot2.z * iot2.f)) - * (1 - iot2.F - (g_z.sum() - np.cumsum(g_z))) + ((1 - iot2.mtr) / (iot2.z * iot2.f * iot2.mtr)) * + (1 - iot2.F - (g_z.sum() - np.cumsum(g_z))) ) # derivation from JJZ analytical solution that doesn't involved integration eti_beliefs_jjz = (g_z - 1) / ( From f069f4df9e8bde64b30c83f1fd3abf53681a891b Mon Sep 17 00:00:00 2001 From: jdebacker Date: Sat, 28 Oct 2023 21:13:04 -0400 Subject: [PATCH 6/9] format --- iot/inverse_optimal_tax.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 58f3c3b..05fe292 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -305,9 +305,8 @@ def find_eti(iot1, iot2, g_z_type="g_z"): else: g_z = iot1.g_z_numerical # The equation below is a simplication of the above to make the integration easier - eti_beliefs_lw = ( - ((1 - iot2.mtr) / (iot2.z * iot2.f * iot2.mtr)) * - (1 - iot2.F - (g_z.sum() - np.cumsum(g_z))) + eti_beliefs_lw = ((1 - iot2.mtr) / (iot2.z * iot2.f * iot2.mtr)) * ( + 1 - iot2.F - (g_z.sum() - np.cumsum(g_z)) ) # derivation from JJZ analytical solution that doesn't involved integration eti_beliefs_jjz = (g_z - 1) / ( From 3f4a627566004c89eedb39ebc21dcd950846485a Mon Sep 17 00:00:00 2001 From: jdebacker Date: Tue, 21 Nov 2023 15:27:41 -0500 Subject: [PATCH 7/9] fix typo --- iot/inverse_optimal_tax.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 05fe292..e5aff57 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -287,7 +287,7 @@ def find_eti(iot1, iot2, g_z_type="g_z"): social welfare function inferred from the policies of iot1. .. math:: - \varepilon_{z} = \frac{(1-T'(z))}{T'(z)}\frac{(1-F(z))}{zf(z)}\int_{z}^{\infty}\frac{1-g_{\tilde{z}}{1-F(y)}dF(\tilde{z}) + \varepsilon_{z} = \frac{(1-T'(z))}{T'(z)}\frac{(1-F(z))}{zf(z)}\int_{z}^{\infty}\frac{1-g_{\tilde{z}}{1-F(y)}dF(\tilde{z}) Args: iot1 (IOT): IOT class instance representing baseline policy From 944ebae1f58ce66c0af2ac944da5f3ec95454cb9 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Tue, 21 Nov 2023 15:27:59 -0500 Subject: [PATCH 8/9] format --- iot/generate_data.py | 4 +--- iot/inverse_optimal_tax.py | 24 +++++++----------------- 2 files changed, 8 insertions(+), 20 deletions(-) diff --git a/iot/generate_data.py b/iot/generate_data.py index b0d2df6..27c00d0 100644 --- a/iot/generate_data.py +++ b/iot/generate_data.py @@ -63,9 +63,7 @@ def gen_microdata( baseline = tc.Policy.read_json_reform(s) else: baseline = s - pol1.implement_reform( - baseline, print_warnings=False, raise_errors=False - ) + pol1.implement_reform(baseline, print_warnings=False, raise_errors=False) else: pol1 = tc.Policy() diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index e5aff57..32fbe4c 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -130,9 +130,7 @@ def compute_mtr_dist( for each income bin """ bins = 1000 # number of equal-width bins - data.loc[:, ["z_bin"]] = pd.cut( - data[income_measure], bins, include_lowest=True - ) + 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"]) @@ -191,14 +189,11 @@ def compute_income_dist( # drop zero income observations data = data[data[income_measure] > 0] if dist_type == "log_normal": - mu = ( - np.log(data[income_measure]) * data[weight_var] - ).sum() / data[weight_var].sum() + mu = (np.log(data[income_measure]) * data[weight_var]).sum() / data[ + weight_var + ].sum() sigmasq = ( - ( - ((np.log(data[income_measure]) - mu) ** 2) - * data[weight_var] - ).values + (((np.log(data[income_measure]) - mu) ** 2) * data[weight_var]).values / data[weight_var].sum() ).sum() # F = st.lognorm.cdf(z_line, s=(sigmasq) ** 0.5, scale=np.exp(mu)) @@ -209,10 +204,7 @@ def compute_income_dist( # analytical derivative of lognormal sigma = np.sqrt(sigmasq) F = (1 / 2) * ( - 1 - + scipy.special.erf( - (np.log(z_line) - mu) / (np.sqrt(2) * sigma) - ) + 1 + scipy.special.erf((np.log(z_line) - mu) / (np.sqrt(2) * sigma)) ) f = ( (1 / (sigma * np.sqrt(2 * np.pi))) @@ -269,9 +261,7 @@ def sw_weights(self): ) # use Lockwood and Weinzierl formula, which should be equivalent but using numerical differentiation bracket_term = ( - 1 - - self.F - - (self.mtr / (1 - self.mtr)) * self.eti * self.z * self.f + 1 - 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) From 50efedffda9b5fbd2ecc9e91cd17936e9fca2fe3 Mon Sep 17 00:00:00 2001 From: jdebacker Date: Wed, 22 Nov 2023 11:36:58 -0500 Subject: [PATCH 9/9] format --- iot/generate_data.py | 4 +++- iot/inverse_optimal_tax.py | 24 +++++++++++++++++------- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/iot/generate_data.py b/iot/generate_data.py index 27c00d0..b0d2df6 100644 --- a/iot/generate_data.py +++ b/iot/generate_data.py @@ -63,7 +63,9 @@ def gen_microdata( baseline = tc.Policy.read_json_reform(s) else: baseline = s - pol1.implement_reform(baseline, print_warnings=False, raise_errors=False) + pol1.implement_reform( + baseline, print_warnings=False, raise_errors=False + ) else: pol1 = tc.Policy() diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py index 32fbe4c..e5aff57 100644 --- a/iot/inverse_optimal_tax.py +++ b/iot/inverse_optimal_tax.py @@ -130,7 +130,9 @@ def compute_mtr_dist( for each income bin """ bins = 1000 # number of equal-width bins - data.loc[:, ["z_bin"]] = pd.cut(data[income_measure], bins, include_lowest=True) + 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"]) @@ -189,11 +191,14 @@ def compute_income_dist( # drop zero income observations data = data[data[income_measure] > 0] if dist_type == "log_normal": - mu = (np.log(data[income_measure]) * data[weight_var]).sum() / data[ - weight_var - ].sum() + mu = ( + np.log(data[income_measure]) * data[weight_var] + ).sum() / data[weight_var].sum() sigmasq = ( - (((np.log(data[income_measure]) - mu) ** 2) * data[weight_var]).values + ( + ((np.log(data[income_measure]) - mu) ** 2) + * data[weight_var] + ).values / data[weight_var].sum() ).sum() # F = st.lognorm.cdf(z_line, s=(sigmasq) ** 0.5, scale=np.exp(mu)) @@ -204,7 +209,10 @@ def compute_income_dist( # analytical derivative of lognormal sigma = np.sqrt(sigmasq) F = (1 / 2) * ( - 1 + scipy.special.erf((np.log(z_line) - mu) / (np.sqrt(2) * sigma)) + 1 + + scipy.special.erf( + (np.log(z_line) - mu) / (np.sqrt(2) * sigma) + ) ) f = ( (1 / (sigma * np.sqrt(2 * np.pi))) @@ -261,7 +269,9 @@ def sw_weights(self): ) # use Lockwood and Weinzierl formula, which should be equivalent but using numerical differentiation bracket_term = ( - 1 - self.F - (self.mtr / (1 - self.mtr)) * self.eti * self.z * self.f + 1 + - 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)