Tax function estimation

[jdebacker]

Beginning with PR #73, which updated the default calibration of OG-USA, we have observed some odd results related to the estimated tax functions. This issue will document what we've noticed in the hopes that we can address any issues with the tax function estimation routines or with the microsimulation model used to calibrate OG-USA (or both).

Things that haven't seemed quite right:

• In PR Update calibration #73, I noted that using DEP tax functions estimated using the most recent Tax-Calculator at the time (v 3.4.1) resulted in tax function parameters that, when used in OG-USA, resulted in an inability for the model SS to solve.
• Also, noted in PR Update calibration #73, when trying to estimate the mono and mono2D functional form for the tax functions, there were failures in the estimation (e.g., no minimum found) (again, using Tax-Calculator 3.4.1)
• In OG-USA simulations since October 2023, we've used GS functional forms for the tax functions (with these, the model solve), but we've noticed significant garbage collection and reductions in computational performance when solving the model (noted in OG-USA Discussions Analysis of Dask distributed workloads #83). Times to solve the model SS have gone up from about 45 seconds to 15 minutes. Note that when using the tax functions parameters in ogusa_default_parameters.json, the warnings and performance reductions pretty much disappear.

[jdebacker]

Some plots:

DEP functions estimated on Tax-Calculator 3.4.1 (each line is a different age- blues for younger, red for older):

GS functions estimated on Tax-Calculator 3.5.1:

[jdebacker]

Some key questions:

1. Are these odd functions and artifact of the microsimulation model output or txfunc.py (there have been changes to both)?
2. Are these functions "correct" (i.e., do they fit the data best)?
3. What is preventing estimation of the mono and mono2D functions?

[jdebacker]

Re (2) above, I don't see how these could be the best fit (albeit, the scatter plot dots do not reflect sampling weights):

ETRs for 40 year olds, DEP functions (tax year 2024):

MTR on labor income for 40 year olds, DEP functions (tax year 2024):

[jdebacker]

I've started looking into the estimation of the tax functions. Some questions I have:

1. Does the numerical optimization method in our minimization of the non-linear least squares estimator matter for our estimates? In particular, as seen above, it is not uncommon to see the fitted functions and see clear room for a different parameterization to fit better than what is returned from the optimizer.
2. Are we using the good starting values in our optimization? Do "better" starting value help reduce variation across age?
3. Is it particularly difficult to estimate MTRs since they display much more variation than ETRs? And if so, is it better to infer the MTRs from the ETRs? But then, how much of the variation in the data are we missing?

[jdebacker]

Re the method of numerical optimization, I'm seeing significant differences across the numerical algorithm used to minimize the nonlinear least squares function. Here are the tax functions for each age estimated using a few different algorithms:

DEP functional form:

L-BFGS-B method

CPS data

PUF data

SLSQP method

CPS data

PUF data

Nelder-Mead method

CPS data

PUF data

GS functional form:

L-BFGS-B method

CPS data

PUF data

SLSQP

CPS data

PUF data

Nelder-Mead

CPS data

PUF data

Summary:

• Quite of bit of variation across datasets (CPS vs PUF) and age -- both suggest that the estimates are very sensitive to the underlying data because there's shouldn't be that much variation in tax rates across the two data files (but we can confirm this).
• The variation across methods to minimize the statistical objective function also suggests parameter estimates that are very sensitive to initial values and algorithms and therefore probably not precisely estimated.

[jdebacker]

ETR function estimation

The above plots are of MTRs on labor income. ETRs seem to be more consistently estimated:

DEP

CPS

PUF

GS

CPS

PUF

GS, Nelder-Mead, PUF

[jdebacker]

Following up, here are comparisons of the GS, DEP, and HSV functional forms (all estimated using the TMD file and on current law policy in 2025):

GS

DEP

HSV

[jdebacker]

Function to plot by age (maybe to be added to ogcore.parameter_plots at some point, but useful here):

# Plot tax functions by age all together
def plot_taxfunc_by_age(tax_params, S, start_year=2023, year_to_plot=2023, tax_func_type="DEP", rate_type="etr", over_labinc=True,
    other_inc_val=1000):
    tp = tax_params["tfunc_" + rate_type + "_params_S"]
    t = year_to_plot - start_year
    # Set income range to plot over (min income value hard coded to 5)
    max_inc_amt = 1000000
    inc_sup = np.exp(np.linspace(np.log(5), np.log(max_inc_amt), 100))
    # Set income value for other income
    inc_fix = other_inc_val
    if over_labinc:
        key1 = "total_labinc"
        X = inc_sup
        Y = inc_fix
    else:
        key1 = "total_capinc"
        X = inc_fix
        Y = inc_sup

    fig, ax = plt.subplots()
    cm = plt.get_cmap("coolwarm")
    ax.set_prop_cycle(color=[cm(1.0 * i / S) for i in range(S)])
    for s in range(S):
        rates = txfunc.get_tax_rates(
            tp[t][s],
            X,
            Y,
            None,
            tax_func_type,
            rate_type,
            for_estimation=False,
        )
        plt.plot(inc_sup, rates, label="Age = " + str(s + 20))
    # plt.legend(loc="center right")
    # if title:
    #     plt.title(title)
    if over_labinc:
        plt.xlabel(r"Labor income")
    else:
        plt.xlabel(r"Capital income")
    plt.ylabel(rate_type.upper())
    # if path is None:
    return fig
    # else:
    #     plt.savefig(path, dpi=300)

cc @rickecon