Skip to content

Commit 5f27f7e

Browse files
authored
Merge pull request #36 from john-p-ryan/ETI_ODE
Solve ODE for ETI beliefs
2 parents 57f417f + bcaf193 commit 5f27f7e

File tree

2 files changed

+173
-64
lines changed

2 files changed

+173
-64
lines changed

examples/Simulate_all_policies.ipynb

Lines changed: 144 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010
},
1111
{
1212
"cell_type": "code",
13-
"execution_count": null,
13+
"execution_count": 1,
1414
"metadata": {},
1515
"outputs": [],
1616
"source": [
@@ -617,7 +617,109 @@
617617
},
618618
{
619619
"cell_type": "code",
620-
"execution_count": 16,
620+
"execution_count": 17,
621+
"metadata": {},
622+
"outputs": [],
623+
"source": [
624+
"eti_dict = {\n",
625+
" \"eti_values\": [0.18, 0.106, 0.567, 1.83, 1.9],\n",
626+
" \"knot_points\": [30000, 75000, 250000, 2000000, 10000000]\n",
627+
"}\n",
628+
"# ETI values from Gruber and Saez (2002) (Table 3) and Saez (2004) (Tables 2, 4, 5)\n",
629+
"# Compute MTR schedule under current law\n",
630+
"iot_2023 = iot_user.iot_comparison(\n",
631+
" policies=[{}],\n",
632+
" baseline_policies=[None],\n",
633+
" labels=[\"2023 Law\"],\n",
634+
" years=[2023],\n",
635+
" data=\"CPS\",\n",
636+
" eti=eti_dict\n",
637+
" )\n",
638+
"fig = px.line(\n",
639+
" x=iot_2023.iot[0].df().z,\n",
640+
" y=iot_2023.iot[0].df().mtr\n",
641+
" )\n",
642+
"fig.update_layout(\n",
643+
" template=template,\n",
644+
" xaxis_title=\"Wages and Salaries\",\n",
645+
" yaxis_title=r\"$T'(z)$\",\n",
646+
")\n",
647+
"fig.write_image(\n",
648+
" os.path.join(path, \"MTR_2023.png\"),\n",
649+
" scale=4\n",
650+
" )"
651+
]
652+
},
653+
{
654+
"cell_type": "code",
655+
"execution_count": null,
656+
"metadata": {},
657+
"outputs": [],
658+
"source": [
659+
"# Thought experiment: What beliefs about ETI do the candidates need to \n",
660+
"# justify their policies if we assume they are utilitarian? \n",
661+
"# If the elasticities are wildly counterfactual, we can reject this hypothesis\n",
662+
"\n",
663+
"eti_utilitarian = np.zeros((len(iot_all.iot), len(iot_all.iot[0].df().z)))\n",
664+
"for i in range(len(iot_all.iot)):\n",
665+
" eti_utilitarian[i, :] = iot.inverse_optimal_tax.find_eti(\n",
666+
" iot_all.iot[i], g_z = np.ones(len(iot_all.iot[i].df().z))\n",
667+
" )\n",
668+
"\n",
669+
"# Convert the 2D array to a DataFrame for plotting\n",
670+
"eti_df = pd.DataFrame(eti_utilitarian.T, columns=labels)\n",
671+
"\n",
672+
"fig_eti = px.line(\n",
673+
" eti_df,\n",
674+
" x=iot_all.iot[0].df().z,\n",
675+
" y=eti_df.columns,\n",
676+
" labels={\"x\": \"Wages and Salaries\", \"value\": r\"$\\varepsilon$\", \"variable\": \"Candidate\"},\n",
677+
")\n",
678+
"fig_eti.update_layout(\n",
679+
" template=template,\n",
680+
")\n",
681+
"\n",
682+
"fig_eti.update_traces(\n",
683+
" line=dict(dash=\"dot\", color=\"blue\"),\n",
684+
" selector=dict(name=\"Obama 2015\")\n",
685+
")\n",
686+
"fig_eti.update_traces(\n",
687+
" line=dict(dash=\"dot\", color=\"red\"),\n",
688+
" selector=dict(name=\"Romney 2012\")\n",
689+
")\n",
690+
"fig_eti.update_traces(\n",
691+
" line=dict(dash=\"dash\", color=\"blue\"),\n",
692+
" selector=dict(name=\"Clinton 2016\")\n",
693+
")\n",
694+
"fig_eti.update_traces(\n",
695+
" line=dict(dash=\"dash\", color=\"red\"),\n",
696+
" selector=dict(name=\"Trump 2016\")\n",
697+
")\n",
698+
"fig_eti.update_traces(\n",
699+
" line=dict(dash=\"dashdot\", color=\"blue\"),\n",
700+
" selector=dict(name=\"Biden 2020\")\n",
701+
")\n",
702+
"fig_eti.update_traces(\n",
703+
" line=dict(dash=\"dashdot\", color=\"red\"),\n",
704+
" selector=dict(name=\"Trump 2020\")\n",
705+
")\n",
706+
"fig_eti.update_traces(\n",
707+
" line=dict(dash=\"solid\", color=\"blue\"),\n",
708+
" selector=dict(name=\"Harris 2024\")\n",
709+
")\n",
710+
"fig_eti.update_traces(\n",
711+
" line=dict(dash=\"solid\", color=\"red\"),\n",
712+
" selector=dict(name=\"Trump 2024\")\n",
713+
")\n",
714+
"fig_eti.write_image(\n",
715+
" os.path.join(path, \"eti_utilitarian.png\"),\n",
716+
" scale=4\n",
717+
" )"
718+
]
719+
},
720+
{
721+
"cell_type": "code",
722+
"execution_count": null,
621723
"metadata": {},
622724
"outputs": [],
623725
"source": [
@@ -626,6 +728,7 @@
626728
"# one plot with epsilon(z) for each candidate\n",
627729
"# Will pick Trump and Clinton (2016) for example\n",
628730
"\n",
731+
"\n",
629732
"# First, plot just their g(z)\n",
630733
"fig = px.line(\n",
631734
" x=iot_all.iot[2].df().z[10:],\n",
@@ -649,56 +752,59 @@
649752
" os.path.join(path, \"trump_clinton_g_z_numerical.png\"),\n",
650753
" scale=4\n",
651754
" )\n",
652-
"# Now find the epsilon(z) that would give Trump's policies the same g(z) as Clinton\n",
653-
"eti_beliefs_lw, eti_beliefs_jjz = iot.inverse_optimal_tax.find_eti(iot_all.iot[2], iot_all.iot[3], g_z_type=\"g_z\")\n",
654-
"idx = np.where(np.absolute(eti_beliefs_jjz[1:]) < 10)[0]\n",
655-
"fig2 = px.line(\n",
656-
" x=iot_all.iot[2].df().z[idx],\n",
657-
" y=eti_beliefs_jjz[idx],\n",
658-
" labels={\"x\": \"Wages and Salaries\", \"y\": r\"$\\text{Implied } \\varepsilon$\"},\n",
755+
"\n",
756+
"\n",
757+
"# What ETI does Clinton need to justify Trump's SWW given her policy?\n",
758+
"eti_clinton = iot.inverse_optimal_tax.find_eti(\n",
759+
" iot_all.iot[2], \n",
760+
" g_z=iot_all.iot[3].df().g_z)\n",
761+
"\n",
762+
"fig_eti_clinton = px.line(\n",
763+
" x=iot_all.iot[2].df().z,\n",
764+
" y=eti_clinton,\n",
765+
" labels={\"x\": \"Wages and Salaries\", \"y\": r\"$\\varepsilon$\"},\n",
659766
" )\n",
660-
"fig2.update_layout(\n",
767+
"fig_eti_clinton.update_layout(\n",
661768
" template=template,\n",
662769
")\n",
663-
"fig2.write_image(\n",
664-
" os.path.join(path, \"trump_eti.png\"),\n",
770+
"fig_eti_clinton.update_traces(\n",
771+
" line=dict(dash=\"dash\", color=\"blue\"),\n",
772+
" selector=dict(name=\"Clinton 2016\")\n",
773+
")\n",
774+
"fig_eti_clinton.write_image(\n",
775+
" os.path.join(path, \"eti_clinton.png\"),\n",
665776
" scale=4\n",
666-
" )"
777+
" )\n"
667778
]
668779
},
669780
{
670781
"cell_type": "code",
671-
"execution_count": 17,
782+
"execution_count": null,
672783
"metadata": {},
673784
"outputs": [],
674785
"source": [
675-
"eti_dict = {\n",
676-
" \"eti_values\": [0.18, 0.106, 0.567, 1.83, 1.9],\n",
677-
" \"knot_points\": [30000, 75000, 250000, 2000000, 10000000]\n",
678-
"}\n",
679-
"# ETI values from Gruber and Saez (2002) (Table 3) and Saez (2004) (Tables 2, 4, 5)\n",
680-
"# Compute MTR schedule under current law\n",
681-
"iot_2023 = iot_user.iot_comparison(\n",
682-
" policies=[{}],\n",
683-
" baseline_policies=[None],\n",
684-
" labels=[\"2023 Law\"],\n",
685-
" years=[2023],\n",
686-
" data=\"CPS\",\n",
687-
" eti=eti_dict\n",
688-
" )\n",
689-
"fig = px.line(\n",
690-
" x=iot_2023.iot[0].df().z,\n",
691-
" y=iot_2023.iot[0].df().mtr\n",
786+
"\n",
787+
"# Converse: What elasticty does Trump need to justify Clinton's weights?\n",
788+
"eti_trump = iot.inverse_optimal_tax.find_eti(\n",
789+
" iot_all.iot[3], \n",
790+
" g_z=iot_all.iot[2].df().g_z)\n",
791+
"\n",
792+
"fig_eti_trump = px.line(\n",
793+
" x=iot_all.iot[2].df().z,\n",
794+
" y=eti_trump,\n",
795+
" labels={\"x\": \"Wages and Salaries\", \"y\": r\"$\\varepsilon$\"},\n",
692796
" )\n",
693-
"fig.update_layout(\n",
797+
"fig_eti_trump.update_layout(\n",
694798
" template=template,\n",
695-
" xaxis_title=\"Wages and Salaries\",\n",
696-
" yaxis_title=r\"$T'(z)$\",\n",
697799
")\n",
698-
"fig.write_image(\n",
699-
" os.path.join(path, \"MTR_2023.png\"),\n",
800+
"fig_eti_trump.update_traces(\n",
801+
" line=dict(dash=\"dash\", color=\"red\"),\n",
802+
" selector=dict(name=\"Trump 2016\")\n",
803+
")\n",
804+
"fig_eti_trump.write_image(\n",
805+
" os.path.join(path, \"eti_trump.png\"),\n",
700806
" scale=4\n",
701-
" )"
807+
" )\n"
702808
]
703809
}
704810
],
@@ -718,7 +824,7 @@
718824
"name": "python",
719825
"nbconvert_exporter": "python",
720826
"pygments_lexer": "ipython3",
721-
"version": "3.12.3"
827+
"version": "3.11.7"
722828
}
723829
},
724830
"nbformat": 4,

iot/inverse_optimal_tax.py

Lines changed: 29 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,8 @@ class IOT:
3030
parametric, if None, then non-parametric bin weights
3131
mtr_smoother (None or str): method used to smooth our mtr
3232
function, if None, then use bin average mtrs
33+
mtr_smooth_param (scalar): parameter for mtr_smoother
34+
kreg_bw (array_like): bandwidth for kernel regression
3335
3436
Returns:
3537
class instance: IOT
@@ -122,6 +124,8 @@ def compute_mtr_dist(
122124
weight_var (str): name of weight measure from data to use
123125
mtr_smoother (None or str): method used to smooth our mtr
124126
function, if None, then use bin average mtrs
127+
mtr_smooth_param (scalar): parameter for mtr_smoother
128+
kreg_bw (array_like): bandwidth for kernel regression
125129
126130
Returns:
127131
tuple:
@@ -207,6 +211,7 @@ def compute_income_dist(
207211
weight_var (str): name of weight measure from data to use
208212
dist_type (None or str): type of distribution to use if
209213
parametric, if None, then non-parametric bin weights
214+
kde_bw (array_like): bandwidth for kernel regression
210215
211216
Returns:
212217
tuple:
@@ -367,7 +372,7 @@ def sw_weights(self):
367372
+ ((self.theta_z * self.eti * self.mtr) / (1 - self.mtr))
368373
+ ((self.eti * self.z * self.mtr_prime) / (1 - self.mtr) ** 2)
369374
)
370-
integral = np.trapz(g_z * self.f, self.z)
375+
integral = np.trapz(g_z * self.f, self.z) # renormalize to integrate to 1
371376
g_z = g_z / integral
372377

373378
# use Lockwood and Weinzierl formula, which should be equivalent but using numerical differentiation
@@ -386,41 +391,39 @@ def sw_weights(self):
386391
return g_z, g_z_numerical
387392

388393

389-
def find_eti(iot1, iot2, g_z_type="g_z"):
394+
def find_eti(iot, g_z = None, eti_0 = 0.25):
390395
"""
391396
This function solves for the ETI that would result in the
392-
policy represented via MTRs in iot2 be consistent with the
393-
social welfare function inferred from the policies of iot1.
397+
policy represented via MTRs in IOT being consistent with the
398+
social welfare function supplied. It solves a first order
399+
ordinary differential equation.
394400
395401
.. math::
396-
\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})
402+
\varepsilon'(z)\left[\frac{zT'(z)}{1-T'(z)}\right] + \varepsilon(z)\left[\theta_z \frac{T'(z)}{1-T'(z)} +\frac{zT''(z)}{(1-T'(z))^2}\right]+ (1-g(z))
397403
398404
Args:
399-
iot1 (IOT): IOT class instance representing baseline policy
400-
iot2 (IOT): IOT class instance representing reform policy
401-
g_z_type (str): type of social welfare function to use
402-
Options are:
403-
* 'g_z' for the analytical formula
404-
* 'g_z_numerical' for the numerical approximation
405+
iot (IOT): instance of the I
406+
g_z (None or array_like): vector of social welfare weights
407+
eti_0 (scalar): guess for ETI at z=0
405408
406409
Returns:
407410
eti_beliefs (array-like): vector of ETI beliefs over z
408411
"""
409-
if g_z_type == "g_z":
410-
g_z = iot1.g_z
411-
else:
412-
g_z = iot1.g_z_numerical
413-
# The equation below is a simplication of the above to make the integration easier
414-
eti_beliefs_lw = ((1 - iot2.mtr) / (iot2.z * iot2.f * iot2.mtr)) * (
415-
1 - iot2.F - (g_z.sum() - np.cumsum(g_z))
416-
)
417-
# derivation from JJZ analytical solution that doesn't involved integration
418-
eti_beliefs_jjz = (g_z - 1) / (
419-
(iot2.theta_z * (iot2.mtr / (1 - iot2.mtr)))
420-
+ (iot2.z * (iot2.mtr_prime / (1 - iot2.mtr) ** 2))
421-
)
422-
423-
return eti_beliefs_lw, eti_beliefs_jjz
412+
413+
if g_z is None:
414+
g_z = iot.g_z
415+
416+
# we solve an ODE of the form f'(z) + P(z)f(z) = Q(z)
417+
P_z = 1/iot.z + iot.f_prime/iot.f + iot.mtr_prime/(iot.mtr * (1-iot.mtr))
418+
# integrating factor for ODE: mu(z) * f'(z) + mu(z) * P(z) * f(z) = mu(z) * Q(z)
419+
mu_z = np.exp(np.cumsum(P_z))
420+
Q_z = (g_z - 1) * (1 - iot.mtr) / (iot.mtr * iot.z)
421+
# integrate Q(z) * mu(z), as we integrate both sides of the ODE
422+
int_mu_Q = np.cumsum(mu_z * Q_z)
423+
424+
eti_beliefs = (eti_0 + int_mu_Q) / mu_z
425+
426+
return eti_beliefs
424427

425428

426429
def wm(value, weight):

0 commit comments

Comments
 (0)