diff --git a/changelog.d/cg-response-functional-form.added.md b/changelog.d/cg-response-functional-form.added.md new file mode 100644 index 00000000000..af676f756ca --- /dev/null +++ b/changelog.d/cg-response-functional-form.added.md @@ -0,0 +1 @@ +Add semi-elasticity and net-of-tax-elasticity forms for the capital gains realization response, selected by setting `gov.simulation.capital_gains_responses.semi_elasticity` or `net_of_tax_elasticity`. The semi-elasticity form is defined at a 0% baseline rate; in a single-rate model with a semi-elasticity above 1 it puts the revenue-maximizing rate at 1 / semi-elasticity, so JCT's 1990 coefficient of 3.51 gives its published 28.5%. The new forms reject parameters that are negative, too large or not finite. The existing log-rate `elasticity` form is unchanged. diff --git a/docs-quarto/_quarto.yml b/docs-quarto/_quarto.yml index 6ef3da5fea4..eacfdb4eeba 100644 --- a/docs-quarto/_quarto.yml +++ b/docs-quarto/_quarto.yml @@ -21,6 +21,7 @@ website: - methodology/index.qmd - methodology/spm-poverty.qmd - methodology/moop-decomposition.qmd + - methodology/capital-gains-responses.qmd - id: programs title: "Programs" contents: diff --git a/docs-quarto/methodology/capital-gains-responses.qmd b/docs-quarto/methodology/capital-gains-responses.qmd new file mode 100644 index 00000000000..843bf710eb5 --- /dev/null +++ b/docs-quarto/methodology/capital-gains-responses.qmd @@ -0,0 +1,82 @@ +--- +title: "Capital gains realization responses" +subtitle: "Three functional forms, how each treats a 0% rate, and how to calibrate them" +--- + +When a reform changes the tax rate on capital gains, people change how much gain they realize. PolicyEngine US models this as a change in each person's long-term capital gains, driven by the change in their marginal tax rate on capital gains. This page describes how that response is computed and the three functional forms it can take. It is maintained with the code in [`capital_gains_responses.py`](https://github.com/PolicyEngine/policyengine-us/blob/main/policyengine_us/variables/gov/simulation/capital_gains_responses.py) and [`capital_gains_response_forms.py`](https://github.com/PolicyEngine/policyengine-us/blob/main/policyengine_us/variables/gov/simulation/capital_gains_response_forms.py). + +## The response + +For each person, + +```text +capital_gains_behavioral_response = long_term_capital_gains_before_response × f(τ0, τ1) +long_term_capital_gains = long_term_capital_gains_before_response + capital_gains_behavioral_response +``` + +so gains after the response are the pre-response gains times the multiplier `1 + f`. The response is zero unless the simulation is a reform with a baseline and one of the three parameters under `gov.simulation.capital_gains_responses` is nonzero. All three default to zero. + +`τ0` and `τ1` are the person's marginal tax rates on capital gains under baseline and reform policy. Both come from `marginal_tax_rate_on_capital_gains`, measured in two branches that hold the same pre-response inputs and have every behavioral response switched off (`get_behavioral_response_measurements`). The rate is a household finite difference, `1 − Δ household_net_income / $1,000`, for an extra $1,000 of gains. It is measured for the two adults in each household with the largest long-term gains and is zero for everyone else, so other members do not respond. Because `household_net_income` includes benefits and refundable credits, the rate can fall outside [0, 1] where the extra $1,000 crosses a benefit cliff or phase-in. + +## The three forms + +Set exactly one parameter. Setting more than one raises an error, because the forms are alternatives and combining them would apply the response more than once. + +| Parameter | Realizations scale as | Accepted values | Single-rate revenue `τ·R` peaks at | Rates used | +|---|---|---|---|---| +| `elasticity` (log-rate) | `τ^ε` | any finite ε; usually negative | nowhere below 100% for −1 < ε ≤ 0 | floored at 0.001 | +| `semi_elasticity` | `exp(−β·τ)` | 0 < β ≤ 50 | `τ = 1/β` when β > 1 | clipped to [0, 1] | +| `net_of_tax_elasticity` | `(1 − τ)^e` | 0 < e ≤ 5 | `τ = 1/(1 + e)` when e ≥ 0.003 | clipped to [0, 0.999] | + +So the multipliers are + +```text +log-rate: 1 + f = exp(ε · [ln τ1 − ln τ0]) +semi: 1 + f = exp(−β · (τ1 − τ0)) +net-of-tax: 1 + f = exp(e · [ln(1 − τ1) − ln(1 − τ0)]) +``` + +The two new forms reject a parameter that is negative, above its maximum or not finite, at the policy level and for any person-level override. Their parameters are positive where realizations fall as the rate rises, the opposite of the log-rate elasticity's usual sign, so a semi-elasticity entered as −3.4 is an error rather than a silent reversal. The maxima keep the dollar response finite in the model's 32-bit arithmetic for any gains below $10 quadrillion; the nonzero semi-elasticities in CRS R41364's Table 1 run from 1.0 to 4.1. A new form's response that is still not finite, which can happen only if a measured rate or a person's gains is not finite or the gains exceed that $10 quadrillion bound, also raises an error. The log-rate form accepts any finite elasticity and keeps its original arithmetic. + +In every form `f = 0` when `τ1 = τ0`, and `1 + f` is never negative, so the response never flips the sign of a person's gains; for a net loss, a smaller multiplier moves the loss toward zero. With a negative log-rate elasticity and the new forms' positive parameters, the multiplier does not rise as `τ1` rises. Multipliers compose along a path of rates up to floating-point rounding, which the model's 32-bit arithmetic makes visible: for a round trip from 0% to 100% and back at β = 3.4, the product differs from 1 by at most 32 times float32 machine epsilon (about 3.8 × 10^−6); at β = 20 the product is 0 because the outbound response factor rounds to −1. In 64-bit arithmetic the same β = 20 round trip differs from 1 by at most 2 × 10^−7. In 64-bit arithmetic the composition holds to one part in 10^9 wherever each multiplier is at least 10^−6. The test suite checks these properties, and each form's single-rate revenue peak, for every input [Hypothesis](https://hypothesis.readthedocs.io/) draws (`policyengine_us/tests/core/test_capital_gains_response_forms.py`). + +These peaks are for a single flat rate on all gains. A progressive schedule, the net investment income tax, the alternative minimum tax, state taxes and benefit interactions move the model's revenue-maximizing statutory rate away from them. + +### Log-rate form + +This was the only form before the other two were added, and it remains the one used whenever `elasticity` is set; its results are unchanged. In a single-rate model revenue is proportional to `τ^(1+ε)`, so with any elasticity between −1 and 0, revenue keeps rising up to a 100% rate. The Joint Committee on Taxation (JCT) staff made the same point about this equation in 1990: it "would predict no tax rate at which revenue is maximized" ([JCS-12-90](https://www.jct.gov/getattachment/9cbc5fb2-3452-4e67-8edf-7b41fa12bdeb/jcs-12-90-3157.pdf#page=66), Appendix B, equation [5]). + +The log of a 0% rate is undefined, so both rates are floored at 0.001. A person whose baseline rate is 0%, such as someone whose gains fall in the 0% bracket in a state without an income tax, therefore gets a log change of `ln(τ1 / 0.001)`. A move from 0% to 15% is a log change of about 5.0, and at ε = −0.7 it removes about 97% of their realizations. In a household test, a person in the 0% bracket facing a new 10% rate loses 96% of their gains under ε = −0.7. + +### Semi-elasticity form + +The response depends only on the change in the rate, so it is defined at a 0% baseline without a floor. The same person loses 29% of their gains at β = 3.4. β is roughly the percent fall in realizations per percentage-point rise in the rate. + +Single-rate revenue `τ · A · exp(−βτ)` has derivative `A · exp(−βτ) · (1 − βτ)`, so it peaks at `τ = 1/β` when that is below 1. Without the clip described below, a fixed state rate `s` that adds to the federal rate `t` multiplies realizations by the constant `exp(−βs)`, so the analytical peak of single-rate federal revenue `t · R` stays at a federal rate (including the net investment income tax) of `1/β`. With the clip, realizations stop falling once `t + s` reaches 100%, and federal revenue then rises with `t`. The analytical peak is the maximum only when `β(1 − s) − ln β > 1`, which for β from 3.1 to 3.51 means a state rate below about 31% to 36%; otherwise revenue at a 100% federal rate is higher. With no state rate the condition holds for every β above 1. + +Rates are clipped to [0, 1] before the change is taken. Without the clip, a household whose measured baseline rate is 10 because of a benefit cliff, and whose reform removes the cliff, would have its realizations multiplied by `exp(9.75β)`. With it, the multiplier `1 + f` stays between `exp(−β)` and `exp(β)`. + +### Net-of-tax form + +Agersnap and Zidar (2021) estimate realization responses to the log net-of-tax rate. The [Yale Budget Lab's](https://budgetlab.yale.edu/research/revenue-potential-taxes-top) model also specifies its realization response on the net-of-tax rate, within a dynamic model of realizing or deferring gains that this constant-elasticity form does not reproduce. The log of the net-of-tax share is undefined at a 100% rate, the mirror of the log-rate form's problem at 0%, so rates are clipped to at most 0.999. Realizations therefore stop falling above a 99.9% rate, and single-rate revenue rises again there; for elasticities below about 0.0027, whose analytical peak is above 99.7%, revenue at 100% is higher than at that peak. A baseline rate at the cap gives a large response to a reform rate well below it. + +## Calibration + +| Source | Form | Value | Single-rate revenue-maximizing rate | +|---|---|---|---| +| JCT staff, 1990 ([JCS-12-90](https://www.jct.gov/getattachment/9cbc5fb2-3452-4e67-8edf-7b41fa12bdeb/jcs-12-90-3157.pdf#page=62), printed pp. 54 and 59) | `ln(realizations) = a + b·MTR`, b = −0.0351 per percentage point | β = 3.51 | 28.5% | +| JCT, as reported by the Congressional Research Service in 2021 ([R41364](https://www.congress.gov/crs_external_products/R/PDF/R41364/R41364.13.pdf#page=10), Table 1) and 2025 ([R48562](https://www.congress.gov/crs_external_products/R/PDF/R48562/R48562.1.pdf#page=11)) | Semi-log | β = 3.1 | 32.3% | +| Treasury, as reported in CRS R41364 | Semi-log | β = 3.25 | 30.8% | +| CBO time-series forecasting model (2010), as reported in CRS R41364 | Semi-log | β = 1.76 | 56.8% | +| Dowd, McClelland and Muthitacharoen ([2015](https://doi.org/10.17310/ntj.2015.3.02)), as reported in CRS R41364 | Semi-log | β = 4.1 | 24.4% | +| Agersnap and Zidar ([2021](https://doi.org/10.1257/aeri.20200535)) | Net-of-tax | | Federal: 33% to 47% across their 10-year specifications; 38% to 47% over years 6 to 10 | + +A target revenue-maximizing rate `τ*` maps to `β = 1/τ*`: 28.5% is 3.51, 29% is 3.45 and 30% is 3.33. JCT's 1990 coefficient gives 28.5%. CRS reported in 2021 that JCT used 3.1, which gives 32.3%, and in 2025 that JCT "has used" 3.1. JCT's staff also cautioned that a revenue-maximizing rate is very sensitive to the functional form chosen and that they did not find it a useful summary of taxpayer responsiveness (JCS-12-90, printed p. 5 and pp. 41–42). + +JCT's coefficient applies to the marginal rate from its individual income tax model. PolicyEngine applies `β` to the change in a household's marginal rate on gains, which also includes state income tax, the net investment income tax and benefit interactions. Because only the change enters, a federal rate change that passes through one-for-one produces the response JCT's equation would. + +## Known limitation: how the rate is measured + +`marginal_tax_rate_on_capital_gains` raises the person's `capital_gains` input by $1,000. Gross income counts `capital_gains`, but the preferential-rate computation reads `long_term_capital_gains`, which that override does not change. So the extra $1,000 reaches adjusted gross income but not net capital gain, and the measured rate can reflect ordinary rates, deductions and how gains stack on other income rather than the long-term rate. In household tests for 2026, the measured rate for a Texas filer with $100 million of long-term gains is 0.296 rather than about 0.24 (20% plus the 3.8% net investment income tax) from raising long-term gains directly, with the alternative minimum tax making up the difference; for a Texas filer with $60,000 of wages and $200,000 of gains it is 0.308 rather than 0.188. Raising the top capital gains rate from 20% to 28% leaves the first filer's measured rate unchanged, so no form produces a response to that reform. All three forms take their rates from this measurement. [PolicyEngine US #10058](https://github.com/PolicyEngine/policyengine-us/pull/10058) raises long-term gains instead. + +Separately, the measured rate stops rising once a top capital gains rate exceeds the top ordinary rate. 26 U.S.C. 1(h)(1) caps the tax on gains at the tax at ordinary rates (`limit_to_tax_at_main_rates`), so raising only the capital gains rates cannot tax large gains above about 37% plus the net investment income tax. diff --git a/docs-quarto/methodology/index.qmd b/docs-quarto/methodology/index.qmd index 16aae236a38..c7249d5ce60 100644 --- a/docs-quarto/methodology/index.qmd +++ b/docs-quarto/methodology/index.qmd @@ -6,5 +6,6 @@ How PolicyEngine US is built and how to read its outputs. - [How PolicyEngine US calculates SPM poverty](spm-poverty.qmd) — population, thresholds from spm-calculator, resources, and what the model does not yet do. - [Medical out-of-pocket expense decomposition](moop-decomposition.qmd) — the architecture for layering rules-based premium computations on top of CPS-imputed medical costs. +- [Capital gains realization responses](capital-gains-responses.qmd) — the log-rate, semi-elasticity and net-of-tax forms, how each treats a 0% rate, and calibration to a revenue-maximizing rate. -Planned sections (not yet written): microdata construction (enhanced CPS), calibration to administrative totals, behavioral responses, uprating, and validation against NBER TAXSIM and IRS SOI. +Planned sections (not yet written): microdata construction (enhanced CPS), calibration to administrative totals, labor supply responses, uprating, and validation against NBER TAXSIM and IRS SOI. diff --git a/policyengine_us/parameters/gov/simulation/capital_gains_responses/elasticity.yaml b/policyengine_us/parameters/gov/simulation/capital_gains_responses/elasticity.yaml index 0afdba34c1b..81445c84ccf 100644 --- a/policyengine_us/parameters/gov/simulation/capital_gains_responses/elasticity.yaml +++ b/policyengine_us/parameters/gov/simulation/capital_gains_responses/elasticity.yaml @@ -1,6 +1,22 @@ -description: Elasticity of capital gains with respect to the capital gains marginal tax rate. +description: Elasticity of long-term capital gains realizations with respect to the capital gains marginal tax rate, applied to the log change in the rate. +documentation: >- + Selects the log-rate form, under which realizations scale as the rate + raised to this power. Negative values reduce realizations when the rate + rises. In a single-rate model revenue is then proportional to the rate + raised to (1 + elasticity), so for any value between -1 and 0 revenue rises + all the way to a 100% rate and the form has no revenue-maximizing rate, as + the Joint Committee on Taxation staff noted of this form in 1990; use + semi_elasticity for a revenue-maximizing rate. Both the baseline and reform + rates are floored at 0.001 before taking logs, so a person whose baseline + rate is 0% gets a log change of ln(reform rate / 0.001). A move from 0% to + 15% is a log change of about 5.0, which at an elasticity of -0.7 removes + about 97% of that person's realizations. At most one of elasticity, + semi_elasticity and net_of_tax_elasticity can be nonzero. values: 2000-01-01: 0 metadata: unit: /1 label: capital gains elasticity + reference: + - title: Joint Committee on Taxation, Explanation of Methodology Used to Estimate Proposals Affecting the Taxation of Income from Capital Gains (JCS-12-90) + href: https://www.jct.gov/getattachment/9cbc5fb2-3452-4e67-8edf-7b41fa12bdeb/jcs-12-90-3157.pdf#page=66 diff --git a/policyengine_us/parameters/gov/simulation/capital_gains_responses/net_of_tax_elasticity.yaml b/policyengine_us/parameters/gov/simulation/capital_gains_responses/net_of_tax_elasticity.yaml new file mode 100644 index 00000000000..0b7d75c2ff6 --- /dev/null +++ b/policyengine_us/parameters/gov/simulation/capital_gains_responses/net_of_tax_elasticity.yaml @@ -0,0 +1,21 @@ +description: Elasticity of long-term capital gains realizations with respect to one minus the capital gains marginal tax rate. +documentation: >- + Selects the net-of-tax form, under which realizations scale as (1 - rate) + raised to this power. Values must be positive, the sign under which + realizations fall as the rate rises, and at most 5. In a single-rate model, + revenue peaks at a rate of 1 / (1 + net_of_tax_elasticity) for values of at + least 0.003. Agersnap and Zidar (2021) estimate realization responses to the + log net-of-tax rate. Both rates are clipped to between 0 and 0.999 before + taking logs, so realizations stop falling above a 99.9% rate, and a + baseline rate at or above 99.9% (possible where the extra gain crosses a + benefit cliff) gives a large log change for a reform rate well below 99.9%. + At most one of elasticity, semi_elasticity and net_of_tax_elasticity can be + nonzero. +values: + 2000-01-01: 0 +metadata: + unit: /1 + label: capital gains net-of-tax elasticity + reference: + - title: Agersnap and Zidar, The Tax Elasticity of Capital Gains and Revenue-Maximizing Rates, American Economic Review Insights 3(4), 2021 + href: https://doi.org/10.1257/aeri.20200535 diff --git a/policyengine_us/parameters/gov/simulation/capital_gains_responses/semi_elasticity.yaml b/policyengine_us/parameters/gov/simulation/capital_gains_responses/semi_elasticity.yaml new file mode 100644 index 00000000000..8482aec8086 --- /dev/null +++ b/policyengine_us/parameters/gov/simulation/capital_gains_responses/semi_elasticity.yaml @@ -0,0 +1,37 @@ +description: Semi-elasticity of long-term capital gains realizations with respect to the capital gains marginal tax rate, about the percent fall in realizations per percentage-point rise in the rate. +documentation: >- + Selects the semi-elasticity form, under which realizations scale by + exp(-semi_elasticity * (reform rate - baseline rate)). Values must be + positive, the sign under which realizations fall as the rate rises (unlike + the negative log-rate elasticity), and at most 50. In a single-rate model, + revenue peaks at a rate of 1 / semi_elasticity when that is below 100%, so + 3.33 to 3.51 puts the peak at 30% to 28.5%; a progressive schedule and other + taxes move the model's peak statutory rate away from that. The Joint + Committee on Taxation staff's 1990 equation, ln(realizations) = a - 0.0351 * + rate in percentage points, is this form with 3.51 and a 28.5% peak. The + Congressional Research Service reported in 2021 that JCT used 3.1 (a 32.3% + peak) and in 2025 that JCT has used 3.1. Because the response depends only + on the change in the rate, the form is defined when the baseline rate is + 0%. Without clipping, a fixed state rate that adds to the federal rate + scales realizations by a constant, so the analytical peak of single-rate + federal revenue stays at a federal rate of 1 / semi_elasticity; once the + combined rate reaches 100% the clip stops realizations falling, and that + peak is the maximum only for a state rate below about 31% to 36% at + semi-elasticities of 3.1 to 3.51. Both rates are clipped to between 0 and + 1 first, because the household marginal rate can fall outside that range + where the extra gain crosses a benefit cliff; this keeps the realization + multiplier between exp(-semi_elasticity) and exp(semi_elasticity). At most + one of elasticity, semi_elasticity and net_of_tax_elasticity can be + nonzero. +values: + 2000-01-01: 0 +metadata: + unit: /1 + label: capital gains semi-elasticity + reference: + - title: "Joint Committee on Taxation, Explanation of Methodology Used to Estimate Proposals Affecting the Taxation of Income from Capital Gains (JCS-12-90) #page 62,67" + href: https://www.jct.gov/getattachment/9cbc5fb2-3452-4e67-8edf-7b41fa12bdeb/jcs-12-90-3157.pdf#page=62 + - title: "Congressional Research Service, Capital Gains Tax Options: Behavioral Responses and Revenues (R41364, 2021) #page 10,16" + href: https://www.congress.gov/crs_external_products/R/PDF/R41364/R41364.13.pdf#page=10 + - title: "Congressional Research Service, Boundaries on the Long-Run Realization Response to Changes in Capital Gains Taxes (R48562, 2025)" + href: https://www.congress.gov/crs_external_products/R/PDF/R48562/R48562.1.pdf#page=11 diff --git a/policyengine_us/tests/core/test_behavioral_response_measurements.py b/policyengine_us/tests/core/test_behavioral_response_measurements.py index dd488b51a7d..92bae6f447e 100644 --- a/policyengine_us/tests/core/test_behavioral_response_measurements.py +++ b/policyengine_us/tests/core/test_behavioral_response_measurements.py @@ -133,6 +133,8 @@ def make_parameters( income_change_bound=0.3, wage_change_bound=0.4, capital_gains_elasticity=-0.62, + capital_gains_semi_elasticity=0.0, + capital_gains_net_of_tax_elasticity=0.0, lsr_income_elasticity=-0.05, lsr_substitution_elasticity=0.25, ): @@ -150,7 +152,9 @@ def make_parameters( ), ), capital_gains_responses=SimpleNamespace( - elasticity=capital_gains_elasticity + elasticity=capital_gains_elasticity, + semi_elasticity=capital_gains_semi_elasticity, + net_of_tax_elasticity=capital_gains_net_of_tax_elasticity, ), ) ) @@ -597,3 +601,62 @@ def test_capital_gains_channel_is_live_end_to_end(): # A 10pp statutory rate increase must show through the measured MTRs. assert reform_mtr[0] - baseline_mtr[0] == pytest.approx(0.10, abs=0.02) assert response[0] < -1_000 + + # With only the elasticity set, the response is the formula main had + # before the semi-elasticity and net-of-tax forms were added (75cdd801), + # bit pattern for bit pattern: both rates floored at 0.001, the log + # change times the person's elasticity, exp of that minus 1, times gains. + elasticity = simulation.calculate("capital_gains_elasticity", year) + log_change = np.log(np.maximum(measurements["reform_capital_gains_mtr"], 0.001)) + log_change -= np.log(np.maximum(measurements["baseline_capital_gains_mtr"], 0.001)) + main_response = simulation.calculate( + "long_term_capital_gains_before_response", year + ) * (np.exp(elasticity * log_change) - 1) + assert response.dtype == main_response.dtype == np.float32 + assert np.array_equal(response.view(np.uint32), main_response.view(np.uint32)) + + +def test_semi_elasticity_capital_gains_response_end_to_end_from_a_zero_rate(): + """The semi-elasticity form through the model, from a 0% baseline rate. + + A Texas filer whose $30,000 of gains sit in the 0% bracket faces a new 10% + rate there. The semi-elasticity response depends only on the change in the + measured rate, so it is about 29% at 3.4, where the log-rate form's 0.001 + floor would remove about 96% (test_capital_gains_response_forms.py covers + the forms' arithmetic without a simulation). + """ + from policyengine_us import Simulation + + year = "2026" + gains = 30_000 + situation = { + "people": {"p": {"age": {year: 50}, "long_term_capital_gains": {year: gains}}}, + "tax_units": {"tu": {"members": ["p"]}}, + "families": {"f": {"members": ["p"]}}, + "spm_units": {"s": {"members": ["p"]}}, + "households": {"h": {"members": ["p"], "state_name": {year: "TX"}}}, + } + reform = { + "gov.irs.capital_gains.rates.1": {"2026-01-01.2026-12-31": 0.10}, + "gov.simulation.capital_gains_responses.semi_elasticity": { + "2026-01-01.2026-12-31": 3.4 + }, + } + + simulation = Simulation( + situation=situation, reform=reform, spm={"geography_kind": "national"} + ) + response = simulation.calculate("capital_gains_behavioral_response", year)[0] + measurements = simulation._behavioral_response_measurements[year] + baseline_rate = float(measurements["baseline_capital_gains_mtr"][0]) + reform_rate = float(measurements["reform_capital_gains_mtr"][0]) + + assert baseline_rate == 0 + assert reform_rate == pytest.approx(0.10, abs=1e-6) + assert response == pytest.approx( + gains * np.expm1(-3.4 * (reform_rate - baseline_rate)), rel=1e-5 + ) + assert response / gains == pytest.approx(-0.288, abs=0.001) + assert simulation.calculate("long_term_capital_gains", year)[0] == pytest.approx( + gains + response, rel=1e-6 + ) diff --git a/policyengine_us/tests/core/test_capital_gains_response_forms.py b/policyengine_us/tests/core/test_capital_gains_response_forms.py new file mode 100644 index 00000000000..079c38a08a5 --- /dev/null +++ b/policyengine_us/tests/core/test_capital_gains_response_forms.py @@ -0,0 +1,758 @@ +"""Functional forms for the capital gains realization response. + +``f(tau0, tau1)`` is the response factor: realizations after the response are +the pre-response gains times the multiplier ``1 + f``. The properties below +hold for every input Hypothesis draws, within the stated domains: + +1. No rate change, no response: ``f(tau, tau) == 0`` for every form, parameter + and rate, including rates outside [0, 1]. +2. Monotonicity: with a log-rate elasticity <= 0 and the new forms' parameters + >= 0, the multiplier is non-increasing in the reform rate, and strictly + decreasing inside each form's unclipped range when the parameter is nonzero + and the multiplier resolves the step. +3. Sign preservation: ``1 + f >= 0``, so the response never flips the sign of + a person's gains, and in the model's float32 arithmetic the new forms' + dollar responses stay finite and sign-preserving across their supported + parameter domains for any gains below $10 quadrillion. +4. Path independence: ``(1 + f(t0, t1)) * (1 + f(t1, t2)) == 1 + f(t0, t2)`` + to 1e-9 in float64 where each multiplier is at least 1e-6. +5. Single-rate revenue: ``tau * (1 + f(tau0, tau))`` peaks at ``1 / beta`` under + the semi-elasticity form for beta > 1 and at ``1 / (1 + e)`` under the + net-of-tax form for e >= 0.003, whatever the baseline rate, up to the + rounding of a multiplier near 0. Under the log-rate form with an + elasticity between -1 and 0 it rises all the way to a 100% rate. +6. Bounds: the semi-elasticity factor lies in + ``[exp(-beta) - 1, exp(beta) - 1]`` for any measured rates. +7. Differential: the log-rate form, which stays the default, equals the + formula on main before the other forms were added, bit pattern for bit + pattern. + +Two documented behaviors are tested as intended rather than as violations: +below a net-of-tax elasticity of about 0.0027 the 99.9% cap makes revenue at +100% exceed the analytical peak, and a float32 round trip at a large +semi-elasticity loses the multiplier entirely. + +Household simulations of these forms are in test_behavioral_response_measurements.py, +which CI runs in the heavy-b group. +""" + +from types import SimpleNamespace + +import numpy as np +import pytest +from hypothesis import assume, example, given, settings +from hypothesis import strategies as st +from hypothesis.extra import numpy as hnp + +import policyengine_us.variables.gov.simulation.capital_gains_responses as capital_gains_module +from policyengine_us.variables.gov.simulation.capital_gains_response_forms import ( + CAPITAL_GAINS_RESPONSE_FORM_VARIABLES, + CAPITAL_GAINS_RESPONSE_FORMS, + CAPITAL_GAINS_RESPONSE_MAXIMUM, + LOG_RATE_ELASTICITY, + NET_OF_TAX_ELASTICITY, + NET_OF_TAX_SHARE_FLOOR, + SEMI_ELASTICITY, + capital_gains_response_factor, + log_rate_elasticity_response_factor, + net_of_tax_elasticity_response_factor, + selected_capital_gains_response_form, + semi_elasticity_response_factor, +) +from policyengine_us.variables.gov.simulation.capital_gains_responses import ( + capital_gains_behavioral_response, + capital_gains_net_of_tax_elasticity, + capital_gains_semi_elasticity, +) + +# Revenue-maximizing rates of 28.5% to 30% correspond to these semi-elasticities. +CALIBRATED_SEMI_ELASTICITIES = (1 / 0.285, 1 / 0.29, 1 / 0.30) + +# Hypothesis input ranges. +RATES = st.floats(min_value=0, max_value=1, allow_nan=False) +WIDE_RATES = st.floats(min_value=-5, max_value=20, allow_nan=False) +LOG_RATE_ELASTICITIES = st.floats(min_value=-3, max_value=0, allow_nan=False) +# The new forms' supported domains. +SEMI_ELASTICITIES = st.floats( + min_value=0, max_value=CAPITAL_GAINS_RESPONSE_MAXIMUM[SEMI_ELASTICITY] +) +NET_OF_TAX_ELASTICITIES = st.floats( + min_value=0, max_value=CAPITAL_GAINS_RESPONSE_MAXIMUM[NET_OF_TAX_ELASTICITY] +) +PARAMETERS_BY_FORM = { + LOG_RATE_ELASTICITY: LOG_RATE_ELASTICITIES, + SEMI_ELASTICITY: SEMI_ELASTICITIES, + NET_OF_TAX_ELASTICITY: NET_OF_TAX_ELASTICITIES, +} +FORM_AND_PARAMETER = st.sampled_from(CAPITAL_GAINS_RESPONSE_FORMS).flatmap( + lambda form: st.tuples(st.just(form), PARAMETERS_BY_FORM[form]) +) +# Rates a tax can take over which each form's factor is strictly monotonic. +UNCLIPPED_RATE_RANGE = { + LOG_RATE_ELASTICITY: (0.001, 1.0), + SEMI_ELASTICITY: (0.0, 1.0), + NET_OF_TAX_ELASTICITY: (0.0, 1 - NET_OF_TAX_SHARE_FLOOR), +} +# Allowance for floating-point rounding in revenue comparisons. The multiplier +# 1 + f is computed as 1 + expm1(x), accurate to about 2.2e-16 in absolute +# terms, so where it is near 0 (realizations a trillionth of their baseline or +# less) neighboring rates can round to the same or reversed values. Revenue, +# a rate of at most 1 times the multiplier, is compared to within a relative +# 1e-12 plus an absolute 1e-15. +REVENUE_RTOL = 1e-12 +REVENUE_ATOL = 1e-15 + + +def at_most(smaller, larger): + return smaller <= larger * (1 + REVENUE_RTOL) + REVENUE_ATOL + + +def main_log_rate_response_factor(baseline_mtr, reform_mtr, elasticity): + """The response factor on main before this change (75cdd801), kept frozen. + + ``capital_gains_behavioral_response`` computed + ``np.exp(elasticity * tax_rate_change) - 1`` with the tax rate change from + ``calculate_relative_capital_gains_mtr_change``, which floored both rates + at 0.001 and took the difference of their logs. + """ + min_rate = 0.001 + baseline = np.maximum(baseline_mtr, min_rate) + reform = np.maximum(reform_mtr, min_rate) + tax_rate_change = np.log(reform) - np.log(baseline) + return np.exp(elasticity * tax_rate_change) - 1 + + +def bits(array): + """An array's bit patterns, so signed zeros and NaN payloads compare too.""" + array = np.asarray(array) + return array.view(np.uint32 if array.dtype == np.float32 else np.uint64) + + +def factor(form, baseline_rate, reform_rate, parameter): + return float( + capital_gains_response_factor( + form, np.float64(baseline_rate), np.float64(reform_rate), parameter + ) + ) + + +def revenue(form, baseline_rate, rate, parameter): + return rate * (1 + factor(form, baseline_rate, rate, parameter)) + + +def capital_gains_responses(**values): + defaults = {form: 0.0 for form in CAPITAL_GAINS_RESPONSE_FORMS} + return SimpleNamespace(**{**defaults, **values}) + + +def make_parameters(**values): + return lambda period: SimpleNamespace( + gov=SimpleNamespace( + simulation=SimpleNamespace( + capital_gains_responses=capital_gains_responses(**values) + ) + ) + ) + + +class FakePerson: + def __init__(self, values): + self.simulation = SimpleNamespace(baseline=object()) + self.values = values + + def __call__(self, variable, period, options=None): + return self.values[variable] + + +# Known values for each form. + + +def test_log_rate_form_matches_hand_calculation(): + # Doubling the rate at an elasticity of -1 halves realizations. + assert factor(LOG_RATE_ELASTICITY, 0.2, 0.4, -1.0) == pytest.approx(-0.5) + # (0.25 / 0.20) ** -0.7 - 1 + assert factor(LOG_RATE_ELASTICITY, 0.2, 0.25, -0.7) == pytest.approx(1.25**-0.7 - 1) + + +def test_log_rate_form_floors_a_zero_baseline_rate_at_one_tenth_of_a_percent(): + # A move from 0% to 15% is a log change of ln(0.15 / 0.001) = ln(150), + # about 5.0, which at an elasticity of -0.7 removes about 97% of gains. + assert np.log(0.15 / 0.001) == pytest.approx(5.01, abs=0.01) + assert factor(LOG_RATE_ELASTICITY, 0.0, 0.15, -0.7) == pytest.approx(150**-0.7 - 1) + assert factor(LOG_RATE_ELASTICITY, 0.0, 0.15, -0.7) == pytest.approx( + -0.970, abs=0.001 + ) + # The floor also applies to negative measured rates, and there is no cap. + assert factor(LOG_RATE_ELASTICITY, -0.3, 0.15, -0.7) == factor( + LOG_RATE_ELASTICITY, 0.0, 0.15, -0.7 + ) + assert factor(LOG_RATE_ELASTICITY, 0.2, 3.0, -0.7) == pytest.approx(15**-0.7 - 1) + + +def test_semi_elasticity_form_matches_hand_calculation(): + assert factor(SEMI_ELASTICITY, 0.2, 0.3, 3.4) == pytest.approx(np.exp(-0.34) - 1) + assert factor(SEMI_ELASTICITY, 0.3, 0.2, 3.4) == pytest.approx(np.exp(0.34) - 1) + + +def test_semi_elasticity_form_is_well_defined_at_a_zero_baseline_rate(): + # Only the change in the rate matters, so a 0% baseline needs no floor: + # a move from 0% to 15% at beta = 3.4 removes 1 - exp(-0.51), about 40%, + # where the log-rate form at -0.7 removes about 97%. + assert factor(SEMI_ELASTICITY, 0.0, 0.15, 3.4) == pytest.approx(np.exp(-0.51) - 1) + assert factor(SEMI_ELASTICITY, 0.0, 0.15, 3.4) == pytest.approx(-0.40, abs=0.001) + assert factor(SEMI_ELASTICITY, 0.0, 0.15, 3.4) == factor( + SEMI_ELASTICITY, 0.10, 0.25, 3.4 + ) + + +def test_semi_elasticity_form_clips_rates_to_zero_and_one(): + # A benefit cliff can put the household finite-difference rate above 1. + assert factor(SEMI_ELASTICITY, 3.0, 0.25, 3.4) == pytest.approx( + np.exp(3.4 * 0.75) - 1 + ) + assert factor(SEMI_ELASTICITY, -0.4, 0.25, 3.4) == pytest.approx( + np.exp(-3.4 * 0.25) - 1 + ) + assert factor(SEMI_ELASTICITY, 3.0, 5.0, 3.4) == 0 + + +def test_net_of_tax_form_matches_hand_calculation(): + # ((1 - 0.4) / (1 - 0.2)) ** 0.5 - 1 + assert factor(NET_OF_TAX_ELASTICITY, 0.2, 0.4, 0.5) == pytest.approx(0.75**0.5 - 1) + assert factor(NET_OF_TAX_ELASTICITY, 0.0, 0.1, 0.7) == pytest.approx(0.9**0.7 - 1) + + +def test_net_of_tax_form_floors_the_net_of_tax_share(): + # Rates at or above 99.9% are treated as 99.9%. + assert factor(NET_OF_TAX_ELASTICITY, 1.0, 0.5, 0.5) == pytest.approx( + (0.5 / NET_OF_TAX_SHARE_FLOOR) ** 0.5 - 1 + ) + assert factor(NET_OF_TAX_ELASTICITY, 4.0, 0.5, 0.5) == factor( + NET_OF_TAX_ELASTICITY, 0.999, 0.5, 0.5 + ) + + +def test_calibrated_semi_elasticities_put_the_revenue_peak_at_28_5_to_30_percent(): + grid = np.linspace(0, 1, 100_001) + for semi_elasticity, peak in zip(CALIBRATED_SEMI_ELASTICITIES, (0.285, 0.29, 0.30)): + revenues = grid * ( + 1 + semi_elasticity_response_factor(0.238, grid, semi_elasticity) + ) + assert grid[np.argmax(revenues)] == pytest.approx(peak, abs=1e-5) + assert CALIBRATED_SEMI_ELASTICITIES[0] == pytest.approx(3.509, abs=0.001) + assert CALIBRATED_SEMI_ELASTICITIES[-1] == pytest.approx(3.333, abs=0.001) + + +def test_vectorized_factors_match_scalar_factors(): + baseline = np.array([0.0, 0.15, 0.238, 0.5, 1.2], dtype=np.float32) + reform = np.array([0.1, 0.25, 0.338, 0.4, 0.3], dtype=np.float32) + for form, parameter in ( + (LOG_RATE_ELASTICITY, -0.7), + (SEMI_ELASTICITY, 3.4), + (NET_OF_TAX_ELASTICITY, 0.7), + ): + vectorized = capital_gains_response_factor(form, baseline, reform, parameter) + assert vectorized.shape == baseline.shape + for i in range(len(baseline)): + assert vectorized[i] == pytest.approx( + factor(form, baseline[i], reform[i], parameter), rel=1e-5 + ) + + +# Selecting the form. + + +def test_no_form_is_selected_when_every_parameter_is_zero(): + assert selected_capital_gains_response_form(capital_gains_responses()) is None + + +@pytest.mark.parametrize("form", CAPITAL_GAINS_RESPONSE_FORMS) +def test_the_nonzero_parameter_selects_its_form(form): + responses = capital_gains_responses(**{form: 0.5}) + assert selected_capital_gains_response_form(responses) == form + + +@pytest.mark.parametrize( + "values", + [ + {"elasticity": -0.7, "semi_elasticity": 3.4}, + {"elasticity": -0.7, "net_of_tax_elasticity": 0.7}, + {"semi_elasticity": 3.4, "net_of_tax_elasticity": 0.7}, + {"elasticity": -0.7, "semi_elasticity": 3.4, "net_of_tax_elasticity": 0.7}, + ], +) +def test_more_than_one_nonzero_parameter_is_an_error(values): + with pytest.raises(ValueError, match="at most one"): + selected_capital_gains_response_form(capital_gains_responses(**values)) + + +@pytest.mark.parametrize( + "form, value", + [ + # A semi-elasticity entered with the log-rate elasticity's sign. + (SEMI_ELASTICITY, -3.4), + (SEMI_ELASTICITY, 50.5), + (NET_OF_TAX_ELASTICITY, -0.7), + (NET_OF_TAX_ELASTICITY, 5.5), + ], +) +def test_new_forms_reject_parameters_outside_their_domain(form, value): + with pytest.raises(ValueError, match="must be between 0 and"): + selected_capital_gains_response_form(capital_gains_responses(**{form: value})) + + +@pytest.mark.parametrize("form", CAPITAL_GAINS_RESPONSE_FORMS) +@pytest.mark.parametrize("value", [np.nan, np.inf, -np.inf]) +def test_every_form_rejects_a_parameter_that_is_not_finite(form, value): + with pytest.raises(ValueError, match="must be finite"): + selected_capital_gains_response_form(capital_gains_responses(**{form: value})) + + +@pytest.mark.parametrize("value", [-5.0, -0.7, 0.5]) +def test_log_rate_form_accepts_any_finite_elasticity(value): + responses = capital_gains_responses(elasticity=value) + assert selected_capital_gains_response_form(responses) == LOG_RATE_ELASTICITY + + +@pytest.mark.parametrize("form", [SEMI_ELASTICITY, NET_OF_TAX_ELASTICITY]) +def test_new_forms_accept_their_maximum(form): + maximum = CAPITAL_GAINS_RESPONSE_MAXIMUM[form] + responses = capital_gains_responses(**{form: maximum}) + assert selected_capital_gains_response_form(responses) == form + + +def test_defaults_select_no_response_in_every_year(): + from policyengine_us.system import system + + parameters = system.parameters + for year in range(2000, 2101, 5): + responses = parameters(f"{year}-01-01").gov.simulation.capital_gains_responses + assert responses.elasticity == 0 + assert responses.semi_elasticity == 0 + assert responses.net_of_tax_elasticity == 0 + assert selected_capital_gains_response_form(responses) is None + + +# The behavioral response variable. + + +MEASUREMENTS = { + "baseline_capital_gains_mtr": np.array([0.0, 0.2, 0.238], dtype=np.float32), + "reform_capital_gains_mtr": np.array([0.1, 0.3, 0.198], dtype=np.float32), +} +GAINS = np.array([100.0, 200.0, -50.0], dtype=np.float32) + + +def response_with(monkeypatch, measurements=MEASUREMENTS, **values): + monkeypatch.setattr( + capital_gains_module, + "get_behavioral_response_measurements", + lambda person, period: measurements, + ) + person = FakePerson( + { + "long_term_capital_gains_before_response": GAINS, + **{ + CAPITAL_GAINS_RESPONSE_FORM_VARIABLES[form]: np.full( + 3, value, dtype=np.float32 + ) + for form, value in values.items() + }, + } + ) + return capital_gains_behavioral_response.formula( + person, 2026, make_parameters(**values) + ) + + +def test_response_uses_the_semi_elasticity_form(monkeypatch): + response = response_with(monkeypatch, semi_elasticity=3.4) + change = ( + MEASUREMENTS["reform_capital_gains_mtr"] + - MEASUREMENTS["baseline_capital_gains_mtr"] + ) + assert np.allclose(response, GAINS * np.expm1(-3.4 * change), rtol=1e-6) + # The first two rates rise, so their gains fall. The third rate falls, so + # that person's realizations grow, which makes their net loss larger. + assert response[0] < 0 and response[1] < 0 and response[2] < 0 + + +def test_response_uses_the_net_of_tax_form(monkeypatch): + response = response_with(monkeypatch, net_of_tax_elasticity=0.5) + baseline_share = 1 - MEASUREMENTS["baseline_capital_gains_mtr"] + reform_share = 1 - MEASUREMENTS["reform_capital_gains_mtr"] + assert np.allclose( + response, GAINS * ((reform_share / baseline_share) ** 0.5 - 1), rtol=1e-5 + ) + + +def test_default_log_rate_response_matches_main_bit_for_bit(monkeypatch): + response = response_with(monkeypatch, elasticity=-0.7) + elasticity = np.full(3, -0.7, dtype=np.float32) + expected = GAINS * main_log_rate_response_factor( + MEASUREMENTS["baseline_capital_gains_mtr"], + MEASUREMENTS["reform_capital_gains_mtr"], + elasticity, + ) + assert response.dtype == expected.dtype + assert np.array_equal(bits(response), bits(expected)) + + +def test_response_rejects_person_level_parameters_outside_the_domain(monkeypatch): + monkeypatch.setattr( + capital_gains_module, + "get_behavioral_response_measurements", + lambda person, period: MEASUREMENTS, + ) + person = FakePerson( + { + "long_term_capital_gains_before_response": GAINS, + "capital_gains_semi_elasticity": np.array( + [3.4, -3.4, 3.4], dtype=np.float32 + ), + } + ) + with pytest.raises(ValueError, match="must be between 0 and 50"): + capital_gains_behavioral_response.formula( + person, 2026, make_parameters(semi_elasticity=3.4) + ) + + +@pytest.mark.parametrize( + "values", [{"semi_elasticity": 3.4}, {"net_of_tax_elasticity": 0.7}] +) +def test_new_forms_raise_on_a_rate_that_is_not_finite(monkeypatch, values): + measurements = { + "baseline_capital_gains_mtr": np.array([0.2, np.nan, 0.2], dtype=np.float32), + "reform_capital_gains_mtr": np.array([0.3, 0.3, 0.3], dtype=np.float32), + } + with pytest.raises(ValueError, match="not finite for 1 people"): + response_with(monkeypatch, measurements, **values) + + +def test_log_rate_form_still_passes_a_rate_that_is_not_finite_through(monkeypatch): + """The default path is unchanged, including its silence on a NaN rate.""" + measurements = { + "baseline_capital_gains_mtr": np.array([0.2, np.nan, 0.2], dtype=np.float32), + "reform_capital_gains_mtr": np.array([0.3, 0.3, 0.3], dtype=np.float32), + } + response = response_with(monkeypatch, measurements, elasticity=-0.7) + assert np.isnan(response[1]) and np.isfinite(response[[0, 2]]).all() + + +def test_response_is_zero_without_measuring_when_no_form_is_selected(monkeypatch): + def fail(person, period): + raise AssertionError("measured although no response form is selected") + + monkeypatch.setattr( + capital_gains_module, "get_behavioral_response_measurements", fail + ) + person = FakePerson({}) + assert ( + capital_gains_behavioral_response.formula(person, 2026, make_parameters()) == 0 + ) + + +def test_response_raises_when_two_forms_are_selected(monkeypatch): + with pytest.raises(ValueError, match="at most one"): + response_with(monkeypatch, elasticity=-0.7, semi_elasticity=3.4) + + +def test_response_is_zero_without_a_baseline(): + person = FakePerson({}) + person.simulation.baseline = None + parameters = make_parameters(semi_elasticity=3.4) + assert capital_gains_behavioral_response.formula(person, 2026, parameters) == 0 + + +def test_form_parameter_variables_read_their_parameters(): + parameters = make_parameters(semi_elasticity=3.4, net_of_tax_elasticity=0.7) + person = FakePerson({}) + assert capital_gains_semi_elasticity.formula(person, 2026, parameters) == 3.4 + assert capital_gains_net_of_tax_elasticity.formula(person, 2026, parameters) == 0.7 + + +# Properties. + + +@given(form_and_parameter=FORM_AND_PARAMETER, rate=WIDE_RATES) +def test_no_rate_change_means_no_response(form_and_parameter, rate): + form, parameter = form_and_parameter + assert factor(form, rate, rate, parameter) == 0 + + +@given( + form_and_parameter=FORM_AND_PARAMETER, + baseline=WIDE_RATES, + low=WIDE_RATES, + high=WIDE_RATES, +) +def test_realizations_do_not_rise_with_the_rate( + form_and_parameter, baseline, low, high +): + form, parameter = form_and_parameter + low, high = sorted((low, high)) + assert factor(form, baseline, low, parameter) >= factor( + form, baseline, high, parameter + ) + # Equivalently, realizations do not fall as the baseline rate rises. + assert factor(form, low, baseline, parameter) <= factor( + form, high, baseline, parameter + ) + + +@given( + form=st.sampled_from(CAPITAL_GAINS_RESPONSE_FORMS), + magnitude=st.floats(min_value=0.05, max_value=5), + baseline=RATES, + low=st.floats(min_value=0, max_value=1), + gap=st.floats(min_value=1e-4, max_value=1), +) +def test_realizations_fall_strictly_inside_the_unclipped_range( + form, magnitude, baseline, low, gap +): + lower_bound, upper_bound = UNCLIPPED_RATE_RANGE[form] + low = lower_bound + low * (upper_bound - lower_bound) + high = low + gap + assume(high <= upper_bound) + parameter = -magnitude if form == LOG_RATE_ELASTICITY else magnitude + at_low = factor(form, baseline, low, parameter) + at_high = factor(form, baseline, high, parameter) + # Near f = -1 a small rate step moves 1 + f by less than f's rounding. + assume(1 + at_high > 1e-6) + assert at_low > at_high + + +@given( + form_and_parameter=FORM_AND_PARAMETER, + baseline=WIDE_RATES, + reform=WIDE_RATES, + gains=st.floats(min_value=-1e10, max_value=1e10), +) +def test_response_never_flips_the_sign_of_gains( + form_and_parameter, baseline, reform, gains +): + form, parameter = form_and_parameter + change = factor(form, baseline, reform, parameter) + assert 1 + change >= 0 + assert gains * (1 + change) * np.sign(gains) >= 0 + + +@given( + form_and_parameter=FORM_AND_PARAMETER, + first=WIDE_RATES, + second=WIDE_RATES, + third=WIDE_RATES, +) +def test_response_factors_compose_along_a_path_of_rates( + form_and_parameter, first, second, third +): + form, parameter = form_and_parameter + legs = ( + 1 + factor(form, first, second, parameter), + 1 + factor(form, second, third, parameter), + ) + direct = 1 + factor(form, first, third, parameter) + # 1 + f is accurate to about 1e-16 in absolute terms, so when it is tiny a + # large second leg magnifies that error; keep factors that carry at least + # ten significant digits. + assume(min(*legs, direct) > 1e-6) + assert legs[0] * legs[1] == pytest.approx(direct, rel=1e-9) + + +@given(semi_elasticity=SEMI_ELASTICITIES, baseline=WIDE_RATES, reform=WIDE_RATES) +def test_semi_elasticity_factor_is_bounded(semi_elasticity, baseline, reform): + change = factor(SEMI_ELASTICITY, baseline, reform, semi_elasticity) + assert np.expm1(-semi_elasticity) <= change <= np.expm1(semi_elasticity) + + +@given( + semi_elasticity=st.floats(min_value=1.05, max_value=50), + baseline=RATES, + rate=RATES, + other_rate=RATES, +) +# Multipliers near 0 that rounding leaves flat or reversed (review of #10054). +@example(semi_elasticity=50, baseline=0, rate=0.6305, other_rate=0.6306) +def test_semi_elasticity_revenue_peaks_at_the_inverse_semi_elasticity( + semi_elasticity, baseline, rate, other_rate +): + peak = 1 / semi_elasticity + peak_revenue = revenue(SEMI_ELASTICITY, baseline, peak, semi_elasticity) + assert peak_revenue > 0 + assert at_most( + revenue(SEMI_ELASTICITY, baseline, rate, semi_elasticity), peak_revenue + ) + # Unimodal: rising below the peak and falling above it. + low, high = sorted((rate, other_rate)) + low_revenue = revenue(SEMI_ELASTICITY, baseline, low, semi_elasticity) + high_revenue = revenue(SEMI_ELASTICITY, baseline, high, semi_elasticity) + if high <= peak: + assert at_most(low_revenue, high_revenue) + if low >= peak: + assert at_most(high_revenue, low_revenue) + + +@given( + net_of_tax_elasticity=st.floats(min_value=0.003, max_value=20), + baseline=st.floats(min_value=0, max_value=0.99), + rate=RATES, + other_rate=RATES, +) +@example(net_of_tax_elasticity=20, baseline=0, rate=0.8, other_rate=0.8001) +def test_net_of_tax_revenue_peaks_at_one_over_one_plus_the_elasticity( + net_of_tax_elasticity, baseline, rate, other_rate +): + peak = 1 / (1 + net_of_tax_elasticity) + peak_revenue = revenue(NET_OF_TAX_ELASTICITY, baseline, peak, net_of_tax_elasticity) + assert at_most( + revenue(NET_OF_TAX_ELASTICITY, baseline, rate, net_of_tax_elasticity), + peak_revenue, + ) + low, high = sorted((rate, other_rate)) + low_revenue = revenue(NET_OF_TAX_ELASTICITY, baseline, low, net_of_tax_elasticity) + high_revenue = revenue(NET_OF_TAX_ELASTICITY, baseline, high, net_of_tax_elasticity) + if high <= peak: + assert at_most(low_revenue, high_revenue) + # Above the 99.9% cap realizations stop falling and revenue rises again. + if peak <= low and high <= 1 - NET_OF_TAX_SHARE_FLOOR: + assert at_most(high_revenue, low_revenue) + + +def test_net_of_tax_cap_puts_revenue_at_100_percent_above_a_tiny_elasticitys_peak(): + """Intended: below an elasticity of about 0.0027 the cap moves the maximum. + + At 0.002 the analytical peak is a 99.8% rate, and revenue at 100%, where + the cap holds realizations at their 99.9% level, is higher. At 0.003 the + analytical peak is higher. + """ + for elasticity, cap_wins in ((0.002, True), (0.003, False)): + peak = 1 / (1 + elasticity) + at_peak = revenue(NET_OF_TAX_ELASTICITY, 0.0, peak, elasticity) + at_100 = revenue(NET_OF_TAX_ELASTICITY, 0.0, 1.0, elasticity) + assert (at_100 > at_peak) is cap_wins + + +@given( + elasticity=st.floats(min_value=-0.99, max_value=0), + baseline=RATES, + rate=st.floats(min_value=0.001, max_value=1), + other_rate=st.floats(min_value=0.001, max_value=1), +) +def test_log_rate_revenue_has_no_interior_peak_above_minus_one( + elasticity, baseline, rate, other_rate +): + """The problem the semi-elasticity form solves: revenue keeps rising.""" + low, high = sorted((rate, other_rate)) + assert at_most( + revenue(LOG_RATE_ELASTICITY, baseline, low, elasticity), + revenue(LOG_RATE_ELASTICITY, baseline, high, elasticity), + ) + + +# No deadline: a wall-clock limit is not part of this compatibility property, +# and generating float32 arrays can exceed Hypothesis's 200 ms default on a +# loaded runner (it did once in review, 253 ms, then 86 ms on replay). +@settings(max_examples=300, deadline=None) +@given( + dtype=st.sampled_from([np.float32, np.float64]), + data=st.data(), +) +def test_log_rate_form_matches_main_bit_for_bit(dtype, data): + shape = data.draw(st.integers(min_value=1, max_value=20)) + rate_elements = st.floats( + min_value=-2, + max_value=5, + allow_nan=False, + width=32 if dtype is np.float32 else 64, + ) + baseline = data.draw(hnp.arrays(dtype, shape, elements=rate_elements)) + reform = data.draw(hnp.arrays(dtype, shape, elements=rate_elements)) + elasticity = data.draw( + hnp.arrays( + dtype, + shape, + elements=st.floats( + min_value=-3, + max_value=0, + allow_nan=False, + width=32 if dtype is np.float32 else 64, + ), + ) + ) + new = log_rate_elasticity_response_factor(baseline, reform, elasticity) + old = main_log_rate_response_factor(baseline, reform, elasticity) + assert new.dtype == old.dtype + assert np.array_equal(bits(new), bits(old)) + + +FLOAT32_RATES = st.floats(min_value=-5, max_value=20, width=32) +# 2 ** 53, about $9 quadrillion, is exact in float32 and inside the documented +# $10 quadrillion limit. +FLOAT32_GAINS = st.floats(min_value=-(2.0**53), max_value=2.0**53, width=32) + + +@settings(deadline=None) +@given( + form=st.sampled_from([SEMI_ELASTICITY, NET_OF_TAX_ELASTICITY]), + data=st.data(), +) +def test_new_forms_give_finite_sign_preserving_float32_responses(form, data): + """Across each new form's domain, as the model computes it in float32.""" + size = data.draw(st.integers(min_value=1, max_value=20)) + parameter = data.draw( + hnp.arrays( + np.float32, + size, + elements=st.floats( + min_value=0, + max_value=CAPITAL_GAINS_RESPONSE_MAXIMUM[form], + width=32, + ), + ) + ) + baseline = data.draw(hnp.arrays(np.float32, size, elements=FLOAT32_RATES)) + reform = data.draw(hnp.arrays(np.float32, size, elements=FLOAT32_RATES)) + gains = data.draw(hnp.arrays(np.float32, size, elements=FLOAT32_GAINS)) + response = gains * capital_gains_response_factor(form, baseline, reform, parameter) + assert response.dtype == np.float32 + assert np.isfinite(response).all() + after = gains + response + assert (after * np.sign(gains) >= 0).all() + + +@pytest.mark.parametrize( + "dtype, semi_elasticity, round_trip, absolute_error", + [ + # At beta = 3.4 the exponents and +1 additions are exact. Allowing + # two ulps per expm1 gives errors of eps and 32*eps in the outbound + # and return factors. The leading product error is bounded by + # (exp(beta) + 32*exp(-beta) + 0.5)*eps, including final rounding. + # With second-order terms this is < 32*eps, across library builds. + (np.float32, 3.4, 1.0, 32 * np.finfo(np.float32).eps), + # Intended: expm1(-20) rounds to -1 in float32, so adding 1 gives + # exact zero and the return leg cannot undo it. + (np.float32, 20, 0.0, 0.0), + # In float64, two outbound ulps amplified by exp(20) contribute + # exp(20)*eps64 ~= 1.08e-7; the other rounding terms are negligible. + (np.float64, 20, 1.0, 2e-7), + ], +) +def test_round_trip_precision_matches_the_documentation( + dtype, semi_elasticity, round_trip, absolute_error +): + zero, one, beta = dtype(0), dtype(1), dtype(semi_elasticity) + out = 1 + semi_elasticity_response_factor(zero, one, beta) + back = 1 + semi_elasticity_response_factor(one, zero, beta) + assert float(out * back) == pytest.approx(round_trip, rel=0, abs=absolute_error) + + +def test_documented_float32_round_trip_bound_rejects_degraded_precision(): + """Rounding the outbound response to float16 must fail the float32 bound.""" + zero, one, beta = np.float32(0), np.float32(1), np.float32(3.4) + response = semi_elasticity_response_factor(zero, one, beta) + out = 1 + np.float32(np.float16(response)) + back = 1 + semi_elasticity_response_factor(one, zero, beta) + assert float(out * back) != pytest.approx( + 1.0, rel=0, abs=32 * np.finfo(np.float32).eps + ) diff --git a/policyengine_us/variables/gov/simulation/capital_gains_response_forms.py b/policyengine_us/variables/gov/simulation/capital_gains_response_forms.py new file mode 100644 index 00000000000..57485c9e2e8 --- /dev/null +++ b/policyengine_us/variables/gov/simulation/capital_gains_response_forms.py @@ -0,0 +1,194 @@ +"""Functional forms for the capital gains realization response. + +Each form maps a person's baseline and reform marginal tax rates on capital +gains (``tau0`` and ``tau1``) to a response factor ``f``. The behavioral +response is ``long_term_capital_gains_before_response * f``, so realizations +after the response are the pre-response gains times the multiplier ``1 + f``. +Every form has ``1 + f >= 0``, so the response never changes the sign of a +person's gains, and ``f = 0`` when ``tau1 == tau0``. For a net loss the +multiplier scales the loss, so a smaller multiplier moves it toward zero. + +The form is chosen by which of three parameters under +``gov.simulation.capital_gains_responses`` is nonzero: + +- ``elasticity`` (log-rate form, the original): + ``1 + f = exp(elasticity * [ln tau1 - ln tau0])``, with both rates floored + at 0.001 by ``calculate_relative_capital_gains_mtr_change``. In a single-rate + model, revenue ``tau * R`` is proportional to ``tau ** (1 + elasticity)``, so + for an elasticity between -1 and 0 it rises all the way to a 100% rate and + the form has no interior revenue-maximizing rate. Because of the floor, a + person whose baseline rate is 0% gets a log change of ``ln(tau1 / 0.001)``: + about 5.0 for a move to 15%, which removes about 97% of their realizations + at an elasticity of -0.7. +- ``semi_elasticity``: ``1 + f = exp(-semi_elasticity * (tau1 - tau0))``, with + both rates clipped to [0, 1]. It is defined at ``tau0 = 0`` with no floor. + In a single-rate model, revenue ``tau * R`` peaks at + ``tau = 1 / semi_elasticity`` when that is below 1 (a semi-elasticity above + 1); otherwise it rises to a 100% rate. +- ``net_of_tax_elasticity``: + ``1 + f = exp(net_of_tax_elasticity * [ln(1 - tau1) - ln(1 - tau0)])``, with + both rates clipped to [0, 1 - NET_OF_TAX_SHARE_FLOOR]. In a single-rate + model, revenue peaks at ``tau = 1 / (1 + net_of_tax_elasticity)`` for a + net-of-tax elasticity of at least 0.003. Above the 99.9% cap realizations + stop falling, so revenue rises again at the very top, and for elasticities + below about 0.0027, whose analytical peak is above 99.7%, revenue at 100% + is higher. + The floor on the net-of-tax share mirrors the log-rate form's floor on the + rate: a baseline rate at or above 99.9% gives a large log change for a + reform rate well below 99.9%. + +A progressive tax schedule, other taxes and benefit interactions move the +model's revenue-maximizing statutory rate away from these single-rate peaks. + +The rates are household finite-difference marginal rates +(``marginal_tax_rate_on_capital_gains``), which can fall outside [0, 1] when the +extra gain crosses a benefit cliff or phase-in. The semi-elasticity and +net-of-tax forms clip them to the range a tax rate can take, which bounds the +semi-elasticity multiplier ``1 + f`` to +[exp(-semi_elasticity), exp(semi_elasticity)]. + +The new forms accept only positive parameters, the sign under which +realizations fall as the rate rises, up to the maxima below; zero leaves the +form unselected. That rejects the likeliest slip, entering a semi-elasticity +with the negative sign the log-rate elasticity uses. The log-rate form keeps +its original arithmetic and accepts any finite elasticity, so scores that use +it do not change. +""" + +import numpy as np + +from policyengine_us.variables.gov.simulation.behavioral_response_measurements import ( + calculate_relative_capital_gains_mtr_change, +) + +LOG_RATE_ELASTICITY = "elasticity" +SEMI_ELASTICITY = "semi_elasticity" +NET_OF_TAX_ELASTICITY = "net_of_tax_elasticity" +CAPITAL_GAINS_RESPONSE_FORMS = ( + LOG_RATE_ELASTICITY, + SEMI_ELASTICITY, + NET_OF_TAX_ELASTICITY, +) +# The person-level variable that carries each form's parameter. +CAPITAL_GAINS_RESPONSE_FORM_VARIABLES = { + LOG_RATE_ELASTICITY: "capital_gains_elasticity", + SEMI_ELASTICITY: "capital_gains_semi_elasticity", + NET_OF_TAX_ELASTICITY: "capital_gains_net_of_tax_elasticity", +} +NET_OF_TAX_SHARE_FLOOR = 0.001 +# Largest parameter each new form accepts. At these limits the multiplier +# 1 + f is at most exp(50), about 5.2e21 (a semi-elasticity of 50 and a rate +# falling from 100% to 0%), or 1000 ** 5 = 1e15 (a net-of-tax elasticity of 5 +# and a net-of-tax share rising from 0.001 to 1), so the response to gains +# below $10 quadrillion stays finite in float32. The semi-elasticities in CRS +# R41364's Table 1 run from 1.0 to 4.1. +CAPITAL_GAINS_RESPONSE_MAXIMUM = { + SEMI_ELASTICITY: 50, + NET_OF_TAX_ELASTICITY: 5, +} + + +def check_capital_gains_response_parameter(form, values): + """Raise ValueError unless ``values`` are in the form's supported domain. + + Every form needs finite values. The semi-elasticity and net-of-tax forms + also need values from 0 to their maximum, 0 meaning no response. + """ + values = np.asarray(values, dtype=np.float64) + name = f"gov.simulation.capital_gains_responses.{form}" + if not np.all(np.isfinite(values)): + raise ValueError(f"{name} must be finite; got {values[~np.isfinite(values)]}.") + if form not in CAPITAL_GAINS_RESPONSE_MAXIMUM: + return + maximum = CAPITAL_GAINS_RESPONSE_MAXIMUM[form] + outside = (values < 0) | (values > maximum) + if np.any(outside): + raise ValueError( + f"{name} must be between 0 and {maximum}; got {values[outside]}. It is " + "positive when realizations fall as the rate rises, unlike the " + "negative log-rate elasticity: a semi-elasticity of 3.4 cuts " + "realizations by about 3.4% per percentage point." + ) + + +def selected_capital_gains_response_form(capital_gains_responses): + """The form whose parameter is nonzero, or None when all are zero. + + Raises ValueError when more than one is nonzero, because the forms are + alternatives and combining them would apply the response more than once, + and when the selected parameter is outside its supported domain. + """ + values = { + form: getattr(capital_gains_responses, form) + for form in CAPITAL_GAINS_RESPONSE_FORMS + } + nonzero = {form: value for form, value in values.items() if value != 0} + if len(nonzero) > 1: + settings = ", ".join(f"{form}={value}" for form, value in nonzero.items()) + raise ValueError( + "Set at most one of gov.simulation.capital_gains_responses." + f"{{{', '.join(CAPITAL_GAINS_RESPONSE_FORMS)}}} to a nonzero value; " + f"got {settings}. Each selects a different functional form for the " + "capital gains realization response." + ) + if not nonzero: + return None + form, value = next(iter(nonzero.items())) + check_capital_gains_response_parameter(form, value) + return form + + +def log_rate_elasticity_response_factor(baseline_rate, reform_rate, elasticity): + """Response factor for realizations proportional to rate ** elasticity.""" + log_rate_change = calculate_relative_capital_gains_mtr_change( + { + "baseline_capital_gains_mtr": baseline_rate, + "reform_capital_gains_mtr": reform_rate, + } + ) + return np.exp(elasticity * log_rate_change) - 1 + + +def semi_elasticity_response_factor(baseline_rate, reform_rate, semi_elasticity): + """Response factor for realizations proportional to exp(-semi_elasticity * rate).""" + baseline = np.clip(baseline_rate, 0, 1) + reform = np.clip(reform_rate, 0, 1) + return np.expm1(-semi_elasticity * (reform - baseline)) + + +def net_of_tax_elasticity_response_factor(baseline_rate, reform_rate, elasticity): + """Response factor for realizations proportional to (1 - rate) ** elasticity.""" + max_rate = 1 - NET_OF_TAX_SHARE_FLOOR + baseline_share = 1 - np.clip(baseline_rate, 0, max_rate) + reform_share = 1 - np.clip(reform_rate, 0, max_rate) + return np.expm1(elasticity * (np.log(reform_share) - np.log(baseline_share))) + + +CAPITAL_GAINS_RESPONSE_FACTORS = { + LOG_RATE_ELASTICITY: log_rate_elasticity_response_factor, + SEMI_ELASTICITY: semi_elasticity_response_factor, + NET_OF_TAX_ELASTICITY: net_of_tax_elasticity_response_factor, +} + + +def capital_gains_response_factor(form, baseline_rate, reform_rate, parameter): + return CAPITAL_GAINS_RESPONSE_FACTORS[form](baseline_rate, reform_rate, parameter) + + +def check_capital_gains_response_is_finite(form, response): + """Raise ValueError if a new form's response is NaN or infinite. + + With parameters in their supported domain this happens only when a measured + rate or a person's gains is not finite, or gains exceed $10 quadrillion. + The log-rate form is left unchecked, as it was before the other forms were + added. + """ + if form == LOG_RATE_ELASTICITY: + return + finite = np.isfinite(response) + if not np.all(finite): + raise ValueError( + f"The {form} capital gains response is not finite for " + f"{np.count_nonzero(~finite)} people; check their measured " + "marginal_tax_rate_on_capital_gains and pre-response gains." + ) diff --git a/policyengine_us/variables/gov/simulation/capital_gains_responses.py b/policyengine_us/variables/gov/simulation/capital_gains_responses.py index bc9bfc51d2c..f9b47d1aa0c 100644 --- a/policyengine_us/variables/gov/simulation/capital_gains_responses.py +++ b/policyengine_us/variables/gov/simulation/capital_gains_responses.py @@ -6,6 +6,13 @@ calculate_relative_capital_gains_mtr_change, get_behavioral_response_measurements, ) +from policyengine_us.variables.gov.simulation.capital_gains_response_forms import ( + CAPITAL_GAINS_RESPONSE_FORM_VARIABLES, + capital_gains_response_factor, + check_capital_gains_response_is_finite, + check_capital_gains_response_parameter, + selected_capital_gains_response_form, +) class relative_capital_gains_mtr_change(Variable): @@ -32,10 +39,40 @@ def formula(person, period, parameters): return gov.simulation.capital_gains_responses.elasticity +class capital_gains_semi_elasticity(Variable): + value_type = float + entity = Person + label = "semi-elasticity of capital gains realizations" + unit = "/1" + definition_period = YEAR + + def formula(person, period, parameters): + gov = parameters(period).gov + return gov.simulation.capital_gains_responses.semi_elasticity + + +class capital_gains_net_of_tax_elasticity(Variable): + value_type = float + entity = Person + label = "net-of-tax elasticity of capital gains realizations" + unit = "/1" + definition_period = YEAR + + def formula(person, period, parameters): + gov = parameters(period).gov + return gov.simulation.capital_gains_responses.net_of_tax_elasticity + + class capital_gains_behavioral_response(Variable): value_type = float entity = Person label = "capital gains behavioral response" + documentation = ( + "Change in long-term capital gains realizations from the change in the" + " person's capital gains marginal tax rate. The functional form is the" + " one whose parameter under gov.simulation.capital_gains_responses is" + " nonzero; see capital_gains_response_forms.py." + ) unit = USD definition_period = YEAR @@ -44,18 +81,24 @@ def formula(person, period, parameters): if simulation.baseline is None: return 0 - if parameters(period).gov.simulation.capital_gains_responses.elasticity == 0: + form = selected_capital_gains_response_form( + parameters(period).gov.simulation.capital_gains_responses + ) + if form is None: return 0 capital_gains = person("long_term_capital_gains_before_response", period) measurements = get_behavioral_response_measurements(person, period) - tax_rate_change = calculate_relative_capital_gains_mtr_change(measurements) - elasticity = person("capital_gains_elasticity", period) - - # Calculate response using log differences - response_factor = np.exp(elasticity * tax_rate_change) - 1 + parameter = person(CAPITAL_GAINS_RESPONSE_FORM_VARIABLES[form], period) + check_capital_gains_response_parameter(form, parameter) + response_factor = capital_gains_response_factor( + form, + measurements["baseline_capital_gains_mtr"], + measurements["reform_capital_gains_mtr"], + parameter, + ) response = capital_gains * response_factor - + check_capital_gains_response_is_finite(form, response) return response