[MaxGhenis]

Follows #10663. North Dakota counts a renter's income "from all sources" for its renter's refund but excludes "any amount excluded from income by federal or state law with the exception of income from social security benefits" (N.D.C.C. 57-02-08.1(5)(c)). On main, nd_renters_refund_income sums market_income, which counts wages before payroll deductions and the tax-exempt parts of interest, pensions, retirement distributions and alimony in full. So $5,000 of municipal bond interest raises a renter's income by $5,000. Main subtracts only pre-tax health premiums and, since #10663, disability benefits from a plan the recipient paid for.

What changes

nd_renters_refund_income now subtracts every part of market_income that federal law leaves out of gross income, in two pieces, and floors the result at zero as before.

nd_renters_refund_federally_excluded_income (new, person) adds up the excluded part of each of a person's market_income components. It takes each part at zero or more and at most the component it comes from, so an input that reports only a component's total keeps all of it as income.

Component	Excluded part	Law
employment_income	the payroll deductions in pre_tax_contributions: traditional 401(k) and 403(b) deferrals, pre-tax health premiums, HSA payroll contributions, and dependent care assistance up to the IRC 129 limit (employment_income less irs_employment_income and taxable_dependent_care_benefits)	IRC 106(a) and (d), 125(a), 129(a), 402(e)(3), 402(g)(1), 403(b)(1)
interest_income	tax_exempt_interest_income	IRC 103(a)
pension_income	tax_exempt_pension_income	IRC 72(b)(1), 402(c)(1), 402A(d)(1)
retirement_distributions	tax_exempt_retirement_distributions	IRC 408(d), 408A(d)(1)
disability_benefits	tax_exempt_disability_benefits	IRC 104(a)(3)
alimony_income	alimony under an instrument executed after 2018 (alimony_income less taxable_alimony_income)	Pub. L. 115-97 §11051; IRS Topic 452

gov.states.nd.tax.property.renters_refund.excluded_deductions (new parameter) lists four federal exclusions that the federal model takes as above-the-line deductions, so market_income still holds the income they exclude. It names savings bond interest spent on higher education (us_bonds_for_higher_ed, IRC 135(a)), employer adoption assistance (qualified_adoption_assistance_expense, IRC 137(a)), and income from Guam, American Samoa, the Northern Mariana Islands and Puerto Rico (specified_possession_income and puerto_rico_income, IRC 931(a) and 933). Households may record the first two on a dependent (gov.irs.ald.filer_amounts_recorded_on_dependents), and the model takes the other two as tax unit inputs. So nd_renters_refund_income takes them out of the tax unit's market income left after the person-level parts, and never out of social security.

These replace main's separate subtractions of pre-tax premiums and excluded disability benefits, which the person-level variable now covers.

What stays in, and why

• These components count in full, as they do in the model's federal gross income: self-employment, partnership and S corporation income, dividends, capital gains, rental and farm income, estate income, illicit and miscellaneous income, military retirement pay, the Alaska Permanent Fund Dividend and other gains. Federal law excludes parts of some of them, such as gain on the sale of a main home (IRC 121) and some military disability retirement pay (IRC 104(a)(4)), but no model input holds those parts.
• The model records wage and self-employment income after the foreign earned income and housing exclusions (IRC 911). foreign_earned_income_exclusion documents that "Other income inputs are the amounts left after the exclusion", so nothing comes out for them.
• The federal model's gross income leaves out strike_benefits, but IRS Publication 525 counts strike benefits as income unless the union meant them as gifts.
• A guaranteed income payment (gi_cash_assistance) counts as a gift, which IRC 102 excludes, or as public assistance, which North Dakota counts, depending on who pays it, and the input doesn't record the payer.
• Social security counts in full, by the statute's own exception.
• Unemployment compensation counts in full: the only federal exclusion, ARPA's $10,200, applied to 2020, and the North Dakota parameters start in 2025.

Readings

This PR logs these readings in source_conflicts/candidates.yaml as us-nd:renters-refund-federally-excluded-income and us-nd:renters-refund-state-income-tax-reductions.

• The Tax Commissioner lists pensions "including military pension and I.R.A.", salaries and wages, dividends and interest, and alimony as income in a guideline (item 10). On the application (SFN 24777, lines 5 to 7), the Commissioner asks for income from salary and wages and from interest. Neither makes an exception for the parts federal law excludes. The statute names the exclusion, so it governs over both, as in us-nd:renters-refund-insurance-proceeds (Count every taxed cash source in SPM market income #10663).
• No section of IRC Part III, which lists the exclusions from gross income, names post-2018 alimony or elective deferrals, so these two make the closest calls. Congress repealed IRC 71 and 61(a)(8) in section 11051 of Pub. L. 115-97, and the IRS says the recipient's gross income doesn't include alimony under an instrument executed after 2018. Congress also treats elective deferrals as neither distributed nor made available to the employee (IRC 402(e)(3)). This PR reads "excluded from income" in its ordinary sense, as any amount federal law leaves out of income, and the statute's exception for social security, part of which IRC 86 leaves out of gross income, fits that reading.
• North Dakota computes taxable income by adjusting federal taxable income, mostly through amounts it lists after "Reduced by", such as interest on U.S. obligations, military pay, retired military benefits and 40 percent of net long-term capital gains and qualified dividends (N.D.C.C. 57-38-30.3(2)). The list also names plain deductions, such as organ donation expenses. The refund's income clause doesn't say whether those reductions exclude amounts from income. The Legislative Assembly added the clause's exception for social security in House Bill 1174 (2019), the Act that created the social security reduction, and the Tax Commissioner's office tied the two together in its fiscal note. The exception also keeps in income the part of benefits that IRC 86 leaves out of gross income, so it has effect even if the reductions exclude nothing. N.D.C.C. 1-02-39 lets a court weigh legislative history and administrative construction when a statute is ambiguous. The Commissioner, who administers the refund, kept "Pensions, including military pension" and "Net gains from the sale of property" in item 10 of the guideline after both 2019 Acts. So this PR adds no state-law subtraction, and a YAML case pins that result.

Invariants

Each holds for every input, and policyengine_us/tests/test_nd_renters_refund_income_exclusions.py tests it with Hypothesis.

1. Every source has a stated treatment. North Dakota counts each market_income component in full, subtracts the part federal law excludes, or counts it although the federal model's gross income leaves it out (with the reason). The test also sorts each above-the-line deduction in the federal model into a deduction, which leaves its income in, or an exclusion, which North Dakota subtracts. A new component or deduction fails the test until someone classifies it, so a wider market_income can't bring excluded income back, as Count every taxed cash source in SPM market income #10663's did.
2. Bounds. Each excluded part lies between zero and the component it comes from. The deducted exclusions come only out of the market income left after the others, so social security always counts in full.
3. Differential. For a filer whose income comes only from sources market_income and federal gross income share, renter's refund income equals irs_gross_income less the deducted exclusions.
4. Dollar for dollar. Moving a dollar from the taxable form of a source to its excluded form lowers renter's refund income by that dollar. A dollar of social security raises it by a dollar, however much of the benefit federal law taxes.

These tests use Python because they draw random amounts over many households at once and read lists from the source files, while YAML pins one household per case. Locally (macOS, Python 3.14, a shared host under load), the file runs 7 tests in 19 s (50 s wall with uv run), with a peak RSS of 1.0 GiB. It lands in the remaining Rest group.

Against main, each of the four tests that compute incomes fails. The differential test gets $74,000 on its explicit example, where the statute gives $42,000, and moving a dollar of taxable interest to tax-exempt interest leaves income unchanged. Alimony under a 2019 or 2024 instrument counts in full ($22,000, where the statute gives $10,000), and the bounds test stops because main has no excluded-income variable. Against this PR's first commit, which didn't yet subtract the four exclusions the model deducts above the line, the differential test gets $48,000. At that commit, the model still counted the example's $2,000 of savings bond interest and $4,000 of adoption assistance. The shift test fails there on savings bond interest, and the bounds and alimony tests pass.

The PR adds 13 YAML cases to nd_renters_refund_income.yaml and 5 to nd_renters_refund_federally_excluded_income.yaml, with expected values from the statute and the IRC sections each case cites.

2026 impact

We ran the microsimulation once at main (5e7e8f3) and once on this PR merged with that main (1eef8ca), each on the default dataset, Populace US 2024 (populace-us-2024-spm-20260909), for 2026. Every figure below comes from those two runs, and a comment below holds the scripts that ran and compared them. North Dakota pays the refund to tax units that rent and have a head or spouse aged 65 or over or permanently and totally disabled, if their income allows. We estimate 26,243 of these older or disabled renters in the state (69 sample records).

	Main	This PR	Change
North Dakota renter's refund	$5.73 million	$6.63 million	+$0.90 million
Tax units receiving it	10,745	12,625	+1,880
Renter's refund income of older or disabled renters	$3,276 million	$3,205 million	−$71.5 million
Household net income in North Dakota	$40,792.6 million	$40,793.5 million	+$0.90 million
SPM net income in North Dakota	$37,294.8 million	$37,295.7 million	+$0.90 million
People in SPM poverty in North Dakota	50,890	50,890	0

• The refund rises for 1,880 tax units (6 sample records) and falls for none. Renter's refund income falls for 9,689 older or disabled renters (31 records) and rises for none.
• Count property tax and rent refunds paid outside the income tax in household and SPM benefits #10725 counts the refund in household and SPM benefits, so household and SPM net income rise for the same 6 records by the same $0.90 million. No household or SPM unit outside North Dakota changes, and no SPM unit crosses its poverty threshold.
• One sample record, weighted at 1,878 tax units, accounts for 99.9% of the dollar change. The model now subtracts its $15,693 of nontaxable pension income, so its renter's refund income falls from $73,337 to $57,644, under the $70,000 limit, and its refund rises from zero to $477.
• Four records start receiving the refund. On main, two had income over the $70,000 limit, and the other two got nothing because 4% of their income exceeded 20% of their rent. Two records that already received it get $1 and $308 more, the second reaching the $600 cap.
• Among older or disabled renters, the model now subtracts nontaxable pensions ($54.5 million, 23 records), tax-exempt interest ($8.8 million, 4), traditional 401(k) deferrals ($5.3 million, 5) and nontaxable IRA distributions ($2.9 million, 4). None of them reports excluded disability benefits, post-2018 alimony, savings bond interest, adoption assistance, possession income or Puerto Rico income, or the pre-tax health premiums that main subtracted.
• These figures rest on few sample records (69 older or disabled renters, 6 that gain, and one that accounts for nearly all of the dollar change), so treat them as rough.

Checks run locally

At 1eef8ca, which merges main at 5e7e8f3, on macOS with Python 3.14:

• ruff format --check and ruff check on the three changed Python files: clean.
• policyengine_us/tests/test_nd_renters_refund_income_exclusions.py: 7 passed.
• policyengine-core test policyengine_us/tests/policy/baseline/gov/states/nd: 467 passed.
• policyengine_us/tests/test_household_spm_resource_parity.py: 68 passed.
• policyengine_us/tests/code_health: 533 passed, 1 skipped.

One of two runs of the parity file failed test_cash_income_never_lowers_spm_net_income, on the Connecticut draw that #10725 records as flaky on main. Main at 5e7e8f3 fails that draw with the same numbers.

Follow-ups

• The model leaves SSI and TANF out of renter's refund income (case 3 of nd_renters_refund_income.yaml keeps $6,000 of SSI out). The statute counts "any county, state, or federal public assistance benefits", and the application lists S.S.I. on line 7, while the guideline's item 11 lists TANF as not income.
• The federal gross income sources leave out strike_benefits, which IRS Publication 525 counts as income unless the union meant them as gifts.
• The federal model deducts us_bonds_for_higher_ed and qualified_adoption_assistance_expense in full, without the limits and income phase-outs of IRC 135(b) and 137(b), and nd_renters_refund_income subtracts the same full amounts (Savings bond (26 USC 135) and employer adoption assistance (26 USC 137) exclusions are deducted as supplied, with no cap or phase-out #10523).

axiom: TheAxiomFoundation/rulespec-us#2042 queued (North Dakota's income definition)

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (aae7f98) to head (1eef8ca).
⚠️ Report is 83 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff            @@
##              main    #10697   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            1         2    +1     
  Lines           37        40    +3     
=========================================
+ Hits            37        40    +3     

Flag	Coverage Δ
unittests	100.00% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[MaxGhenis]

These scripts produced the 2026 impact in the description. Each nd_impact.py run loads the default dataset, Populace US 2024 (populace-us-2024-spm-20260909), from a checkout of one commit:

PYTHONPATH=<checkout at 5e7e8f35fa> python nd_impact.py v3_base_2026.json 2026
PYTHONPATH=<checkout at 1eef8ca2f0> python nd_impact.py v3_head_2026.json 2026
python compare_impact.py v3_base_2026.json v3_head_2026.json 5e7e8f35fa 1eef8ca2f0 compare_v3_2026.json
python gainers.py v3_base_2026.json v3_head_2026.json

compare_impact.py checks that both runs measured the right commit on the same dataset, year, tax units and weights before it compares them. gainers.py lists each record whose refund changes, with the households and SPM units that change alongside it.

nd_impact.py

"""North Dakota's renter's refund on the whole pinned Populace file, for one tree.

Usage: PYTHONPATH=<policyengine-us tree> python nd_impact.py OUT.json YEAR

Writes OUT.json with North Dakota totals (refund, recipients, eligibility,
income, and the market income parts federal law excludes) and OUT.npz with
each tax unit's rows, so two runs can be compared tax unit by tax unit.
Every number comes from a Microsimulation of the default dataset.
"""

import json
import os
import subprocess
import sys
from pathlib import Path

out, year = Path(sys.argv[1]), int(sys.argv[2])
# One run per output: the first run to start claims it, and a duplicate queued
# for the same output exits at once, before it loads the model.
claim = out.with_suffix(".claimed")
try:
    fd = os.open(claim, os.O_CREAT | os.O_EXCL | os.O_WRONLY)
except FileExistsError:
    print(f"{out} already claimed: {claim.read_text().strip()}; exiting", flush=True)
    sys.exit(0)
os.write(fd, f"pid {os.getpid()} PYTHONPATH {os.environ.get('PYTHONPATH', '')}\n".encode())
os.close(fd)

import numpy as np

import policyengine_us
from policyengine_us import Microsimulation
from policyengine_us.system import DEFAULT_DATASET

sim = Microsimulation()
variables = sim.tax_benefit_system.variables
tree = Path(policyengine_us.__file__).parent.parent


def calc(name, map_to=None):
    return np.asarray(sim.calculate(name, year, map_to=map_to).values)


def unit_sum(person_values):
    """The tax unit's total of a person array."""
    return np.asarray(
        sim.map_result(np.asarray(person_values, dtype=float), "person", "tax_unit", how="sum"),
        dtype=float,
    )


person_weight = calc("person_weight").astype(float)
unit_weight = calc("tax_unit_weight").astype(float)
person_nd = calc("state_code_str", map_to="person").astype(str) == "ND"
unit_nd = unit_sum(person_nd) > 0

refund = calc("nd_renters_refund").astype(float)
income = calc("nd_renters_refund_income").astype(float)
eligible = calc("nd_renters_refund_eligible").astype(bool)
rent = unit_sum(calc("rent").astype(float))
age = calc("age").astype(float)
disabled = calc("is_permanently_and_totally_disabled").astype(bool)
head_or_spouse = calc("is_tax_unit_head_or_spouse").astype(bool)
p = sim.tax_benefit_system.parameters(f"{year}-01-01").gov.states.nd.tax.property.renters_refund
age_or_disability = unit_sum(((age >= p.age_threshold) | disabled) & head_or_spouse) > 0
# The population the refund can reach: ND tax units with a qualifying head or
# spouse who pay rent.
reachable = unit_nd & age_or_disability & (rent > 0)


def weighted(values, mask, weight=unit_weight):
    return float((np.asarray(values, dtype=float) * weight * mask).sum())


row = {
    "commit": subprocess.run(
        ["git", "-C", str(tree), "rev-parse", "HEAD"], capture_output=True, text=True
    ).stdout.strip(),
    "tree": str(tree),
    "dataset": DEFAULT_DATASET,
    "year": year,
    "nd_tax_units": weighted(np.ones_like(unit_weight), unit_nd),
    "nd_people": float((person_weight * person_nd).sum()),
    "reachable_tax_units": weighted(np.ones_like(unit_weight), reachable),
    "nd_renters_refund": weighted(refund, unit_nd),
    "recipients": weighted(refund > 0, unit_nd),
    "eligible": weighted(eligible, unit_nd),
    "nd_renters_refund_income": weighted(income, unit_nd),
    "nd_renters_refund_income_reachable": weighted(income, reachable),
    "income_limit": float(p.income_limit),
}
# The parts of market income federal law excludes, among ND people and among
# people in reachable tax units. Read from their inputs, so both trees report
# them.
person_reachable = np.asarray(sim.map_result(reachable.astype(float), "tax_unit", "person"), dtype=float) > 0
PERSON_TOTALS = (
    "market_income",
    "employment_income",
    "irs_employment_income",
    "pre_tax_contributions",
    "traditional_401k_contributions",
    "traditional_403b_contributions",
    "pre_tax_health_insurance_premiums",
    "health_savings_account_payroll_contributions",
    "dependent_care_employer_benefits",
    "taxable_dependent_care_benefits",
    "interest_income",
    "tax_exempt_interest_income",
    "pension_income",
    "tax_exempt_pension_income",
    "retirement_distributions",
    "tax_exempt_retirement_distributions",
    "disability_benefits",
    "tax_exempt_disability_benefits",
    "taxable_disability_benefits",
    "alimony_income",
    "taxable_alimony_income",
    "strike_benefits",
    "gi_cash_assistance",
    "social_security",
    "nd_renters_refund_federally_excluded_income",
    # Federal exclusions the model deducts above the line.
    "us_bonds_for_higher_ed",
    "qualified_adoption_assistance_expense",
)
# Tax unit inputs of the same kind.
UNIT_TOTALS = ("specified_possession_income", "puerto_rico_income")
person_rows = {}
for name in PERSON_TOTALS:
    if name not in variables:
        row[name] = None
        continue
    values = calc(name).astype(float)
    person_rows[name] = values
    row[name] = float((values * person_weight * person_nd).sum())
    row[name + "_reachable"] = float((values * person_weight * person_reachable).sum())
unit_inputs = {}
for name in UNIT_TOTALS:
    values = calc(name).astype(float)
    unit_inputs[name] = values
    row[name] = weighted(values, unit_nd)
    row[name + "_reachable"] = weighted(values, reachable)
divorce_year = calc("divorce_year")
row["alimony_after_2018_reachable"] = float(
    (person_rows["alimony_income"] * (divorce_year >= 2019) * person_weight * person_reachable).sum()
)
row["people_with_divorce_year_after_2018"] = float(((divorce_year >= 2019) * person_weight).sum())


def any_nd(entity):
    """Whether each row of an entity has a North Dakota member."""
    return (
        np.asarray(
            sim.map_result(person_nd.astype(float), "person", entity, how="sum"),
            dtype=float,
        )
        > 0
    )


# Net income, which counts the refund since #10725. These totals run the whole
# model on every record.
household_weight = calc("household_weight").astype(float)
household_nd = any_nd("household")
household_net = calc("household_net_income").astype(float)
spm_weight = calc("spm_unit_weight").astype(float)
spm_nd = any_nd("spm_unit")
spm_net = calc("spm_unit_net_income").astype(float)
spm_benefits = calc("spm_unit_benefits").astype(float)
spm_poor = calc("spm_unit_is_in_spm_poverty").astype(bool)
person_poor = calc("spm_unit_is_in_spm_poverty", map_to="person").astype(bool)
row["household_net_income_nd"] = float((household_net * household_weight * household_nd).sum())
row["household_net_income_us"] = float((household_net * household_weight).sum())
row["spm_unit_net_income_nd"] = float((spm_net * spm_weight * spm_nd).sum())
row["spm_unit_benefits_nd"] = float((spm_benefits * spm_weight * spm_nd).sum())
row["people_in_spm_poverty_nd"] = float((person_poor * person_weight * person_nd).sum())
row["people_in_spm_poverty_us"] = float((person_poor * person_weight).sum())
row["people_us"] = float(person_weight.sum())
np.savez_compressed(
    out.with_suffix(".net.npz"),
    household_id=calc("household_id").astype(np.int64),
    household_weight=household_weight,
    household_nd=household_nd,
    household_net=household_net,
    spm_unit_id=calc("spm_unit_id").astype(np.int64),
    spm_weight=spm_weight,
    spm_nd=spm_nd,
    spm_net=spm_net,
    spm_benefits=spm_benefits,
    spm_poor=spm_poor,
)
out.write_text(json.dumps(row, indent=1) + "\n")

# Tax unit rows for comparing two runs.
unit_rows = dict(
    tax_unit_id=calc("tax_unit_id").astype(np.int64),
    weight=unit_weight,
    nd=unit_nd,
    reachable=reachable,
    refund=refund,
    income=income,
    eligible=eligible,
    rent=rent,
)
for name, values in person_rows.items():
    unit_rows["unit_" + name] = unit_sum(values)
for name, values in unit_inputs.items():
    unit_rows["unit_" + name] = values
np.savez_compressed(out.with_suffix(".npz"), **unit_rows)
print(json.dumps(row, indent=1), flush=True)

compare_impact.py

"""Compare two nd_impact.py runs, tax unit by tax unit.

Usage: python compare_impact.py BASE.json HEAD.json BASE_SHA HEAD_SHA OUT.json

Checks that each run measured the commit it should, on the same dataset and
year, with the same tax units and weights. Then reports the change in North
Dakota's renter's refund and in renter's refund income, which tax units gain
or lose, and the market income components behind the income change. Every
number comes from the two runs' outputs; nothing is scaled or imputed.
"""

import json
import sys
from pathlib import Path

import numpy as np

base_json, head_json = Path(sys.argv[1]), Path(sys.argv[2])
want_base, want_head, out = sys.argv[3], sys.argv[4], Path(sys.argv[5])
base = json.loads(base_json.read_text())
head = json.loads(head_json.read_text())
assert base["commit"].startswith(want_base), ("base commit", base["commit"])
assert head["commit"].startswith(want_head), ("head commit", head["commit"])
assert base["dataset"] == head["dataset"], (base["dataset"], head["dataset"])
assert base["year"] == head["year"]

b = np.load(base_json.with_suffix(".npz"))
h = np.load(head_json.with_suffix(".npz"))
assert np.array_equal(b["tax_unit_id"], h["tax_unit_id"]), "tax units differ"
assert np.allclose(b["weight"], h["weight"]), "weights differ"
assert np.array_equal(b["nd"], h["nd"]), "ND tax units differ"
assert np.array_equal(b["reachable"], h["reachable"]), "reachable units differ"
w = b["weight"].astype(float)
nd = b["nd"].astype(bool)
reach = b["reachable"].astype(bool)


def weighted(values, mask):
    return float((np.asarray(values, dtype=float) * w * mask).sum())


def count(mask):
    return {"weighted": weighted(np.ones_like(w), mask), "records": int(mask.sum())}


refund_b, refund_h = b["refund"].astype(float), h["refund"].astype(float)
income_b, income_h = b["income"].astype(float), h["income"].astype(float)
gain = refund_h - refund_b
TOL = 0.5
up = nd & (gain > TOL)
down = nd & (gain < -TOL)
new = nd & (refund_b <= TOL) & (refund_h > TOL)
lost = nd & (refund_b > TOL) & (refund_h <= TOL)
income_down = nd & (income_h < income_b - TOL)
income_up = nd & (income_h > income_b + TOL)
outside = nd & ~reach & (np.abs(gain) > TOL)


def unit(name, arrays=h):
    key = "unit_" + name
    return arrays[key].astype(float) if key in arrays.files else None


# The head's excluded parts, by component, from the tax unit totals of the
# inputs. The head variable caps each part per person; these totals use the
# unit sums, so the script also reports the head variable's own total.
wages = unit("employment_income")
counted_wages = unit("irs_employment_income") + unit("taxable_dependent_care_benefits")
components = {
    "payroll deductions": wages - counted_wages,
    "tax-exempt interest": unit("tax_exempt_interest_income"),
    "nontaxable pensions": unit("tax_exempt_pension_income"),
    "nontaxable retirement distributions": unit("tax_exempt_retirement_distributions"),
    "disability benefits from a plan the recipient paid for": unit(
        "tax_exempt_disability_benefits"
    ),
    "alimony under a post-2018 instrument": unit("alimony_income")
    - unit("taxable_alimony_income"),
    "savings bond interest for higher education": unit("us_bonds_for_higher_ed"),
    "employer adoption assistance": unit("qualified_adoption_assistance_expense"),
    "possession income (IRC 931)": unit("specified_possession_income"),
    "Puerto Rico income (IRC 933)": unit("puerto_rico_income"),
}
# What main subtracted: pre-tax health premiums and excluded disability benefits.
main_subtractions = {
    "pre-tax health premiums": unit("pre_tax_health_insurance_premiums", b),
    "disability benefits from a plan the recipient paid for": unit(
        "disability_benefits", b
    )
    - unit("taxable_disability_benefits", b),
}


def by_component(table, mask):
    return {
        name: {
            "weighted_total": weighted(values, mask),
            "units_with_any": count(mask & (values > TOL)),
        }
        for name, values in table.items()
    }


report = {
    "dataset": base["dataset"],
    "year": base["year"],
    "base_commit": base["commit"],
    "head_commit": head["commit"],
    "nd_tax_units": count(nd),
    "reachable_tax_units": count(reach),
    "refund": {
        "base": weighted(refund_b, nd),
        "head": weighted(refund_h, nd),
        "change": weighted(gain, nd),
    },
    "recipients": {"base": count(nd & (refund_b > TOL)), "head": count(nd & (refund_h > TOL))},
    "eligible": {
        "base": count(nd & b["eligible"].astype(bool)),
        "head": count(nd & h["eligible"].astype(bool)),
    },
    "units_refund_up": count(up),
    "units_refund_changed_outside_reach": count(outside),
    "units_refund_down": count(down),
    "units_newly_receiving": count(new),
    "units_no_longer_receiving": count(lost),
    "gain_among_gainers": {
        "total": weighted(gain, up),
        "mean": weighted(gain, up) / max(weighted(np.ones_like(w), up), 1e-9),
        "max": float(gain[up].max()) if up.any() else 0.0,
    },
    "loss_among_losers": {"total": weighted(gain, down)},
    "income_reachable": {
        "base": weighted(income_b, reach),
        "head": weighted(income_h, reach),
        "change": weighted(income_h - income_b, reach),
    },
    "income_nd": {
        "base": weighted(income_b, nd),
        "head": weighted(income_h, nd),
        "change": weighted(income_h - income_b, nd),
    },
    "units_income_down": count(nd & income_down),
    "units_income_up": count(nd & income_up),
    "reachable_units_income_down": count(reach & income_down),
    "reachable_units_income_up": count(reach & income_up),
    "head_excluded_income_variable_reachable": head.get(
        "nd_renters_refund_federally_excluded_income_reachable"
    ),
    "head_exclusions_reachable": by_component(components, reach),
    "head_exclusions_among_refund_gainers": by_component(components, up),
    "main_subtractions_reachable": by_component(main_subtractions, reach),
    "main_subtractions_among_refund_losers": by_component(main_subtractions, down),
    "income_limit": head["income_limit"],
}

# Net income, which counts the refund since #10725.
bn = np.load(base_json.with_suffix(".net.npz"))
hn = np.load(head_json.with_suffix(".net.npz"))
for key in ("household_id", "spm_unit_id", "household_nd", "spm_nd"):
    assert np.array_equal(bn[key], hn[key]), key
for key in ("household_weight", "spm_weight"):
    assert np.allclose(bn[key], hn[key]), key
hw = bn["household_weight"].astype(float)
sw = bn["spm_weight"].astype(float)
hnd = bn["household_nd"].astype(bool)
snd = bn["spm_nd"].astype(bool)
household_change = hn["household_net"].astype(float) - bn["household_net"].astype(float)
spm_change = hn["spm_net"].astype(float) - bn["spm_net"].astype(float)
benefit_change = hn["spm_benefits"].astype(float) - bn["spm_benefits"].astype(float)
left_poverty = bn["spm_poor"].astype(bool) & ~hn["spm_poor"].astype(bool)
entered_poverty = ~bn["spm_poor"].astype(bool) & hn["spm_poor"].astype(bool)


def total(values, mask, weights):
    return float((np.asarray(values, dtype=float) * weights * mask).sum())


def cnt(mask, weights):
    return {"weighted": float((weights * mask).sum()), "records": int(mask.sum())}


report["net_income"] = {
    "household_net_income_nd": {
        "base": base["household_net_income_nd"],
        "head": head["household_net_income_nd"],
        "change": total(household_change, hnd, hw),
    },
    "household_net_income_us": {
        "base": base["household_net_income_us"],
        "head": head["household_net_income_us"],
        "change": total(household_change, np.ones_like(hw), hw),
    },
    "households_up_nd": cnt(hnd & (household_change > TOL), hw),
    "households_down_nd": cnt(hnd & (household_change < -TOL), hw),
    "largest_absolute_household_change_outside_nd": float(np.abs(household_change[~hnd]).max()),
    "refund_change_less_household_net_income_change_nd": weighted(gain, nd)
    - total(household_change, hnd, hw),
    "spm_unit_net_income_nd": {
        "base": base["spm_unit_net_income_nd"],
        "head": head["spm_unit_net_income_nd"],
        "change": total(spm_change, snd, sw),
    },
    "spm_unit_benefits_nd": {
        "base": base["spm_unit_benefits_nd"],
        "head": head["spm_unit_benefits_nd"],
        "change": total(benefit_change, snd, sw),
    },
    "spm_units_up_nd": cnt(snd & (spm_change > TOL), sw),
    "spm_units_down_nd": cnt(snd & (spm_change < -TOL), sw),
    "largest_absolute_spm_change_outside_nd": float(np.abs(spm_change[~snd]).max()),
    "spm_units_leaving_poverty": cnt(left_poverty, sw),
    "spm_units_entering_poverty": cnt(entered_poverty, sw),
    "people_in_spm_poverty_nd": {
        "base": base["people_in_spm_poverty_nd"],
        "head": head["people_in_spm_poverty_nd"],
        "change": head["people_in_spm_poverty_nd"] - base["people_in_spm_poverty_nd"],
    },
    "spm_poverty_rate_nd": {
        "base": base["people_in_spm_poverty_nd"] / base["nd_people"],
        "head": head["people_in_spm_poverty_nd"] / head["nd_people"],
    },
    "people_in_spm_poverty_us": {
        "base": base["people_in_spm_poverty_us"],
        "head": head["people_in_spm_poverty_us"],
        "change": head["people_in_spm_poverty_us"] - base["people_in_spm_poverty_us"],
    },
}
out.write_text(json.dumps(report, indent=1) + "\n")
print(json.dumps(report, indent=1))

gainers.py

"""List the tax units whose North Dakota renter's refund changes between two nd_impact.py runs.

Usage: python gainers.py BASE.json HEAD.json

Prints each gaining unit's weight, income, rent and refund on both runs and the
excluded amounts behind its income change. Then compares the changed households and
SPM units with the changed tax units, so a gap between the totals can be traced
to weights or to amounts. Reads only the two runs' outputs.
"""
import sys
from pathlib import Path

import numpy as np

base_json, head_json = Path(sys.argv[1]), Path(sys.argv[2])
b = np.load(base_json.with_suffix(".npz"))
h = np.load(head_json.with_suffix(".npz"))
bn = np.load(base_json.with_suffix(".net.npz"))
hn = np.load(head_json.with_suffix(".net.npz"))
w = b["weight"].astype(float)
nd = b["nd"].astype(bool)
gain = h["refund"].astype(float) - b["refund"].astype(float)
idx = np.where(nd & (np.abs(gain) > 0.005))[0]
idx = idx[np.argsort(-(gain[idx] * w[idx]))]
total = float((gain[idx] * w[idx]).sum())
print(f"total refund change among listed units: {total:,.2f}")
components = [
    "unit_tax_exempt_pension_income",
    "unit_tax_exempt_interest_income",
    "unit_tax_exempt_retirement_distributions",
    "unit_nd_renters_refund_federally_excluded_income",
]
for i in idx:
    extras = {c.replace("unit_", ""): round(float(h[c][i]), 2) for c in components if c in h.files}
    print(
        f"unit {int(b['tax_unit_id'][i])} weight {w[i]:,.4f} share {gain[i]*w[i]/total:7.3%} "
        f"income {b['income'][i]:,.2f}->{h['income'][i]:,.2f} rent {b['rent'][i]:,.2f} "
        f"refund {b['refund'][i]:,.2f}->{h['refund'][i]:,.2f} gain {gain[i]:,.2f} "
        f"4%inc {0.04*b['income'][i]:,.2f} 20%rent {0.2*b['rent'][i]:,.2f} {extras}"
    )
hw = bn["household_weight"].astype(float)
sw = bn["spm_weight"].astype(float)
dh = hn["household_net"].astype(float) - bn["household_net"].astype(float)
ds = hn["spm_net"].astype(float) - bn["spm_net"].astype(float)
db = hn["spm_benefits"].astype(float) - bn["spm_benefits"].astype(float)
print("--- households that change")
for i in np.where(np.abs(dh) > 0.005)[0]:
    print(f"household {int(bn['household_id'][i])} weight {hw[i]:,.4f} change {dh[i]:,.4f} weighted {dh[i]*hw[i]:,.2f}")
print("--- SPM units that change")
for i in np.where(np.abs(ds) > 0.005)[0]:
    print(f"spm {int(bn['spm_unit_id'][i])} weight {sw[i]:,.4f} net {ds[i]:,.4f} benefits {db[i]:,.4f} weighted {ds[i]*sw[i]:,.2f}")
print("--- sums")
print(f"refund   {(gain*w*nd).sum():,.2f}")
print(f"household {(dh*hw).sum():,.2f}")
print(f"spm net  {(ds*sw).sum():,.2f}")
print(f"spm ben  {(db*sw).sum():,.2f}")

[MaxGhenis]

#10725 counts nd_renters_refund in household and SPM benefits, so this change now moves both net incomes in North Dakota.