From 7968e033871988da3bb46047c0bc1ec83e6c8974 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 11 Jun 2026 15:29:31 +0200 Subject: [PATCH] populace-us-2024-9f1260b-20260611: eCPS-free build provenance Primary-source-only construction (CPS ASEC, IRS PUF, Fed SCF 2022, SIPP, CPS-ORG, MEPS-IC, ACS 2022); enhanced CPS benchmark-only. Gates: parity 0, exported-nonzero pass (309 cols), calibration 95.09% within 10%, smoke in band (STCG -$77.5B vs -$76.8B signed target). Score (matched 41,314, symmetric refit): train 0.190 vs 1.089, holdout 0.038 vs 0.317, full 0.228 vs 1.406. Card + manifest + chain snapshots. --- packages/populace-data/build/us/README.md | 31 +- .../populace-data/build/us/build_dataset.py | 255 +++ .../build/us/build_populace_us_dataset.py | 186 ++- .../build/us/build_us_candidate.py | 1400 +++++++++++++++++ .../populace-data/build/us/check_parity.py | 181 +++ .../populace-data/build/us/enrich_artifact.py | 150 ++ .../build/us/extract_target_surface.py | 51 +- .../populace-data/build/us/hf_dataset_card.md | 115 +- .../build/us/primary_source_impute.py | 601 +++++++ .../build/us/release_manifest.json | 58 + packages/populace-data/build/us/run_chain.sh | 30 + 11 files changed, 2885 insertions(+), 173 deletions(-) create mode 100644 packages/populace-data/build/us/build_dataset.py create mode 100644 packages/populace-data/build/us/build_us_candidate.py create mode 100644 packages/populace-data/build/us/check_parity.py create mode 100644 packages/populace-data/build/us/enrich_artifact.py create mode 100644 packages/populace-data/build/us/primary_source_impute.py create mode 100644 packages/populace-data/build/us/release_manifest.json create mode 100755 packages/populace-data/build/us/run_chain.sh diff --git a/packages/populace-data/build/us/README.md b/packages/populace-data/build/us/README.md index 7fbb2c86..26f527b1 100644 --- a/packages/populace-data/build/us/README.md +++ b/packages/populace-data/build/us/README.md @@ -1,18 +1,17 @@ -# build/us — provenance snapshot +# populace-us build provenance -The scripts that produced `populace_us_2024.h5`, kept verbatim for audit. They -are a **snapshot, not a turnkey pipeline**: they reference build-machine paths -(worktrees, cached surfaces, a scoring-harness checkout via -`SCORING_HARNESS_SRC`) and ran against `populace` installed from this -repository plus `policyengine-us==1.723.0`. +Verbatim snapshots of the scripts that produced the published artifact +(build `populace-us-2024-9f1260b-20260611`, HF revision `4a8e7d39eb9e`), +plus its release manifest. The chain (`run_chain.sh`) runs: full build +(`build_us_candidate.py`, with the primary-source imputation stages in +`primary_source_impute.py`) -> surface extraction +(`extract_target_surface.py`) -> calibration + artifact +(`build_dataset.py`) -> simulation-dependent enrichment +(`enrich_artifact.py`) -> acceptance gates (`check_parity.py`, using +`populace.build` gates). -- `extract_target_surface.py` — extracts the raw (unscaled) PE-native target - surface (3,704 targets) for the pool. -- `build_populace_us_dataset.py` — calibrates the pool's household weights to - that surface with `populace.calibrate` (hard `max_weight_ratio=50`) and - writes the published `USSingleYearDataset`, including the entity-table - surgeries it documents. -- `bounded_recal_experiment.py` + `bounded_recal_results.json` — the cap sweep - that selected the 50× bound. -- `hf_dataset_card.md` — the dataset card published to - `policyengine/populace-us` on the Hugging Face Hub. +Every donor is a primary source (CPS ASEC, IRS PUF, Fed SCF 2022, Census +SIPP, CPS-ORG, MEPS-IC parameters, Census ACS 2022); the enhanced CPS is +the scoring benchmark only. These are audit copies — the living code is +in the build worktree until the populace.build.us port lands, after which +the port is canonical and these snapshots freeze. diff --git a/packages/populace-data/build/us/build_dataset.py b/packages/populace-data/build/us/build_dataset.py new file mode 100644 index 00000000..1d134e7a --- /dev/null +++ b/packages/populace-data/build/us/build_dataset.py @@ -0,0 +1,255 @@ +"""Build the publishable populace-US v2 dataset: calibrate the v2 pool's +household weights to its raw PE-native surface (ratio-50 bound), write them +into the USSingleYearDataset copy and the timeperiod export, verify.""" + +import shutil +import sys +import time + +import h5py +import numpy as np +import pandas as pd + +from populace.calibrate import Target, TargetSet, calibrate +from populace.frame import EntitySchema, Frame, WeightKind, Weights + +ART = "/Users/maxghenis/.claude-worktrees/microplex-spec-build/artifacts" +POOL = f"{ART}/spec_candidate_full_2024/candidate_policyengine_us.h5" +TP = f"{ART}/spec_candidate_full_2024/candidate_timeperiod.h5" +SURFACE = f"{ART}/target_surface_raw.npz" +OUT = f"{ART}/populace_us_2024.h5" +OUT_TP = f"{ART}/populace_us_2024_timeperiod.h5" + + +def log(*a): + print(f"[{time.strftime('%H:%M:%S')}]", *a, flush=True) + + +def main(): + surf = np.load(SURFACE, allow_pickle=True) + A = surf["A"].astype(np.float64) + b = surf["b"].astype(np.float64) + w0 = surf["w0"].astype(np.float64) + names = [str(x) for x in surf["names"]] + n_hh, n_t = A.shape + log(f"{n_t} targets x {n_hh} households; w0 sum {w0.sum()/1e6:.1f}M") + + household = pd.DataFrame({"household_id": np.arange(n_hh, dtype=np.int64)}) + person = pd.DataFrame( + { + "person_id": np.arange(n_hh, dtype=np.int64), + "person_household_id": np.arange(n_hh, dtype=np.int64), + } + ) + frame = Frame( + {"person": person, "household": household}, + EntitySchema(group_entities=("household",)), + {"household": Weights(values=w0.copy(), kind=WeightKind.DESIGN)}, + ) + targets = TargetSet( + tuple( + Target( + name=names[t], + entity="household", + aggregation="sum", + value=float(b[t]), + measure=(lambda _f, col=A[:, t].copy(): col), + ) + for t in range(n_t) + ) + ) + # Signed heavy-tail targets (architecture review): calibration was free + # to amplify loss-heavy records in dimensions absent from the surface. + # Net short-term capital gains is anchored to the PUF donor's own + # weighted, uprated total — primary-source, computed, never hand-typed. + import h5py as _h5py + + _puf = pd.read_csv( + "/Users/maxghenis/.cache/microplex/puf_2015.csv", + usecols=["P22250", "S006"], + ) + _stcg_value = float( + ( + pd.to_numeric(_puf["P22250"], errors="coerce").fillna(0) + * pd.to_numeric(_puf["S006"], errors="coerce").fillna(0) + / 100.0 + ).sum() + * 1.8 # microplex puf.py uprating factor for short_term_capital_gains + ) + with _h5py.File(TP) as _f: + _stcg_p = _f["short_term_capital_gains"]["2024"][:].astype(np.float64) + _phh = _f["person_household_id"]["2024"][:] + _hid = _f["household_id"]["2024"][:] + _hidx = {h: i for i, h in enumerate(_hid.tolist())} + _stcg_hh = np.zeros(n_hh, dtype=np.float64) + np.add.at(_stcg_hh, np.fromiter((_hidx[h] for h in _phh.tolist()), dtype=np.int64), _stcg_p) + targets = TargetSet( + tuple(targets) + + ( + Target( + name="puf/net_short_term_capital_gains", + entity="household", + aggregation="sum", + value=_stcg_value, + measure=(lambda _f2, col=_stcg_hh.copy(): col), + ), + ) + ) + log(f"signed STCG target appended: ${_stcg_value/1e9:.1f}B (PUF weighted, uprated)") + + log("calibrating (ratio-50 bound)...") + t0 = time.time() + result = calibrate( + frame, + targets, + weight_entity="household", + epochs=3000, + learning_rate=0.15, + mass="free", + max_weight_ratio=50.0, + seed=0, + ) + cw = result.frame.resolve_weights("household").values.astype(np.float64) + # Telemetry (fail-soft): per-target diagnostics power the observatory's + # live fit tables and cross-run regression checks. + try: + import pathlib as _pathlib + import sys as _sys + + _sys.path.insert(0, str(_pathlib.Path(__file__).resolve().parent)) + import populace_telemetry as _telemetry + + _telemetry.push_target_diagnostics(result.diagnostics) + except Exception as _err: # noqa: BLE001 - telemetry never fails a build + log(f"telemetry skipped: {_err}") + + log( + f"done {time.time()-t0:.0f}s | loss {result.initial_loss:.3f}->" + f"{result.final_loss:.4f} | within10 " + f"{result.fraction_within_10pct*100:.2f}% | max {cw.max():,.0f} | " + f">500k {(cw>5e5).sum()}" + ) + + from policyengine_us.data import USSingleYearDataset + + shutil.copy(POOL, OUT) + ds = USSingleYearDataset(file_path=OUT) + assert len(ds.household) == n_hh + ds.household["household_weight"] = cw + for ent in ("person", "household"): + tbl = getattr(ds, ent) + if "year" in tbl.columns: + del tbl["year"] + log(f"dropped year from {ent}") + ds.save(OUT) + # No all-zero stored layers, period (the exported_nonzero gate's + # invariant): an all-zero column either masks a PE formula (a stored + # input supersedes computation) or is dead scaffolding shadowing the + # engine's own default. Drop every all-zero numeric/bool column whose + # PE default is itself zero/False — the artifact then says exactly what + # it knows and nothing else. Structural id/weight columns are kept. + from policyengine_us.system import system as _pe + + dropped_masks = [] + kept_zero = [] + for ent in ("person", "household", "tax_unit", "spm_unit", "family", + "marital_unit"): + tbl = getattr(ds, ent) + structural = {f"{ent}_id", f"{ent}_weight"} | { + c for c in tbl.columns if c.startswith("person_") + } + for c in list(tbl.columns): + if c in structural: + continue + vals = tbl[c].to_numpy() + if vals.dtype.kind not in "fiub" or np.any(vals): + continue + var = _pe.variables.get(c) + default = getattr(var, "default_value", 0) if var is not None else 0 + default_is_zero = ( + default in (0, 0.0, False) or default is None + ) + if default_is_zero: + del tbl[c] + dropped_masks.append(f"{ent}.{c}") + else: + # Dropping would CHANGE semantics (engine default != 0): + # the stored zeros are a real statement. Keep + report. + kept_zero.append(f"{ent}.{c} (default {default!r})") + if dropped_masks: + log(f"dropped {len(dropped_masks)} all-zero columns: {dropped_masks}") + ds.save(OUT) + if kept_zero: + log(f"kept all-zero columns with nonzero engine defaults: {kept_zero}") + + # other_health_insurance_premiums: usdata's decomposition — reported + # non-Medicare premiums minus the baseline-computed CHIP/marketplace/ + # Medicaid premiums, floored at zero (mirrors derive_other_health_ + # insurance_premiums; runs a baseline sim on the artifact itself). + from policyengine_us import Microsimulation as _Msim + + _sim = _Msim(dataset=USSingleYearDataset(file_path=OUT)) + # All terms mapped to person grain explicitly — the premium variables + # live at different entities (person vs tax unit). + _reported = np.asarray( + _sim.calculate( + "health_insurance_premiums_without_medicare_part_b", + 2024, + map_to="person", + ).values, + dtype=np.float64, + ) + _modeled = sum( + np.asarray( + _sim.calculate(_v, 2024, map_to="person").values, dtype=np.float64 + ) + for _v in ("chip_premium", "marketplace_net_premium", "medicaid_premium") + ) + _other = np.maximum(_reported - _modeled, 0.0) + _ds2 = USSingleYearDataset(file_path=OUT) + _ds2.person["other_health_insurance_premiums"] = _other + _ds2.save(OUT) + log( + f"other_health_insurance_premiums decomposed: nz {(_other>0).mean()*100:.1f}%" + ) + + chk = USSingleYearDataset(file_path=OUT) + assert np.array_equal( + np.asarray(chk.household["household_weight"], dtype=np.float64), cw + ) + # No misplaced-entity columns expected in v2 (fixed at export). + from collections import Counter + + tables = { + e: getattr(chk, e) + for e in ("person", "household", "tax_unit", "spm_unit", "family", "marital_unit") + } + cnt = Counter(c for t in tables.values() for c in t.columns) + dups = {c: n for c, n in cnt.items() if n > 1} + assert not dups, f"duplicate columns across entities: {dups}" + log("USSingleYearDataset verified (weights byte-identical, no dup columns)") + + shutil.copy(TP, OUT_TP) + with h5py.File(OUT_TP, "r+") as f: + assert f["household_weight/2024"].shape == cw.shape + f["household_weight/2024"][...] = cw + with h5py.File(OUT_TP) as f: + assert np.array_equal(f["household_weight/2024"][:], cw) + log("timeperiod export verified") + + np.savez_compressed( + f"{ART}/populace_us_2024_calibration.npz", + calibrated_weights=cw, + max_weight_ratio=50.0, + final_loss=result.final_loss, + within_10pct=result.fraction_within_10pct, + epochs=3000, + learning_rate=0.15, + seed=0, + ) + log("BUILD COMPLETE") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/packages/populace-data/build/us/build_populace_us_dataset.py b/packages/populace-data/build/us/build_populace_us_dataset.py index 6001fc06..58eadf6d 100644 --- a/packages/populace-data/build/us/build_populace_us_dataset.py +++ b/packages/populace-data/build/us/build_populace_us_dataset.py @@ -1,16 +1,12 @@ -"""Build the publishable populace-US dataset: calibrate the v5 pool's household -weights to the full PE-native target surface via populace.calibrate, then write -the calibrated weights back into the USSingleYearDataset H5. - -The pool (support) already beats eCPS in the symmetric-refit comparison; this -bakes in its own calibrated weights so the published dataset hits the -administrative targets out of the box. -""" +"""Build the publishable populace-US v2 dataset: calibrate the v2 pool's +household weights to its raw PE-native surface (ratio-50 bound), write them +into the USSingleYearDataset copy and the timeperiod export, verify.""" import shutil import sys import time +import h5py import numpy as np import pandas as pd @@ -18,9 +14,11 @@ from populace.frame import EntitySchema, Frame, WeightKind, Weights ART = "/Users/maxghenis/.claude-worktrees/microplex-spec-build/artifacts" -POOL = f"{ART}/spec_candidate_full_2024_v5_beats/candidate_policyengine_us.h5" -SURFACE = f"{ART}/v5_target_surface_raw.npz" -OUT = f"{ART}/populace_us_2024.h5" +POOL = f"{ART}/spec_candidate_full_2024/candidate_policyengine_us.h5" +TP = f"{ART}/spec_candidate_full_2024/candidate_timeperiod.h5" +SURFACE = f"{ART}/v2_target_surface_raw.npz" +OUT = f"{ART}/populace_us_2024_v2.h5" +OUT_TP = f"{ART}/populace_us_2024_v2_timeperiod.h5" def log(*a): @@ -28,19 +26,14 @@ def log(*a): def main(): - log("loading cached PE-native target surface...") surf = np.load(SURFACE, allow_pickle=True) - A = surf["A"].astype(np.float64) # (n_households, n_targets), RAW units - b = surf["b"].astype(np.float64) # (n_targets,), RAW units — the scaled - # surface (targets ~0.01) breaks the eCPS +1 loss and collapses weights + A = surf["A"].astype(np.float64) + b = surf["b"].astype(np.float64) w0 = surf["w0"].astype(np.float64) names = [str(x) for x in surf["names"]] - n_hh, n_targets = A.shape - log(f" {n_targets} targets over {n_hh} households; w0 sum {w0.sum()/1e6:.1f}M") + n_hh, n_t = A.shape + log(f"{n_t} targets x {n_hh} households; w0 sum {w0.sum()/1e6:.1f}M") - # Minimal Frame: one person per household, the household entity weighted by - # the pool's current weights. Calibration moves the HOUSEHOLD weights; the - # target rows are the (precompiled, PE-derived) per-household contributions. household = pd.DataFrame({"household_id": np.arange(n_hh, dtype=np.int64)}) person = pd.DataFrame( { @@ -53,10 +46,6 @@ def main(): EntitySchema(group_entities=("household",)), {"household": Weights(values=w0.copy(), kind=WeightKind.DESIGN)}, ) - - # Each PE-native target is a sum constraint whose per-household measure is the - # cached A column (aligned to the household row order, verified identical w0). - log("building target set from the PE-native surface...") targets = TargetSet( tuple( Target( @@ -66,11 +55,49 @@ def main(): value=float(b[t]), measure=(lambda _f, col=A[:, t].copy(): col), ) - for t in range(n_targets) + for t in range(n_t) + ) + ) + # Signed heavy-tail targets (architecture review): calibration was free + # to amplify loss-heavy records in dimensions absent from the surface. + # Net short-term capital gains is anchored to the PUF donor's own + # weighted, uprated total — primary-source, computed, never hand-typed. + import h5py as _h5py + + _puf = pd.read_csv( + "/Users/maxghenis/.cache/microplex/puf_2015.csv", + usecols=["P22250", "S006"], + ) + _stcg_value = float( + ( + pd.to_numeric(_puf["P22250"], errors="coerce").fillna(0) + * pd.to_numeric(_puf["S006"], errors="coerce").fillna(0) + / 100.0 + ).sum() + * 1.8 # microplex puf.py uprating factor for short_term_capital_gains + ) + with _h5py.File(TP) as _f: + _stcg_p = _f["short_term_capital_gains"]["2024"][:].astype(np.float64) + _phh = _f["person_household_id"]["2024"][:] + _hid = _f["household_id"]["2024"][:] + _hidx = {h: i for i, h in enumerate(_hid.tolist())} + _stcg_hh = np.zeros(n_hh, dtype=np.float64) + np.add.at(_stcg_hh, np.fromiter((_hidx[h] for h in _phh.tolist()), dtype=np.int64), _stcg_p) + targets = TargetSet( + tuple(targets) + + ( + Target( + name="puf/net_short_term_capital_gains", + entity="household", + aggregation="sum", + value=_stcg_value, + measure=(lambda _f2, col=_stcg_hh.copy(): col), + ), ) ) + log(f"signed STCG target appended: ${_stcg_value/1e9:.1f}B (PUF weighted, uprated)") - log("calibrating (torch APG over log-weights, full surface)...") + log("calibrating (ratio-50 bound)...") t0 = time.time() result = calibrate( frame, @@ -79,65 +106,84 @@ def main(): epochs=3000, learning_rate=0.15, mass="free", - # Hard per-record bound (the landmine guard). Unbounded calibration - # reached a 1.29M max household weight (16 records > 500k); ratio=50 - # costs ~0.1pt of within-10% and caps the max at ~318k with none - # above 500k (see bounded_recal_results.json for the sweep). max_weight_ratio=50.0, seed=0, ) - log(f" done in {time.time()-t0:.0f}s") - log(f" loss {result.initial_loss:.4f} -> {result.final_loss:.4f}") - log(f" targets within 10%: {result.fraction_within_10pct*100:.1f}%") - cw = result.frame.resolve_weights("household").values - log(f" calibrated weight sum {cw.sum()/1e6:.1f}M (was {w0.sum()/1e6:.1f}M)") - log(f" weight range: {cw.min():.2f} .. {cw.max():.0f}") - - # Write calibrated weights back into the pool H5 (via PolicyEngine's loader so - # the output stays a valid USSingleYearDataset). - log("writing calibrated dataset...") + cw = result.frame.resolve_weights("household").values.astype(np.float64) + log( + f"done {time.time()-t0:.0f}s | loss {result.initial_loss:.3f}->" + f"{result.final_loss:.4f} | within10 " + f"{result.fraction_within_10pct*100:.2f}% | max {cw.max():,.0f} | " + f">500k {(cw>5e5).sum()}" + ) + from policyengine_us.data import USSingleYearDataset shutil.copy(POOL, OUT) ds = USSingleYearDataset(file_path=OUT) - assert len(ds.household) == n_hh, "household count mismatch on write-back" + assert len(ds.household) == n_hh ds.household["household_weight"] = cw - # The pool carries a bookkeeping `year` column in several entity tables; - # it is not a PolicyEngine variable, and Microsimulation's loader refuses - # to flatten a column that appears in more than one entity. - for entity_name in ("person", "household"): - table = getattr(ds, entity_name) - if "year" in table.columns: - del table["year"] - # interest_deduction is a tax-unit variable but the pool stores it on - # persons (head carries the value, other members zero — at most one - # nonzero contributor per unit, so a group sum is exact). Microsimulation - # rejects a known variable stored at the wrong entity length. - if "interest_deduction" in ds.person.columns: - unit_sum = ( - pd.Series(ds.person["interest_deduction"].to_numpy(dtype=np.float64)) - .groupby(ds.person["person_tax_unit_id"].to_numpy()) - .sum() - ) - ds.tax_unit["interest_deduction"] = ( - unit_sum.reindex(ds.tax_unit["tax_unit_id"].to_numpy()) - .fillna(0.0) - .to_numpy(dtype=np.float64) - ) - del ds.person["interest_deduction"] + for ent in ("person", "household"): + tbl = getattr(ds, ent) + if "year" in tbl.columns: + del tbl["year"] + log(f"dropped year from {ent}") ds.save(OUT) - log(f" saved {OUT}") + # All-zero stored copies of PE FORMULA variables mask the formulas (a + # stored input supersedes computation): drop them so PE computes live — + # e.g. traditional_401k_contributions from the *_desired inputs. + from policyengine_us.system import system as _pe + + dropped_masks = [] + for ent in ("person", "household", "tax_unit", "spm_unit", "family", + "marital_unit"): + tbl = getattr(ds, ent) + for c in list(tbl.columns): + var = _pe.variables.get(c) + if var is not None and var.formulas: + vals = tbl[c].to_numpy() + if vals.dtype.kind in "fiu" and not np.any(vals): + del tbl[c] + dropped_masks.append(f"{ent}.{c}") + if dropped_masks: + log(f"dropped formula-masking zero columns: {dropped_masks}") + ds.save(OUT) + + chk = USSingleYearDataset(file_path=OUT) + assert np.array_equal( + np.asarray(chk.household["household_weight"], dtype=np.float64), cw + ) + # No misplaced-entity columns expected in v2 (fixed at export). + from collections import Counter + + tables = { + e: getattr(chk, e) + for e in ("person", "household", "tax_unit", "spm_unit", "family", "marital_unit") + } + cnt = Counter(c for t in tables.values() for c in t.columns) + dups = {c: n for c, n in cnt.items() if n > 1} + assert not dups, f"duplicate columns across entities: {dups}" + log("USSingleYearDataset verified (weights byte-identical, no dup columns)") + + shutil.copy(TP, OUT_TP) + with h5py.File(OUT_TP, "r+") as f: + assert f["household_weight/2024"].shape == cw.shape + f["household_weight/2024"][...] = cw + with h5py.File(OUT_TP) as f: + assert np.array_equal(f["household_weight/2024"][:], cw) + log("timeperiod export verified") - # persist diagnostics np.savez_compressed( - f"{ART}/populace_us_2024_calibration.npz", + f"{ART}/populace_us_2024_v2_calibration.npz", calibrated_weights=cw, - initial_loss=result.initial_loss, + max_weight_ratio=50.0, final_loss=result.final_loss, within_10pct=result.fraction_within_10pct, - loss_trajectory=result.loss_trajectory, + epochs=3000, + learning_rate=0.15, + seed=0, ) - log("CALIBRATION COMPLETE") + log("V2 BUILD COMPLETE") return 0 diff --git a/packages/populace-data/build/us/build_us_candidate.py b/packages/populace-data/build/us/build_us_candidate.py new file mode 100644 index 00000000..b9da4b46 --- /dev/null +++ b/packages/populace-data/build/us/build_us_candidate.py @@ -0,0 +1,1400 @@ +"""Build a spec-driven US eCPS-replacement candidate and optionally score it. + +Pipeline (v1 architecture, see _MISSION_JOURNAL.md): + +1. Load ASEC persons/households with raw pointer columns. +2. Construct the six-entity unit structure (microunit tax engine). +3. Aggregate persons to tax units -> the spine base frame. +4. run_spec: seeded 50/50 support spine + PUF donor imputation + (steps lifted from packs/us/specs/us-2024.yaml). +5. Assign block geography per household. +6. Re-attach persons; allocate tax-unit-imputed amounts to heads. +7. Export a USSingleYearDataset H5 gated by the eCPS export contract. +8. (--score) Run the sound eCPS-replacement comparison via the legacy + harness in ~/CosilicoAI/microplex-us. + +Smoke: .venv/bin/python scripts/build_us_candidate.py --mode smoke +Full: .venv/bin/python scripts/build_us_candidate.py --mode full --score +""" + +from __future__ import annotations + +import argparse +import json +import subprocess +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +import yaml + +REPO = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(REPO / "src")) + +from microplex.data_sources.cps import load_cps_asec # noqa: E402 +from microplex.data_sources.us_registry import ( # noqa: E402 + create_us_asec_puf_source_registry, +) +from microplex.export import ( # noqa: E402 + ExportContract, + export_policyengine_us_dataset, +) +from microplex.run import run_spec # noqa: E402 +from microplex.spec import load_spec_dict # noqa: E402 +from microplex.units import assign_us_unit_structure # noqa: E402 + +US_DATA_STORAGE = Path( + "~/PolicyEngine/policyengine-us-data/policyengine_us_data/storage" +).expanduser() +OLD_WORKTREE = Path("~/CosilicoAI/microplex-us").expanduser() +BLOCK_CROSSWALK = Path( + "~/CosilicoAI/microplex/data/block_probabilities.parquet" +).expanduser() + +MICROUNIT_RAW_COLUMNS = ( + "A_LINENO", + "A_AGE", + "A_MARITL", + "A_SPOUSE", + "PEPAR1", + "PEPAR2", + "A_EXPRRP", + "SPM_ID", + "PF_SEQ", + "A_HSCOL", +) + +# Raw ASEC columns backing CPS-threadable contract variables (mapping rules +# mirror the eCPS loader in policyengine-us-data cps.py). +EXTRA_ASEC_COLUMNS = ( + # v2 parity additions (rules mirror usdata cps.py; see V2_PLAN.md) + "ED_VAL", + "FIN_VAL", + "SRVS_VAL", + "VET_VAL", + "WC_VAL", + "OI_VAL", + "OI_OFF", + "A_HRS1", + "WKSWORK", + "POCCU2", + "PEIOOCC", + "RETCB_VAL", + "WSAL_VAL", + "SEMP_VAL", + "CSP_VAL", + "CHSP_VAL", + "DIS_VAL1", + "DIS_VAL2", + "DIS_SC1", + "DIS_SC2", + "NOW_GRP", + "NOW_MRK", + "WICYN", + "PHIP_VAL", + "POTC_VAL", + "PMED_VAL", + "PEDISDRS", + "PEDISEAR", + "PEDISEYE", + "PEDISOUT", + "PEDISPHY", + "PEDISREM", + "RESNSS1", + "RESNSS2", + "SPM_ENGVAL", + "SPM_CHILDCAREXPNS", + "NOW_OWNGRP", + "NOW_HIPAID", + "NOW_GRPFTYP", + "PERIDNUM", + "LKWEEKS", + "DST_SC1", + "DST_SC2", + "DST_SC1_YNG", + "DST_SC2_YNG", + "DST_VAL1", + "DST_VAL2", + "DST_VAL1_YNG", + "DST_VAL2_YNG", + "SPM_CAPWKCCXPNS", + # v3 parity additions (ASEC raw sources for eCPS-populated layers) + "A_FTPT", + "P_SEQ", + "NOW_NONM", + "NOW_CAID", + "NOW_OTHMT", + "NOW_MIL", + "NOW_CHAMPVA", + "NOW_VACARE", + "NOW_IHSFLG", + "SPM_CAPHOUSESUB", +) + +# Contract-required columns eCPS sources from PUF/SIPP/SCF detail our v1 +# donor surface does not carry. Explicit zero/false defaults, recorded in +# the export gate's `defaulted` accounting. Iterate in v2. +V1_ZERO_DEFAULTS: dict[str, object] = { + **{ + c: 0.0 + for c in ( + "bank_account_assets", + "bond_assets", + "stock_assets", + "partnership_se_income", + "tip_income", + "employer_sponsored_insurance_premiums", + "roth_401k_contributions", + "roth_ira_contributions", + "traditional_401k_contributions", + "self_employed_pension_contributions", + ) + }, + "receives_housing_assistance": False, + "takes_up_housing_assistance_if_eligible": False, + "takes_up_medicare_if_eligible": True, + "reported_owns_employer_sponsored_health_insurance_at_interview": False, + "is_surviving_spouse": False, +} + + +def _derive_person_columns(person: pd.DataFrame) -> pd.DataFrame: + """Derive CPS-threadable contract columns (eCPS loader rules).""" + p = person.copy() + num = lambda c: pd.to_numeric(p.get(c, 0), errors="coerce").fillna(0) # noqa: E731 + p["child_support_received"] = num("CSP_VAL").astype(float) + p["child_support_expense"] = num("CHSP_VAL").astype(float) + p["disability_benefits"] = ( + num("DIS_VAL1") * (num("DIS_SC1") != 1) + + num("DIS_VAL2") * (num("DIS_SC2") != 1) + ).astype(float) + p["has_esi"] = (num("NOW_GRP") == 1).astype(bool) + p["has_marketplace_health_coverage"] = (num("NOW_MRK") == 1).astype(bool) + p["receives_wic"] = (num("WICYN") == 1).astype(bool) + p["health_insurance_premiums_without_medicare_part_b"] = num( + "PHIP_VAL" + ).astype(float) + p["over_the_counter_health_expenses"] = num("POTC_VAL").astype(float) + p["other_medical_expenses"] = num("PMED_VAL").astype(float) + dis_flags = ["PEDISDRS", "PEDISEAR", "PEDISEYE", "PEDISOUT", "PEDISPHY", "PEDISREM"] + p["is_disabled"] = ( + pd.concat([num(c) == 1 for c in dis_flags], axis=1).any(axis=1) + ).astype(bool) + p["cps_race"] = num("race").astype(int) + p["is_female"] = (num("sex") == 2).astype(bool) + p["is_hispanic"] = (num("hispanic") == 1).astype(bool) + p["is_household_head"] = num("A_EXPRRP").isin([1, 2]).astype(bool) + p["is_separated"] = (num("A_MARITL") == 6).astype(bool) + p["is_unmarried_partner_of_household_head"] = ( + num("A_EXPRRP") == 13 + ).astype(bool) + # own children: count persons naming me as parent within the household. + key = p["household_id"].astype(str) + me = key + ":" + num("A_LINENO").astype(int).astype(str) + par1 = key + ":" + num("PEPAR1").astype(int).astype(str) + par2 = key + ":" + num("PEPAR2").astype(int).astype(str) + counts = par1.value_counts().add(par2.value_counts(), fill_value=0) + p["own_children_in_household"] = me.map(counts).fillna(0).astype(float) + # household child counts (eCPS groups by household). + u18 = (num("A_AGE") < 18).groupby(p["household_id"]).transform("sum") + u6 = (num("A_AGE") < 6).groupby(p["household_id"]).transform("sum") + p["count_under_18"] = u18.astype(float) + p["count_under_6"] = u6.astype(float) + # Social security split by RESNSS reason codes, age-62 fallback. + ss = pd.to_numeric(p.get("social_security", 0), errors="coerce").fillna(0) + r1, r2 = num("RESNSS1"), num("RESNSS2") + retired = (r1 == 1) | (r2 == 1) + disabled = ((r1 == 2) | (r2 == 2)) & ~retired + survivor = (r1.isin([3, 5]) | r2.isin([3, 5])) & ~retired & ~disabled + dependent = ( + (r1.isin([4, 6, 7]) | r2.isin([4, 6, 7])) + & ~retired + & ~disabled + & ~survivor + ) + unclassified = (ss > 0) & ~(retired | disabled | survivor | dependent) + age = num("A_AGE") + retired = retired | (unclassified & (age >= 62)) + disabled = disabled | (unclassified & (age < 62)) + p["social_security_retirement"] = (ss * retired).astype(float) + p["social_security_disability"] = (ss * disabled).astype(float) + p["social_security_survivors"] = (ss * survivor).astype(float) + p["social_security_dependents"] = (ss * dependent).astype(float) + # ---- v2 parity derivations (rules mirror usdata cps.py; V2_PLAN.md) ---- + p["educational_assistance"] = num("ED_VAL").astype(float) + p["financial_assistance"] = num("FIN_VAL").astype(float) + p["survivor_benefits"] = num("SRVS_VAL").astype(float) + p["veterans_benefits"] = num("VET_VAL").astype(float) + p["workers_compensation"] = num("WC_VAL").astype(float) + oi_val, oi_off = num("OI_VAL"), num("OI_OFF") + strike = oi_off == 12 + alimony_oi = oi_off == 20 + p["strike_benefits"] = (oi_val * strike).astype(float) + p["miscellaneous_income"] = ( + oi_val * ~(strike | alimony_oi) + ).astype(float) + p["hours_worked_last_week"] = num("A_HRS1").clip(lower=0).astype(float) + p["weeks_worked"] = num("WKSWORK").clip(0, 52).astype(float) + p["detailed_occupation_recode"] = num("POCCU2").astype(float) + # Weeks looking for work (LKWEEKS: -1 NIU -> 0), mirroring usdata. + p["weeks_unemployed"] = num("LKWEEKS").clip(lower=0).astype(float) + # v3 parity: ASEC raw derivations mirroring usdata cps.py exactly. + p["is_blind"] = (num("PEDISEYE") == 1).astype(bool) + p["is_surviving_spouse"] = (num("A_MARITL") == 4).astype(bool) + p["is_full_time_college_student"] = ( + (num("A_HSCOL") == 2) & (num("A_FTPT") == 1) + ).astype(bool) + # Household head: person sequence 1 within the household (usdata P_SEQ rule). + p["is_household_head"] = (num("P_SEQ") == 1).astype(bool) + # Current health coverage at interview (ASEC NOW_* flags; 1 = covered). + for _pe_name, _now in { + "has_marketplace_health_coverage_at_interview": "NOW_MRK", + "has_non_marketplace_direct_purchase_health_coverage_at_interview": "NOW_NONM", + "has_medicaid_health_coverage_at_interview": "NOW_CAID", + "has_other_means_tested_health_coverage_at_interview": "NOW_OTHMT", + "has_tricare_health_coverage_at_interview": "NOW_MIL", + "has_champva_health_coverage_at_interview": "NOW_CHAMPVA", + "has_va_health_coverage_at_interview": "NOW_VACARE", + "has_indian_health_service_coverage_at_interview": "NOW_IHSFLG", + }.items(): + p[_pe_name] = (num(_now) == 1).astype(bool) + # SPM-reported housing assistance (capped housing subsidy > 0). + p["receives_housing_assistance"] = (num("SPM_CAPHOUSESUB") > 0).astype(bool) + # FLSA overtime occupation flags from POCCU2 (codes shipped by the engine). + from policyengine_us.data.cps import ( + CPS_FLSA_EXECUTIVE_ADMINISTRATIVE_PROFESSIONAL_OCCUPATION_CODES as _EXEC, + CPS_FLSA_OVERTIME_OCCUPATION_CODES as _OCC, + ) + + for _flag, _code in _OCC.items(): + p[_flag] = (num("POCCU2") == _code).astype(bool) + p["is_executive_administrative_professional"] = ( + num("POCCU2").isin(list(_EXEC)).astype(bool) + ) + # RETCB proportional split (usdata cps.py:1505-1552; shares from + # imputation_parameters.yaml — BEA/FRED + IRS SOI administrative shares). + retcb = num("RETCB_VAL").clip(lower=0) + has_wages = num("WSAL_VAL") > 0 + has_se = num("SEMP_VAL") > 0 + has_earned = has_wages | has_se + se_pension = retcb * 0.046 * has_se + p["self_employed_pension_contributions_desired"] = se_pension.astype(float) + remaining = (retcb - se_pension).clip(lower=0) + dc_pool = remaining * 0.908 * has_wages + ira_pool = (remaining - dc_pool) * has_earned + p["traditional_401k_contributions_desired"] = (dc_pool * 0.85).astype(float) + p["roth_401k_contributions_desired"] = (dc_pool * 0.15).astype(float) + p["traditional_ira_contributions_desired"] = (ira_pool * 0.392).astype(float) + p["roth_ira_contributions_desired"] = (ira_pool * 0.608).astype(float) + return p + +# Person-level ASEC harmonized income columns that sum to tax-unit totals. +PERSON_INCOME_COLUMNS = ( + "employment_income", + "self_employment_income", + "taxable_interest_income", + "rental_income", + "social_security", + "taxable_pension_income", + "unemployment_compensation", +) + + +# Donor-named variables the PUF source actually carries, imputed at +# tax-unit grain. CPS-measured ones (also on the spine base) are listed in +# CPS_MEASURED; the rest are PUF-only detail imputed onto both halves. +PUF_IMPUTE_VARS = ( + "employment_income", + "self_employment_income", + "social_security", + "taxable_pension_income", + "taxable_interest_income", + "unemployment_compensation", + "rental_income", + "partnership_s_corp_income", + "farm_income", + "tax_exempt_interest_income", + "qualified_dividend_income", + "ordinary_dividend_income", + "short_term_capital_gains", + "long_term_capital_gains", + "taxable_pension_income", + "total_pension_income", + "ira_distributions", + "alimony_received", + "charitable_cash", + "charitable_noncash", + "mortgage_interest_paid", + "real_estate_tax_paid", + "student_loan_interest", + "ira_deduction", + "farm_income", + # v2 parity fields (donor names from microplex puf.py field map) + "alimony_expense", + "casualty_loss", + "domestic_production_ald", + "educator_expense", + "estate_income", + "health_savings_account_ald", + "long_term_capital_gains_on_collectibles", + "unreimbursed_business_employee_expenses", + "qualified_tuition_expenses", + "business_is_sstb", + "sstb_self_employment_income_would_be_qualified", + "farm_rent_income", + "self_employed_pension_contribution_ald", + "unrecaptured_section_1250_gain", + "puf_miscellaneous_income", + "salt_refund_income", + "investment_income_elected_form_4952", + "capital_gains_distributions", + # v3 QBI/partnership block (derived in the PUF loader from usdata rules) + "partnership_se_income", + "w2_wages_from_qualified_business", + "unadjusted_basis_qualified_property", + "sstb_self_employment_income", + "sstb_w2_wages_from_qualified_business", + "sstb_unadjusted_basis_qualified_property", + "qualified_bdc_income", + "qualified_reit_and_ptp_income", +) +CPS_MEASURED = ( + "employment_income", + "self_employment_income", + "social_security", + "taxable_pension_income", + "taxable_interest_income", + "unemployment_compensation", + "rental_income", +) + +# donor/common name -> PolicyEngine contract name at person allocation. +DONOR_TO_PE = { + "ira_distributions": "taxable_ira_distributions", + "alimony_received": "alimony_income", + "charitable_cash": "charitable_cash_donations", + "charitable_noncash": "charitable_non_cash_donations", + "mortgage_interest_paid": "interest_deduction", + "real_estate_tax_paid": "real_estate_taxes", + "ordinary_dividend_income": "non_qualified_dividend_income", + "ira_deduction": "traditional_ira_contributions", + "puf_miscellaneous_income": "miscellaneous_income", + "capital_gains_distributions": "non_sch_d_capital_gains", + "sstb_self_employment_income": "sstb_self_employment_income_before_lsr", +} + +SHARED_PREDICTORS = ("age", "is_joint", "n_people", "n_children") + + +def _build_imputation_steps(*, weighted: bool = True) -> list[dict]: + """ASEC+PUF imputation graph over donor-available variables. + + weighted=True fits donors with the PUF design weight — the verified + landmine fix (microplex#76): the PUF oversamples high incomes, so + unweighted fits draw from the sample distribution (measured here: + full-loss 130.1 vs eCPS 1.41, candidate LTCG $200T vs $257B). + """ + puf_vars = list(dict.fromkeys(PUF_IMPUTE_VARS)) + puf_only = [v for v in puf_vars if v not in CPS_MEASURED] + w = {"weights": "weight"} if weighted else {} + return [ + {"onto": "synthetic_puf", "from": "puf", "vars": puf_vars, + "order": "spine_first", **w}, + {"onto": "cps_keep", "from": "puf", "vars": puf_only, + "condition_on": ["demographics", *CPS_MEASURED], **w}, + ] + + +def _aggregate_tax_units(person: pd.DataFrame, tax_unit: pd.DataFrame) -> pd.DataFrame: + """Aggregate persons to the tax-unit-grain spine base frame.""" + g = person.groupby("person_tax_unit_id", sort=True) + base = pd.DataFrame(index=g.size().index) + base["tax_unit_id"] = base.index + base["household_id"] = g["household_id"].first() + base["n_people"] = g.size().astype(float) + is_head = person["tax_unit_role_input"] == "HEAD" + head_age = ( + person.loc[is_head] + .groupby("person_tax_unit_id")["A_AGE"] + .max() + .astype(float) + ) + base["age"] = head_age.reindex(base.index).fillna( + g["A_AGE"].max().astype(float) + ) + base["n_children"] = ( + person.assign(_child=(person["A_AGE"] < 18).astype(float)) + .groupby("person_tax_unit_id")["_child"] + .sum() + .reindex(base.index) + .fillna(0.0) + ) + base["is_joint"] = ( + person.assign(_sp=(person["tax_unit_role_input"] == "SPOUSE").astype(float)) + .groupby("person_tax_unit_id")["_sp"] + .max() + .reindex(base.index) + .fillna(0.0) + ) + for col in PERSON_INCOME_COLUMNS: + if col in person.columns: + base[col] = g[col].sum() + fs = tax_unit.set_index("tax_unit_id")["filing_status_input"] + base["filing_status_input"] = base["tax_unit_id"].map(fs) + base = base.reset_index(drop=True) + return base + + +def _attach_household_columns( + base: pd.DataFrame, households: pd.DataFrame +) -> pd.DataFrame: + keep = ["household_id", "state_fips", "household_weight", "tenure_type"] + have = [c for c in keep if c in households.columns] + return base.merge( + households[have].drop_duplicates("household_id"), + on="household_id", + how="left", + ) + + +def _assign_blocks( + households: pd.DataFrame, crosswalk_path: Path, seed: int +) -> pd.DataFrame: + """Probability-weighted census block assignment per household by state.""" + xw = pd.read_parquet(crosswalk_path) + rng = np.random.default_rng(seed) + households = households.copy() + xw["state_fips"] = xw["state_fips"].astype(int) + block = pd.Series("", index=households.index, dtype=object) + county = pd.Series("", index=households.index, dtype=object) + cd = pd.Series(0, index=households.index, dtype="int64") + for state, idx in households.groupby( + households["state_fips"].astype(int) + ).groups.items(): + pool = xw[xw["state_fips"] == state] + if pool.empty: + pool = xw + p = pool["prob"].to_numpy() + p = p / p.sum() + draw = rng.choice(len(pool), size=len(idx), p=p) + chosen = pool.iloc[draw] + geoid = chosen["geoid"].astype(str).to_numpy() + block.loc[idx] = geoid + county.loc[idx] = [g[:5] for g in geoid] + district = ( + chosen["cd_id"] + .astype(str) + .str.extract(r"(\d+)$")[0] + .fillna("0") + .astype(int) + .to_numpy() + ) + cd.loc[idx] = int(state) * 100 + district + households["block_geoid"] = block.astype(str) + households["county_fips"] = county.astype(str) + households["congressional_district_geoid"] = cd.astype(np.int32) + return households + + +def _build_unit_map( + base: pd.DataFrame, halves: dict[str, pd.DataFrame] +) -> pd.DataFrame: + """Original->engine unit-id mapping recovered via the halves' row index. + + The support spine re-identifies the synthetic half's tax_unit_id (they + are new synthetic units) and zeroes its weights; the per-half frames + keep the base frame's positional index, which is the join key back to + the original ids. + """ + min_id = int(base["tax_unit_id"].min()) + offset = int(base["tax_unit_id"].max()) - min(0, min_id) + 1 + pieces = [] + for half_name, frame in halves.items(): + new_ids = frame["tax_unit_id"].to_numpy() + orig_ids = new_ids - offset if half_name == "synthetic_puf" else new_ids + pieces.append( + pd.DataFrame( + { + "orig_tax_unit_id": orig_ids, + "new_tax_unit_id": new_ids, + "_half": half_name, + } + ) + ) + unit_map = pd.concat(pieces, ignore_index=True) + base_ids = set(base["tax_unit_id"].tolist()) + recovered = set(unit_map["orig_tax_unit_id"].tolist()) + if unit_map["orig_tax_unit_id"].duplicated().any(): + raise ValueError("unit map has duplicated original tax unit ids") + if recovered != base_ids: + raise ValueError( + "recovered original ids do not partition the base " + f"(missing {len(base_ids - recovered)}, " + f"extra {len(recovered - base_ids)})" + ) + return unit_map + + +def _allocate_to_persons( + person: pd.DataFrame, + spine: pd.DataFrame, + imputed_vars: list[str], + unit_map: pd.DataFrame, +) -> pd.DataFrame: + """Re-attach spine tax-unit values to persons via the unit map. + + cps_keep persons keep their ASEC person-level values for CPS-measured + columns; PUF-only imputed amounts go to the unit head. synthetic_puf + persons get every imputed variable head-allocated (others zero). + """ + person = person.copy() + m = unit_map.set_index("orig_tax_unit_id") + person["_half"] = person["person_tax_unit_id"].map(m["_half"]) + person["new_tax_unit_id"] = person["person_tax_unit_id"].map( + m["new_tax_unit_id"] + ) + unmapped = int(person["_half"].isna().sum()) + if unmapped: + raise ValueError(f"{unmapped} persons failed to map to a spine half") + spine_idx = spine.set_index("tax_unit_id") + is_head = person["tax_unit_role_input"] == "HEAD" + synthetic = person["_half"] == "synthetic_puf" + + for var in imputed_vars: + if var not in spine_idx.columns: + continue + unit_value = person["new_tax_unit_id"].map(spine_idx[var]) + head_alloc = np.where(is_head, unit_value.fillna(0.0), 0.0) + if var in person.columns: + # CPS-measured: keep person values on cps_keep, head-allocate + # the synthetic draw on the synthetic half. + person[var] = np.where( + synthetic, head_alloc, person[var].fillna(0.0) + ) + else: + person[var] = head_alloc + return person + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--mode", choices=["smoke", "full"], default="smoke") + ap.add_argument("--asec-year", type=int, default=2025) + ap.add_argument("--calendar-year", type=int, default=2024) + ap.add_argument("--puf-year", type=int, default=2024) + ap.add_argument("--max-tax-units", type=int, default=None) + ap.add_argument("--max-puf-rows", type=int, default=None) + ap.add_argument("--seed", type=int, default=20260529) + ap.add_argument("--output-dir", type=Path, default=None) + ap.add_argument("--score", action="store_true") + ap.add_argument( + "--baseline-h5", + type=Path, + default=OLD_WORKTREE / "artifacts/baselines/enhanced_cps_2024_hf_main.h5", + ) + ap.add_argument( + "--usdata-repo", + type=Path, + default=Path("~/.claude-worktrees/usdata-f7458313").expanduser(), + ) + ap.add_argument( + "--tail-units", + type=int, + default=15000, + help="Top-income PUF returns carried verbatim as a tail-support " + "stratum at design weights (0 disables).", + ) + args = ap.parse_args() + + # The usdata package provides primary-source rules (tipped occupations, + # QBI simulators, SCF/ACS loaders); make it importable for every stage. + sys.path.insert(0, str(args.usdata_repo)) + + smoke = args.mode == "smoke" + max_units = args.max_tax_units or (4000 if smoke else None) + max_puf = args.max_puf_rows or (8000 if smoke else None) + out = args.output_dir or ( + REPO / "artifacts" / f"spec_candidate_{args.mode}_{args.calendar_year}" + ) + out.mkdir(parents=True, exist_ok=True) + log = lambda *a: print("[build]", *a, flush=True) # noqa: E731 + + # ---- Stage A: ASEC persons + households ------------------------------- + log("stage A: loading ASEC persons/households") + ds = load_cps_asec( + year=args.asec_year, + extra_person_columns=( + list(MICROUNIT_RAW_COLUMNS) + ["PH_SEQ"] + list(EXTRA_ASEC_COLUMNS) + ), + extra_household_columns=["H_TENURE", "GTCO"], + ) + persons = ds.persons.to_pandas() + households = ds.households.to_pandas() + if "GTCO" in households.columns: + # NYC boroughs by state*1000+county FIPS (usdata cps.py NYC list). + _scf = ( + pd.to_numeric(households.get("state_fips", 0), errors="coerce").fillna(0) + * 1000 + + pd.to_numeric(households["GTCO"], errors="coerce").fillna(0) + ) + households["in_nyc"] = _scf.isin([36005, 36047, 36061, 36081, 36085]).astype(bool) + if "H_TENURE" in households.columns: + # Same base map as the eCPS loader (usdata cps.py:418). + households["tenure_type"] = ( + pd.to_numeric(households["H_TENURE"], errors="coerce") + .fillna(0) + .map({0: "NONE", 1: "OWNED_WITH_MORTGAGE", 2: "RENTED", 3: "NONE"}) + ) + log(f" persons={len(persons):,} households={len(households):,}") + + # ---- Stage B: unit structure ------------------------------------------ + log("stage B: unit assignment (microunit)") + units = assign_us_unit_structure(persons, year=args.calendar_year) + person = units.person.rename( + columns={ + "wage_income": "employment_income", + "interest_income": "taxable_interest_income", + } + ) + log( + f" tax_units={len(units.tax_unit):,} spm={len(units.spm_unit):,} " + f"families={len(units.family):,} marital={len(units.marital_unit):,}" + ) + + # ---- Stage C: tax-unit spine base -------------------------------------- + log("stage C: aggregating tax-unit spine base") + units_tu = units.tax_unit.rename(columns={"TAX_ID": "tax_unit_id"}) + base = _aggregate_tax_units(person, units_tu) + base = _attach_household_columns(base, households) + if max_units is not None: + keep_units = base["tax_unit_id"].head(max_units) + base = base[base["tax_unit_id"].isin(keep_units)].copy() + person = person[person["person_tax_unit_id"].isin(keep_units)].copy() + log(f" spine base tax units: {len(base):,}") + + # ---- Stage D: run_spec -------------------------------------------------- + log("stage D: run_spec (support spine + PUF imputation)") + steps = _build_imputation_steps() + imputed_vars = sorted({v for s in steps for v in s.get("vars", [])}) + registry = create_us_asec_puf_source_registry( + asec_year=args.asec_year, + calendar_year=args.calendar_year, + puf_year=args.puf_year, + puf_path=US_DATA_STORAGE / "puf_2015.csv", + puf_demographics_path=US_DATA_STORAGE / "demographics_2015.csv", + ) + id_keep = [ + "tax_unit_id", + "household_id", + "household_weight", + "state_fips", + "filing_status_input", + ] + spec = load_spec_dict( + { + "meta": {"country": "us", "model_year": args.calendar_year}, + "sources": { + "cps_asec": { + "dataset": ( + f"cps_asec_{args.asec_year}_calendar_{args.calendar_year}" + ), + "role": "spine", + }, + "puf": {"dataset": f"puf_{args.puf_year}", "role": "donor"}, + }, + "spine": { + "base": "cps_asec", + "method": "support_spine", + "support": {"seed": args.seed}, + "halves": [ + {"name": "cps_keep", "keep": "all"}, + { + "name": "synthetic_puf", + "strip_to": ["demographics", *id_keep], + }, + ], + }, + "imputation": steps, + } + ) + puf = registry.resolve_source(spec, "puf") + if max_puf is not None: + puf = puf.head(max_puf).copy() + # Tail-support stratum: the top returns by income proxy, kept verbatim at + # design weights BEFORE the population resample discards them. Free + # support for top-bracket targets; calibration owns how much to use. + income_proxy = sum( + puf[c].clip(lower=0) + for c in ( + "employment_income", + "self_employment_income", + "partnership_s_corp_income", + "taxable_interest_income", + "qualified_dividend_income", + "ordinary_dividend_income", + "long_term_capital_gains", + "short_term_capital_gains", + "taxable_pension_income", + ) + if c in puf.columns + ) + tail = ( + puf.loc[income_proxy.nlargest(min(args.tail_units, len(puf) // 20)).index].copy() + if args.tail_units + else puf.head(0).copy() + ) + log(f" tail stratum: {len(tail):,} units, weight {tail['weight'].sum()/1e6 if len(tail) else 0:.2f}M") + + # Population-resample the donor by its design weight: the PUF oversamples + # high incomes, and microimpute's canonical Imputer silently ignores + # weight_col (verified 2026-06-10), so importance-resampling is the + # correct way to get population conditionals from the fit. + rng = np.random.default_rng(args.seed) + pw = puf["weight"].to_numpy(dtype=float) + pw = pw / pw.sum() + puf = puf.iloc[ + rng.choice(len(puf), size=len(puf), replace=True, p=pw) + ].reset_index(drop=True) + # Harmonize donor names + derive the shared predictor surface. + puf = puf.rename(columns={"gross_social_security": "social_security"}) + puf["is_joint"] = (puf["filing_status"] == "JOINT").astype(float) + puf["n_people"] = puf["exemptions_count"].clip(lower=1).astype(float) + puf["n_children"] = puf["ctc_children"].fillna(0).astype(float) + puf["age"] = puf["age"].astype(float) + log(f" puf donor rows: {len(puf):,}") + result = run_spec( + spec, + {"cps_asec": base, "puf": puf}, + demographic_columns=SHARED_PREDICTORS, + spine_keywords=( + "employment_income", + "self_employment_income", + "social_security", + "taxable_pension_income", + "taxable_interest_income", + ), + ) + spine = result.frame + half_col = [c for c in spine.columns if c.startswith("_spine")][0] + log( + f" spine rows={len(spine):,} halves=" + f"{spine[half_col].value_counts().to_dict()}" + ) + + # ---- Stage E: geography ------------------------------------------------- + log("stage E: block geography") + households = _assign_blocks( + households[households["household_id"].isin(person["household_id"])], + BLOCK_CROSSWALK, + args.seed, + ) + + # ---- Stage F: entity assembly ------------------------------------------ + log("stage F: person re-attach + entity tables") + unit_map = _build_unit_map(base, dict(result.halves)) + orig_household_size = person.groupby("household_id").size() + person = _allocate_to_persons(person, spine, imputed_vars, unit_map) + # Donor/common names -> PolicyEngine contract names; derived splits. + person = person.rename(columns=DONOR_TO_PE) + if {"total_pension_income", "taxable_pension_income"} <= set(person.columns): + person["tax_exempt_pension_income"] = ( + person["total_pension_income"] - person["taxable_pension_income"] + ).clip(lower=0.0) + person = person.drop(columns=["total_pension_income"]) + if { + "non_qualified_dividend_income", + "qualified_dividend_income", + } <= set(person.columns): + person["non_qualified_dividend_income"] = ( + person["non_qualified_dividend_income"] + - person["qualified_dividend_income"] + ).clip(lower=0.0) + person = _derive_person_columns(person) + # Tipped-occupation codes via the eCPS's own mapping (usdata module). + from policyengine_us_data.datasets.cps.tipped_occupation import ( + derive_is_tipped_occupation, + derive_treasury_tipped_occupation_code, + ) + + ttoc = derive_treasury_tipped_occupation_code(person["PEIOOCC"]) + person["treasury_tipped_occupation_code"] = ttoc.astype(float) + person["is_tipped_occupation"] = derive_is_tipped_occupation(ttoc) + + # Retirement-account distributions: ASEC DST_SC/DST_VAL pairs by source + # code, split by the usdata taxable fractions (imputation_parameters.yaml). + import yaml as _yaml + + _ipar = _yaml.safe_load( + (args.usdata_repo / "policyengine_us_data" / "datasets" / "cps" + / "imputation_parameters.yaml").read_text() + ) + _codes = {1: "401k", 2: "403b", 6: "sep"} + # Codes without a taxable split: keogh (5) is its own PE input; roth IRA + # (3) is tax-exempt by assumption (usdata cps.py RETIREMENT_CODES). + for _code, _pe in ((5, "keogh_distributions"), (3, "tax_exempt_ira_distributions")): + _tot = 0 + for _i in ("1", "2", "1_YNG", "2_YNG"): + _sc = pd.to_numeric(person.get(f"DST_SC{_i}", 0), errors="coerce").fillna(0) + _val = pd.to_numeric(person.get(f"DST_VAL{_i}", 0), errors="coerce").fillna(0) + _tot = _tot + (_sc == _code) * _val + person[_pe] = pd.Series(_tot, index=person.index).astype(float) + for _code, _name in _codes.items(): + _tot = 0 + for _i in ("1", "2", "1_YNG", "2_YNG"): + _sc = pd.to_numeric(person.get(f"DST_SC{_i}", 0), errors="coerce").fillna(0) + _val = pd.to_numeric(person.get(f"DST_VAL{_i}", 0), errors="coerce").fillna(0) + _tot = _tot + (_sc == _code) * _val + _frac = float(_ipar[f"taxable_{_name}_distribution_fraction"]) + person[f"taxable_{_name}_distributions"] = (_tot * _frac).astype(float) + person[f"tax_exempt_{_name}_distributions"] = (_tot * (1 - _frac)).astype(float) + log(" retirement distributions split (DST codes x usdata taxable fractions)") + # Pre-response copies and aliases the contract requires alongside the + # base variables. + person["employment_income_before_lsr"] = person["employment_income"] + person["self_employment_income_before_lsr"] = person[ + "self_employment_income" + ] + person["long_term_capital_gains_before_response"] = person[ + "long_term_capital_gains" + ] + person["taxable_unemployment_compensation"] = person[ + "unemployment_compensation" + ] + person["farm_operations_income"] = person["farm_income"] + person["taxable_private_pension_income"] = person["taxable_pension_income"] + person["tax_exempt_private_pension_income"] = person[ + "tax_exempt_pension_income" + ] + # Re-key every unit system per (original id, half): the synthetic half's + # units are new synthetic entities, and a multi-unit household whose + # units land in different halves splits into per-half export households. + person["person_tax_unit_id"] = person["new_tax_unit_id"].astype(np.int64) + for unit in ("spm_unit", "family", "marital_unit"): + key = ( + person[f"person_{unit}_id"].astype(str) + "|" + person["_half"] + ) + person[f"person_{unit}_id"] = (pd.factorize(key)[0] + 1).astype( + np.int64 + ) + person["_orig_household_id"] = person["household_id"] + hh_key = person["household_id"].astype(str) + "|" + person["_half"] + person["person_household_id"] = (pd.factorize(hh_key)[0] + 1).astype( + np.int64 + ) + person["person_id"] = np.arange(1, len(person) + 1, dtype=np.int64) + person["age"] = person["A_AGE"].astype(float) + # eCPS source flags: synthetic_puf half maps to the PUF-clone marker the + # loss surface's nation/source/* household-count targets read. + person["person_is_puf_clone"] = (person["_half"] == "synthetic_puf").astype( + bool + ) + + # ---- Tail-support stratum persons ------------------------------------ + if len(tail): + t = tail.reset_index(drop=True).copy() + t = t.rename(columns={"gross_social_security": "social_security"}) + t = t.rename(columns=DONOR_TO_PE) + if {"total_pension_income", "taxable_pension_income"} <= set(t.columns): + t["tax_exempt_pension_income"] = ( + t["total_pension_income"] - t["taxable_pension_income"] + ).clip(lower=0.0) + if { + "non_qualified_dividend_income", + "qualified_dividend_income", + } <= set(t.columns): + t["non_qualified_dividend_income"] = ( + t["non_qualified_dividend_income"] + - t["qualified_dividend_income"] + ).clip(lower=0.0) + n = len(t) + joint = (t["filing_status"] == "JOINT").to_numpy() + base_ids = { + "tax": int(person["person_tax_unit_id"].max()), + "hh": int(person["person_household_id"].max()), + "spm": int(person["person_spm_unit_id"].max()), + "fam": int(person["person_family_id"].max()), + "mar": int(person["person_marital_unit_id"].max()), + "pid": int(person["person_id"].max()), + } + head = pd.DataFrame(index=range(n)) + head["person_tax_unit_id"] = base_ids["tax"] + 1 + np.arange(n) + head["person_household_id"] = base_ids["hh"] + 1 + np.arange(n) + head["person_spm_unit_id"] = base_ids["spm"] + 1 + np.arange(n) + head["person_family_id"] = base_ids["fam"] + 1 + np.arange(n) + head["person_marital_unit_id"] = base_ids["mar"] + 1 + np.arange(n) + head["age"] = ( + pd.to_numeric(t["age"], errors="coerce").fillna(50).clip(18, 85) + ) + head["A_AGE"] = head["age"] + head["tax_unit_role_input"] = "HEAD" + head["_half"] = "tail_puf" + value_cols = [ + c + for c in t.columns + if c in set(person.columns) + and c + not in ("age", "weight", "tax_unit_id", "filing_status", "_survey") + and pd.api.types.is_numeric_dtype(t[c]) + ] + for c in value_cols: + head[c] = pd.to_numeric(t[c], errors="coerce").fillna(0.0) + head["own_children_in_household"] = ( + pd.to_numeric(t.get("ctc_children", 0), errors="coerce") + .fillna(0) + .to_numpy() + ) + head["count_under_18"] = head["own_children_in_household"] + spouse = head[joint].copy() + for c in value_cols: + spouse[c] = 0.0 + spouse["tax_unit_role_input"] = "SPOUSE" + spouse["own_children_in_household"] = 0.0 + tail_person = pd.concat([head, spouse], ignore_index=True) + tail_person["person_id"] = ( + base_ids["pid"] + 1 + np.arange(len(tail_person)) + ) + tail_person["person_is_puf_clone"] = True + for col in ( + "employment_income_before_lsr", + "self_employment_income_before_lsr", + ): + src = col.replace("_before_lsr", "") + if src in tail_person.columns: + tail_person[col] = tail_person[src] + if "long_term_capital_gains" in tail_person.columns: + tail_person["long_term_capital_gains_before_response"] = ( + tail_person["long_term_capital_gains"] + ) + if "unemployment_compensation" in tail_person.columns: + tail_person["taxable_unemployment_compensation"] = tail_person[ + "unemployment_compensation" + ] + if "taxable_pension_income" in tail_person.columns: + tail_person["taxable_private_pension_income"] = tail_person[ + "taxable_pension_income" + ] + if "tax_exempt_pension_income" in tail_person.columns: + tail_person["tax_exempt_private_pension_income"] = tail_person[ + "tax_exempt_pension_income" + ] + if "farm_income" in tail_person.columns: + tail_person["farm_operations_income"] = tail_person["farm_income"] + tail_person = tail_person.reindex(columns=person.columns) + for c in person.columns: + if tail_person[c].isna().all(): + dt = person[c].dtype + if pd.api.types.is_bool_dtype(dt): + tail_person[c] = False + elif pd.api.types.is_numeric_dtype(dt): + tail_person[c] = 0 + else: + tail_person[c] = "" + tail_person["_orig_household_id"] = np.nan + log( + f" tail persons: {len(tail_person):,} " + f"({int(joint.sum()):,} spouses added)" + ) + else: + tail_person = person.head(0) + + + # Export households = (original household, half) pieces. Attributes come + # from the original household; weight is prorated by the piece's person + # share so person-mass totals are preserved exactly (full ASEC scale). + piece = ( + person.groupby("person_household_id") + .agg( + _orig=("_orig_household_id", "first"), + _half=("_half", "first"), + _n=("person_id", "size"), + ) + .reset_index() + .rename(columns={"person_household_id": "household_id"}) + ) + hh_attrs = households.drop_duplicates("household_id").rename( + columns={"household_id": "_orig", "household_weight": "_orig_weight"} + ) + hh = piece.merge(hh_attrs, on="_orig", how="left") + hh["household_weight"] = ( + hh["_orig_weight"] + * hh["_n"] + / hh["_orig"].map(orig_household_size).to_numpy() + ).astype(float) + hh["tract_geoid"] = hh["block_geoid"].astype(str).str[:11] + hh["household_is_puf_clone"] = (hh["_half"] == "synthetic_puf").astype(bool) + hh = hh.drop(columns=["_orig", "_half", "_n", "_orig_weight"]) + + # Tail-stratum households: design weights, geography sampled from the + # main household distribution (the PUF has no state). + if len(tail): + rng_t = np.random.default_rng(args.seed + 1) + donor_geo = hh.sample( + n=len(tail), + replace=True, + weights=hh["household_weight"].clip(lower=1e-9), + random_state=int(rng_t.integers(0, 2**31)), + ).reset_index(drop=True) + tail_hh = donor_geo.copy() + tail_hh["household_id"] = ( + tail_person["person_household_id"].unique() + ) + tail_hh["household_weight"] = tail["weight"].to_numpy(dtype=float) + tail_hh["household_is_puf_clone"] = True + hh = pd.concat([hh, tail_hh], ignore_index=True) + person = pd.concat([person, tail_person], ignore_index=True) + log(f" households incl. tail: {len(hh):,}") + + # ---- Stage F2: primary-source imputation (v3, eCPS-free) --------------- + # Wealth from Fed SCF, tips from SIPP, wages from CPS-ORG. The enhanced + # CPS is never a build input; it remains only the scoring benchmark. + log("stage F2: primary-source imputation (SCF / SIPP / ORG)") + sys.path.insert(0, str(args.usdata_repo)) + import primary_source_impute as psi + + person, hh = psi.add_scf_wealth(person, hh, seed=args.seed, log=log) + person = psi.add_sipp_tips(person, log=log) + person = psi.add_org_wages(person, hh, args.calendar_year, log=log) + person = psi.add_meps_esi_premiums(person, log=log) + person = psi.add_prior_year_income(person, args.asec_year, log=log) + person = psi.add_mortgage_conversion(person, hh, args.calendar_year, log=log) + person, hh = psi.add_acs_rent(person, hh, seed=args.seed, log=log) + person, hh = psi.add_vehicle_assets(person, hh, log=log) + + def _group_clone_flag(id_col: str) -> pd.Series: + share = person.groupby(person[id_col])["person_is_puf_clone"].mean() + return share > 0.5 + + clone_flags = { + c: _group_clone_flag(f"person_{c.split('_is_')[0]}_id") + for c in ( + "household_is_puf_clone", + "tax_unit_is_puf_clone", + "spm_unit_is_puf_clone", + "family_is_puf_clone", + ) + } + + def unit_table(id_col: str, source: pd.DataFrame | None = None) -> pd.DataFrame: + ids = np.sort(person[f"person_{id_col}"].unique()) + t = pd.DataFrame({id_col: ids}) + if source is not None: + extra = source.rename(columns={"TAX_ID": id_col}) + t = t.merge(extra, on=id_col, how="left") + flag = f"{id_col.rsplit('_id', 1)[0]}_is_puf_clone" + if flag in clone_flags: + t[flag] = t[id_col].map(clone_flags[flag]).fillna(False).astype(bool) + return t + + spm = unit_table("spm_unit_id") + # SPM-record childcare expenses: raw ASEC columns, max per SPM unit + # (constant within unit on the SPM record), mirroring the energy subsidy. + for raw, pe_name in ( + ("SPM_CHILDCAREXPNS", "spm_unit_pre_subsidy_childcare_expenses"), + ("SPM_CAPWKCCXPNS", "spm_unit_capped_work_childcare_expenses"), + ): + if raw in person.columns: + agg = ( + pd.to_numeric(person[raw], errors="coerce") + .fillna(0.0) + .groupby(person["person_spm_unit_id"]) + .max() + ) + spm[pe_name] = ( + spm["spm_unit_id"].map(agg).fillna(0.0).astype(float) + ) + if "SPM_ENGVAL" in person.columns: + eng = ( + pd.to_numeric(person["SPM_ENGVAL"], errors="coerce") + .fillna(0.0) + .groupby(person["person_spm_unit_id"]) + .max() + ) + spm["spm_unit_energy_subsidy"] = ( + spm["spm_unit_id"].map(eng).fillna(0.0).astype(float) + ) + + class _Key: + def __init__(self, value: str): + self.value = value + + # Group-owned id columns must exist only on their group tables; the + # person table carries them as join keys until this point. + person = person.drop( + columns=[ + c + for c in ( + "household_id", + "tax_unit_id", + "spm_unit_id", + "family_id", + "marital_unit_id", + "household_weight", + ) + if c in person.columns + ] + ) + + units_tu_new = ( + unit_map.rename(columns={"orig_tax_unit_id": "tax_unit_id"}) + .merge(units_tu, on="tax_unit_id", how="left") + .drop(columns=["tax_unit_id", "_half"]) + .rename(columns={"new_tax_unit_id": "TAX_ID"}) + ) + # ---- v3: support guard anchored to the PUF's OWN realized ranges ------- + # (The v2 guard clipped to the eCPS baseline's ranges — an eCPS + # contamination and, per the architecture review, the wrong reference. + # Structural heavy-tail control belongs to signed calibration targets; + # this clip only enforces the donor's own support.) + from microplex.data_sources.puf import PUF_VARIABLE_MAP as PUF_FIELD_MAP + from microplex.data_sources.puf import UPRATING_FACTORS + + _puf_csv = Path.home() / ".cache" / "microplex" / "puf_2015.csv" + _header = pd.read_csv(_puf_csv, nrows=0).columns + _by_upper = {c.upper(): c for c in _header} + _want = { + _by_upper[k.upper()]: v + for k, v in PUF_FIELD_MAP.items() + if k != "S006" and k.upper() in _by_upper + } + _puf_raw = pd.read_csv(_puf_csv, usecols=list(_want)) + _puf_raw.columns = [c for c in _puf_raw.columns] + PUF_FIELD_MAP = {c: _want[c] for c in _want} + _ranges: dict[str, tuple[float, float]] = {} + for _raw, _donor in PUF_FIELD_MAP.items(): + if _raw == "S006" or _raw not in _puf_raw.columns: + continue + _up = UPRATING_FACTORS.get(_donor, 1.0) + _v = pd.to_numeric(_puf_raw[_raw], errors="coerce").dropna() + _pe = DONOR_TO_PE.get(_donor, _donor) + _lo, _hi = float(_v.min()) * _up, float(_v.max()) * _up + _ranges[_pe] = (min(_lo, _hi), max(_lo, _hi)) + if "_estate_income_gross" in _puf_raw.columns or True: + # estate_income = E26390 - E26400 (uprated): bound by the rowwise diff. + try: + _est = pd.read_csv( + Path.home() / ".cache" / "microplex" / "puf_2015.csv", + usecols=["E26390", "E26400"], + ) + _diff = ( + pd.to_numeric(_est["E26390"], errors="coerce").fillna(0) + - pd.to_numeric(_est["E26400"], errors="coerce").fillna(0) + ) * UPRATING_FACTORS.get("estate_income", 1.0) + _ranges["estate_income"] = (float(_diff.min()), float(_diff.max())) + except Exception: + pass + for _c in list(person.columns): + if _c.startswith("person_") or _c not in _ranges: + continue + _vals = pd.to_numeric(person[_c], errors="coerce") + if _vals.isna().all(): + continue + _lo, _hi = _ranges[_c] + _clipped = _vals.clip(_lo, _hi) + _n = int((_clipped != _vals).sum()) + if _n: + log(f" puf-support-guard {_c}: clipped {_n} to [{_lo:,.0f}, {_hi:,.0f}]") + person[_c] = _clipped + + # ---- v3: AOTC factual inputs from the PUF tuition signal ---------------- + # usdata extended_cps: with no credit signal, the AOTC student mask is + # simply qualified_tuition_expenses > 0; the five factual eligibility + # flags are set for those students. + _aotc = ( + pd.to_numeric( + person.get("qualified_tuition_expenses", 0), errors="coerce" + ).fillna(0) + > 0 + ) + for _flag in ( + "is_pursuing_credential_for_american_opportunity_credit", + "attends_eligible_educational_institution_for_american_opportunity_credit", + "is_enrolled_at_least_half_time_for_american_opportunity_credit", + "has_american_opportunity_credit_1098_t_or_exception", + "has_american_opportunity_credit_institution_ein", + ): + person[_flag] = _aotc.astype(bool) + log(f" AOTC factual inputs: {_aotc.mean()*100:.1f}% of persons flagged") + # QRF-drawn boolean flags can land as mixed-object columns; normalize to + # clean bools so HDF serialization and PE casting are deterministic. + for _bcol in ( + "business_is_sstb", + "sstb_self_employment_income_would_be_qualified", + "self_employment_income_would_be_qualified", + ): + if _bcol in person.columns: + person[_bcol] = ( + pd.to_numeric(person[_bcol], errors="coerce").fillna(0) > 0.5 + ).astype(bool) + + # ---- v2: place variables at their PolicyEngine entity ------------------ + # The PUF and donor stages leave tax-unit and SPM-entity amounts on the + # person/household frames (head-carried); PE rejects inputs stored at the + # wrong entity length, so move each to its owning table. + from policyengine_us.system import system as _pe_system + + def _pe_entity(col: str) -> str | None: + var = _pe_system.variables.get(col) + return var.entity.key if var is not None else None + + tu_moves = [ + c for c in person.columns + if not c.startswith("person_") and _pe_entity(c) == "tax_unit" + ] + # Aggregate person-stored tax-unit amounts per unit id; applied onto the + # BUILT tax-unit table below (whose ids derive from persons, so the + # aggregation covers every unit — attaching to units_tu_new would leave + # NaNs for tail-stratum units absent from the spine unit map). + tu_agg = { + c: ( + pd.to_numeric(person[c], errors="coerce") + .fillna(0.0) + .groupby(person["person_tax_unit_id"]) + .sum() + ) + for c in tu_moves + } + person = person.drop(columns=tu_moves) + if tu_moves: + log(f" moved to tax_unit entity: {tu_moves}") + person_spm_moves = [ + c for c in person.columns + if not c.startswith("person_") and _pe_entity(c) == "spm_unit" + ] + spm_agg = { + c: person[c].groupby(person["person_spm_unit_id"]).first() + for c in person_spm_moves + } + person = person.drop(columns=person_spm_moves) + if person_spm_moves: + log(f" moved person->spm_unit entity: {person_spm_moves}") + spm_moves = [c for c in hh.columns if _pe_entity(c) == "spm_unit"] + for c in spm_moves: + val = dict(zip(hh["household_id"], hh[c])) + per_person = person["person_household_id"].map(val) + agg = per_person.groupby(person["person_spm_unit_id"]).first() + spm[c] = spm["spm_unit_id"].map(agg).fillna(0.0).astype(float) + hh = hh.drop(columns=[c]) + if spm_moves: + log(f" moved to spm_unit entity: {spm_moves}") + for c, agg in spm_agg.items(): + mapped = spm["spm_unit_id"].map(agg) + if pd.api.types.is_bool_dtype(person.dtypes.get(c, bool)) or mapped.dtype == object: + spm[c] = mapped.fillna(False).astype(bool) + else: + spm[c] = pd.to_numeric(mapped, errors="coerce").fillna(0.0) + + tax_unit_tbl = unit_table("tax_unit_id", units_tu_new) + for c, agg in tu_agg.items(): + tax_unit_tbl[c] = ( + tax_unit_tbl["tax_unit_id"].map(agg).fillna(0.0).astype(float) + ) + entity_frames = { + _Key("person"): person, + _Key("household"): hh, + _Key("tax_unit"): tax_unit_tbl, + _Key("spm_unit"): spm, + _Key("family"): unit_table("family_id"), + _Key("marital_unit"): unit_table("marital_unit_id"), + } + + # ---- Stage G: export ---------------------------------------------------- + log("stage G: export") + contract = ExportContract.from_path( + REPO / "packs/us/manifests/ecps_export_contract.json" + ) + defaults = json.loads( + (REPO / "packs/us/manifests/export_defaults.json").read_text() + ) + defaults.pop("_source", None) + for column, value in V1_ZERO_DEFAULTS.items(): + defaults.setdefault(column, value) + candidate_h5 = out / "candidate_policyengine_us.h5" + gate = export_policyengine_us_dataset( + entity_frames, + period=args.calendar_year, + output_path=candidate_h5, + contract=contract, + defaults=defaults, + allow_incomplete=smoke, + ) + (out / "export_gate.json").write_text(json.dumps(gate.to_dict(), indent=2)) + # Sibling export in the eCPS time-period layout ({variable}/{period} + # datasets) — the format the comparison harness and HF artifacts use. + import h5py + + tp_h5 = out / "candidate_timeperiod.h5" + allowed = set(contract.required) | set(contract.optional) + from policyengine_us.data import USSingleYearDataset + + saved = USSingleYearDataset(file_path=str(candidate_h5)) + saved_tables = [ + saved.person, + saved.household, + saved.tax_unit, + saved.spm_unit, + saved.family, + saved.marital_unit, + ] + with h5py.File(tp_h5, "w") as handle: + seen: set[str] = set() + for frame in saved_tables: + if frame is None or len(frame) == 0: + continue + for column in frame.columns: + if column in seen or column not in allowed: + continue + seen.add(column) + values = frame[column].to_numpy() + if values.dtype.kind in {"U", "O"}: + values = values.astype("S") + elif values.dtype.kind == "b": + values = values.astype(bool) + grp = handle.create_group(column) + grp.create_dataset(str(args.calendar_year), data=values) + missing_tp = sorted(set(contract.required) - seen) + log(f" time-period export: {tp_h5.name} cols={len(seen)} missing={len(missing_tp)}") + log( + f" gate passed={gate.passed} missing={len(gate.missing_required)} " + f"defaulted={len(gate.defaulted)} dropped={len(gate.dropped)}" + ) + if gate.missing_required: + log(f" missing (first 25): {list(gate.missing_required)[:25]}") + + # ---- Stage H: score ----------------------------------------------------- + if args.score and candidate_h5.exists(): + log("stage H: sound eCPS comparison (legacy harness)") + cmd = [ + str(OLD_WORKTREE / ".venv/bin/python"), + "-m", + "microplex_us.pipelines.ecps_replacement_comparison", + "--candidate-dataset", + str(out / "candidate_timeperiod.h5"), + "--baseline-dataset", + str(args.baseline_h5), + "--output-dir", + str(out / "sound_comparison"), + "--period", + str(args.calendar_year), + "--force", + "--policyengine-us-data-repo", + str(args.usdata_repo), + "--policyengine-us-data-python", + str(args.usdata_repo / ".venv/bin/python"), + ] + log(" " + " ".join(cmd)) + proc = subprocess.run(cmd, cwd=OLD_WORKTREE, capture_output=True, text=True) + (out / "score_stdout.log").write_text(proc.stdout) + (out / "score_stderr.log").write_text(proc.stderr) + log(f" harness exit: {proc.returncode}") + result_json = out / "sound_comparison" / "sound_ecps_replacement_comparison.json" + if result_json.exists(): + payload = json.loads(result_json.read_text()) + log(json.dumps(payload.get("headline", payload), indent=2)[:2000]) + return proc.returncode + + return 0 if (gate.passed or smoke) else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/packages/populace-data/build/us/check_parity.py b/packages/populace-data/build/us/check_parity.py new file mode 100644 index 00000000..f4f2aa43 --- /dev/null +++ b/packages/populace-data/build/us/check_parity.py @@ -0,0 +1,181 @@ +"""Parity + smoke gate for the built artifact, via populace.build gates. + +Parity is judged at SIMULATION level over the reference's stored input +layers: for every variable the reference (enhanced CPS) stores and populates, +the candidate's simulation must produce a non-zero layer too (stored or +computed by engine formulas — dropping a formula-masking zero column and +letting the engine compute is parity, not a gap). + +Smoke prints the headline aggregates (population, SNAP, net worth, net STCG, +investment interest) for the publish decision. +""" + +import json +import sys +from pathlib import Path + +import h5py +import numpy as np + +ART = Path.home() / ".claude-worktrees" / "microplex-spec-build" / "artifacts" +CANDIDATE = ART / "populace_us_2024.h5" +REFERENCE = Path.home() / "populace-score-work" / "enhanced_cps_2024_hf_main.h5" +OUT = ART / "parity_gate.json" +YEAR = 2024 + + +def stored_layers(path: Path) -> dict[str, float]: + """entity-unprefixed variable -> nonzero share among stored columns.""" + shares: dict[str, float] = {} + with h5py.File(path) as f: + def visit(name, node): + if not isinstance(node, h5py.Dataset): + return + var = name.split("/")[-2] if name.endswith(str(YEAR)) else name.split("/")[-1] + # layouts: "entity/var" (single-year) or "var/2024" (timeperiod) + parts = name.split("/") + if len(parts) == 2 and parts[1] == str(YEAR): + var = parts[0] + elif len(parts) == 2: + var = parts[1] + else: + var = parts[-1] + values = node[:] + if values.dtype.kind in ("S", "O", "U"): + return + shares[var] = float((np.asarray(values, dtype=np.float64) != 0).mean()) + f.visititems(visit) + return shares + + +def candidate_stored_shares() -> dict[str, float]: + """entity.column -> nonzero share over every stored candidate column.""" + from policyengine_us.data import USSingleYearDataset + + ds = USSingleYearDataset(file_path=str(CANDIDATE)) + shares: dict[str, float] = {} + for ent in ("person", "household", "tax_unit", "spm_unit", "family", + "marital_unit"): + tbl = getattr(ds, ent) + structural = {f"{ent}_id", f"{ent}_weight"} | { + c for c in tbl.columns if c.startswith("person_") + } + for c in tbl.columns: + if c in structural: + continue + vals = tbl[c].to_numpy() + if vals.dtype.kind not in "fiub": + continue + shares[f"{ent}.{c}"] = float( + (np.asarray(vals, dtype=np.float64) != 0).mean() + ) + return shares + + +def main() -> int: + from policyengine_us import Microsimulation + from populace.build import exported_nonzero_gate, parity_gate + + # gate 0: every stored column carries signal (populate it or drop it) + nonzero = exported_nonzero_gate(candidate_stored_shares()) + print( + f"exported_nonzero: passed={nonzero.passed} " + f"({nonzero.details['columns_checked']} stored columns)" + ) + for line in nonzero.failures: + print(f" ZERO {line}") + + ref_layers = stored_layers(REFERENCE) + print(f"reference stored layers: {len(ref_layers)}") + + sim = Microsimulation(dataset=str(CANDIDATE)) + tbs = sim.tax_benefit_system + candidate_shares: dict[str, float] = {} + reference_shares: dict[str, float] = {} + skipped: list[str] = [] + for var, ref_share in sorted(ref_layers.items()): + if var not in tbs.variables: + skipped.append(var) + continue + if tbs.variables[var].definition_period not in ("year",): + skipped.append(var) + continue + try: + values = np.asarray( + sim.calculate(var, YEAR).values, dtype=np.float64 + ) + except Exception as error: # noqa: BLE001 - report, never mask + candidate_shares[var] = 0.0 + reference_shares[var] = ref_share + print(f" calc failed {var}: {type(error).__name__} {error}") + continue + candidate_shares[var] = float((values != 0).mean()) + reference_shares[var] = ref_share + + # weights are structural, not layers + for structural in ("household_weight", "person_weight"): + candidate_shares.pop(structural, None) + reference_shares.pop(structural, None) + + result = parity_gate(candidate_shares, reference_shares) + print( + f"parity: passed={result.passed} gaps={result.details['gaps']} " + f"(checked {result.details['reference_populated_layers']} populated " + f"reference layers, skipped {len(skipped)} non-annual/unknown)" + ) + for line in result.failures: + print(f" GAP {line}") + + # ---- smoke aggregates ------------------------------------------------- + def total(var: str) -> float: + return float(sim.calculate(var, YEAR).sum()) + + smoke = { + # weighted person count via microdf's weighted .count() + "people_m": float(sim.calculate("age", YEAR).count()) / 1e6, + "snap_b": total("snap") / 1e9, + "net_worth_t": total("net_worth") / 1e12, + "net_stcg_b": total("short_term_capital_gains") / 1e9, + "investment_interest_expense_b": total("investment_interest_expense") + / 1e9, + "tips_b": total("tip_income") / 1e9, + "pre_subsidy_rent_b": total("pre_subsidy_rent") / 1e9, + } + print("smoke:", json.dumps(smoke, indent=1)) + + # Telemetry (fail-soft): gate verdicts as queryable rows. + try: + import sys as _sys + + _sys.path.insert(0, str(Path(__file__).resolve().parent)) + import populace_telemetry as _telemetry + + _telemetry.push_gate_result(nonzero) + _telemetry.push_gate_result(result) + except Exception as _err: # noqa: BLE001 + print(f"telemetry skipped: {_err}") + + def _gate_dict(gate): + return { + "passed": gate.passed, + "failures": list(gate.failures), + "details": dict(gate.details), + } + + OUT.write_text( + json.dumps( + { + "exported_nonzero": _gate_dict(nonzero), + "parity": _gate_dict(result), + "smoke": smoke, + "skipped_layers": skipped, + }, + indent=1, + ) + ) + print(f"wrote {OUT}") + return 0 if (result.passed and nonzero.passed) else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/packages/populace-data/build/us/enrich_artifact.py b/packages/populace-data/build/us/enrich_artifact.py new file mode 100644 index 00000000..8263ba76 --- /dev/null +++ b/packages/populace-data/build/us/enrich_artifact.py @@ -0,0 +1,150 @@ +"""Post-calibration enrichment: simulation-dependent factual inputs. + +Three eCPS layers need a baseline simulation on the finished artifact (their +rules read eligibility or credit aggregates), so they are written after the +calibrated artifact exists — mirroring usdata's own derivations exactly: + +- ``takes_up_housing_assistance_if_eligible`` (SPM grain): reported + recipients first, then a seeded fill to the national receipt rate among + the simulation-eligible. +- ``is_pregnant`` (person): seeded state-rate draws for women 15-44 + (CDC/Census rates via usdata's etl). +- ``would_file_taxes_voluntarily`` (tax unit): usdata's demographic rate + table over children/wage/age bins, excluding EITC claimants. + +Weights are untouched; the script verifies them byte-identical after save. +""" + +import sys +from pathlib import Path + +import numpy as np + +USDATA = Path.home() / ".claude-worktrees" / "usdata-populace" +sys.path.insert(0, str(USDATA)) + +ART = Path.home() / ".claude-worktrees" / "microplex-spec-build" / "artifacts" +OUT = ART / "populace_us_2024.h5" +OUT_TP = ART / "populace_us_2024_timeperiod.h5" +YEAR = 2024 + + +def main() -> int: + from policyengine_us import Microsimulation + from policyengine_us.data import USSingleYearDataset + from policyengine_us_data.datasets.cps.cps import ( + _voluntary_filing_age_bin, + _voluntary_filing_children_bin, + _voluntary_filing_rate_by_tax_unit, + _voluntary_filing_wage_income_bin, + ) + from policyengine_us_data.datasets.cps.takeup import ( + prioritize_reported_recipients, + ) + from policyengine_us_data.db.etl_pregnancy import get_state_pregnancy_rates + from policyengine_us_data.parameters import load_take_up_rate + from policyengine_us_data.utils.randomness import seeded_rng + + ds = USSingleYearDataset(file_path=str(OUT)) + weights_before = np.asarray(ds.household["household_weight"], dtype=np.float64) + sim = Microsimulation(dataset=str(OUT)) + + # --- housing assistance take-up (SPM grain) ---------------------------- + rate = load_take_up_rate("housing_assistance", YEAR) + reported = np.asarray( + ds.spm_unit["receives_housing_assistance"], dtype=bool + ) + eligible = np.asarray( + sim.calculate("is_eligible_for_housing_assistance", YEAR).values, + dtype=bool, + ) + rng = seeded_rng("takes_up_housing_assistance_if_eligible") + takes_up = prioritize_reported_recipients( + reported, rate, rng.random(len(reported)), eligible_mask=eligible + ) + ds.spm_unit["takes_up_housing_assistance_if_eligible"] = np.asarray( + takes_up, dtype=bool + ) + print( + f"takes_up_housing_assistance_if_eligible: rate={rate:.3f} " + f"reported nz={reported.mean()*100:.1f}% -> takeup nz={np.mean(takes_up)*100:.1f}%" + ) + + # --- pregnancy (person) ------------------------------------------------- + rates = get_state_pregnancy_rates(cdc_year=YEAR, acs_year=YEAR) + national = 0.041 + state = np.asarray( + sim.calculate("state_fips", YEAR, map_to="person").values, dtype=int + ) + by_person = np.array([rates.get(int(s), national) for s in state]) + age = np.asarray(ds.person["age"], dtype=float) + is_female = np.asarray(ds.person["is_female"], dtype=bool) + eligible_preg = is_female & (age >= 15) & (age <= 44) + rng = seeded_rng("is_pregnant") + ds.person["is_pregnant"] = eligible_preg & ( + rng.random(len(age)) < by_person + ) + print(f"is_pregnant: nz={np.mean(ds.person['is_pregnant'])*100:.2f}%") + + # --- voluntary filing (tax unit) ---------------------------------------- + voluntary_rates = load_take_up_rate("voluntary_filing", YEAR) + takes_up_eitc = np.asarray(ds.tax_unit["takes_up_eitc"], dtype=bool) + eitc = np.asarray(sim.calculate("eitc", YEAR).values, dtype=float) + claims_eitc = takes_up_eitc & (eitc > 0) + children = np.asarray( + sim.calculate("tax_unit_child_dependents", YEAR).values, dtype=float + ) + wage = np.asarray( + sim.calculate("employment_income", YEAR, map_to="tax_unit").values, + dtype=float, + ) + age_head = np.asarray(sim.calculate("age_head", YEAR).values, dtype=float) + rate_by_unit = _voluntary_filing_rate_by_tax_unit( + voluntary_rates, + _voluntary_filing_children_bin(children), + _voluntary_filing_wage_income_bin(wage), + _voluntary_filing_age_bin(age_head), + ) + rng = seeded_rng("would_file_taxes_voluntarily") + ds.tax_unit["would_file_taxes_voluntarily"] = ~claims_eitc & ( + rng.random(len(rate_by_unit)) < np.asarray(rate_by_unit, dtype=float) + ) + print( + "would_file_taxes_voluntarily: " + f"nz={np.mean(ds.tax_unit['would_file_taxes_voluntarily'])*100:.2f}%" + ) + + ds.save(str(OUT)) + check = USSingleYearDataset(file_path=str(OUT)) + weights_after = np.asarray( + check.household["household_weight"], dtype=np.float64 + ) + assert np.array_equal(weights_before, weights_after), "weights changed!" + + # propagate the three new layers into the timeperiod export + import h5py + + with h5py.File(OUT_TP, "a") as f: + for var, values in ( + ( + "takes_up_housing_assistance_if_eligible", + np.asarray( + check.spm_unit["takes_up_housing_assistance_if_eligible"] + ), + ), + ("is_pregnant", np.asarray(check.person["is_pregnant"])), + ( + "would_file_taxes_voluntarily", + np.asarray(check.tax_unit["would_file_taxes_voluntarily"]), + ), + ): + key = f"{var}/{YEAR}" + if key in f: + del f[key] + f.create_dataset(key, data=values) + print("enriched + verified (weights byte-identical)") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/packages/populace-data/build/us/extract_target_surface.py b/packages/populace-data/build/us/extract_target_surface.py index e349ca8b..e7097d61 100644 --- a/packages/populace-data/build/us/extract_target_surface.py +++ b/packages/populace-data/build/us/extract_target_surface.py @@ -1,44 +1,27 @@ -"""Re-extract the PE-native surface, saving the UNSCALED (raw) matrix + targets -so calibration uses the eCPS +1 loss on raw magnitudes (not the comparison's -normalized ~0.01 scale, where the +1 regularizer breaks).""" -import os -import sys - +import os, sys import numpy as np - -# The scoring-harness checkout that provides the extraction entrypoint; set -# SCORING_HARNESS_SRC to its src/ directory to re-run this snapshot. -sys.path.insert(0, os.environ["SCORING_HARNESS_SRC"]) +sys.path.insert(0, os.path.expanduser("~/CosilicoAI/microplex-us/src")) from microplex_us.pipelines.ecps_replacement_comparison import _extract_pe_native_loss_inputs from pathlib import Path - -cand = "/Users/maxghenis/.claude-worktrees/microplex-spec-build/artifacts/spec_candidate_full_2024_v5_beats/candidate_timeperiod.h5" -repo = Path("/Users/maxghenis/.claude-worktrees/usdata-f7458313") -print("extracting (raw)...", flush=True) +ART = Path.home()/".claude-worktrees"/"microplex-spec-build"/"artifacts" +cand = ART/"spec_candidate_full_2024"/"candidate_timeperiod.h5" +repo = Path.home()/".claude-worktrees"/"usdata-populace" +print("extracting v2 raw surface...", flush=True) inp = _extract_pe_native_loss_inputs( - input_dataset_path=cand, period=2024, + input_dataset_path=str(cand), period=2024, policyengine_us_data_repo=repo, policyengine_us_data_python=repo/".venv/bin/python", skip_tax_expenditure_targets=False, ) -A_s = np.asarray(inp["scaled_matrix"], np.float64) # (n_hh, n_targets) scaled +A_s = np.asarray(inp["scaled_matrix"], np.float64) b_s = np.asarray(inp["scaled_target"], np.float64) -w0 = np.asarray(inp["initial_weights"], np.float64) -b_raw = inp.get("unscaled_target") -scaling = inp.get("scaling") +w0 = np.asarray(inp["initial_weights"], np.float64) +b_raw = np.asarray(inp["unscaled_target"], np.float64) +scaling = np.asarray(inp["scaling"], np.float64) names = [str(x) for x in inp["metadata"]["target_names"]] -print("has unscaled_target:", b_raw is not None, "| has scaling:", scaling is not None, flush=True) - -if b_raw is not None and scaling is not None: - b_raw = np.asarray(b_raw, np.float64); scaling = np.asarray(scaling, np.float64) - # scaled = raw * scaling => raw_A = scaled_A / scaling (per target/column) - A_raw = A_s / scaling[None, :] - # verify: raw estimate at w0 ~ raw target - est = A_raw.T @ w0 - print("raw b: median %.4g max %.4g | est/b median %.3f" % (np.median(b_raw), b_raw.max(), np.median((est+1)/(b_raw+1))), flush=True) -else: - raise SystemExit("no unscaled arrays returned; need a different extraction path") - -out = "/Users/maxghenis/.claude-worktrees/microplex-spec-build/artifacts/v5_target_surface_raw.npz" +A_raw = A_s / scaling[None, :] +est = A_raw.T @ w0 +print(f"targets {len(names)} | b median {np.median(b_raw):.3g} | est/b median {np.median((est+1)/(b_raw+1)):.3f}", flush=True) +out = ART/"target_surface_raw.npz" np.savez_compressed(out, A=A_raw, b=b_raw, w0=w0, names=np.array(names, dtype=object)) -loss = float(np.mean(((A_raw.T@w0 - b_raw + 1.0)/(b_raw + 1.0))**2)) -print("SAVED", out, "| raw initial loss %.4f (should be small, ~eCPS-comparable)" % loss, flush=True) +loss = float(np.mean(((A_raw.T@w0 - b_raw)/(b_raw + 1.0))**2)) +print(f"SAVED {out} | raw initial loss {loss:.4f}", flush=True) diff --git a/packages/populace-data/build/us/hf_dataset_card.md b/packages/populace-data/build/us/hf_dataset_card.md index ce7f031c..ccf45d6f 100644 --- a/packages/populace-data/build/us/hf_dataset_card.md +++ b/packages/populace-data/build/us/hf_dataset_card.md @@ -13,9 +13,11 @@ tags: The **populace-built US population**: a calibrated synthetic microdataset for [PolicyEngine-US](https://github.com/PolicyEngine/policyengine-us), built by the -[`populace`](https://github.com/PolicyEngine/populace) stack. It loads anywhere -the enhanced CPS loads (an API-compatible alternative population), with its own -calibrated weights — and its own strengths and gaps, both documented below. +[`populace`](https://github.com/PolicyEngine/populace) stack **entirely from +primary sources** — the enhanced CPS appears only as the benchmark this file is +scored against, never as a build input. It loads anywhere the enhanced CPS +loads (an API-compatible alternative population), with its own calibrated +weights — and its own strengths and gaps, both documented below. ## Load it @@ -43,21 +45,47 @@ path = hf_hub_download( ) ``` -## What it is - -One HDF5 `USSingleYearDataset` per year. The population is generated from a -declarative spec: Current Population Survey ASEC provides household structure, -demographics, benefits, and tenure (~half the records are PUF-derived clones, -flagged per record); tax detail is imputed from the IRS Public Use File with -weight-aware quantile-forest models; the wealth, mortgage, vehicle, insurance --premium, and prior-year-income layers are imputed with the same models using -the published enhanced CPS as the donor (those layers are survey-imputed in -the incumbent itself); every imputed value is clipped to the incumbent's -realized per-record range (the support guard); and the result is calibrated to -PolicyEngine's administrative target surface (3,704 IRS/Census/program -targets) with a hard per-record weight bound (`max_weight_ratio=50`), so no -aggregate leans on a handful of super-weighted records. Same source classes -and hosting precedent as PolicyEngine's published enhanced CPS. +## How it is built + +One HDF5 `USSingleYearDataset` per year. Every layer comes from a primary +survey or administrative source: + +| source | provides | +| --- | --- | +| Census CPS ASEC | household structure, demographics, incomes, benefits, tenure, hours, occupation flags, health coverage at interview, retirement distributions (DST codes), childcare, prior-year income (longitudinal PERIDNUM join) | +| IRS SOI Public Use File 2015 (uprated) | tax detail: capital gains, dividends, interest, itemized-deduction inputs, QBI/SSTB components, partnership self-employment, estates, tuition | +| Fed SCF 2022 | wealth: bank/stock/bond assets, net worth, mortgage balance hints | +| Census SIPP | tip income for tipped occupations; household vehicles (count and value) | +| CPS-ORG | hourly wage, paid-hourly status, union coverage | +| MEPS-IC parameters | employer-sponsored insurance premiums | +| Census ACS 2022 | rent for renter households | + +Imputations use weight-aware quantile-forest models fit on each donor's own +records, and every imputed value is clipped to **that donor's** realized range +(the support guard) — nothing is anchored to the enhanced CPS. The result is +calibrated to PolicyEngine's administrative target surface (3,704 IRS/Census/ +program targets, **plus a signed net short-term capital gains target** so the +optimizer cannot silently drive a net-negative aggregate to extremes) with a +hard per-record weight bound (`max_weight_ratio=50`), so no aggregate leans on +a handful of super-weighted records. + +## Acceptance gates + +The build refuses to publish unless every gate passes; this file passed all +of them: + +- **Parity 0**: every PolicyEngine input layer the enhanced CPS populates + non-degenerately, this file's simulation populates (169 reference layers + checked at simulation level). +- **Exported-nonzero**: all 309 stored columns carry signal — no all-zero + scaffolding that would silently mask engine formulas or defaults. +- **Calibration**: 95.09% of 3,704 targets within 10% (loss 0.022); max + household weight 297,651 with **zero records above 500k** (the enhanced CPS + ships 21, max 1.05M). +- **Smoke aggregates** through `Microsimulation`: 332.3M people, $97.8B SNAP, + $176.5T net worth (Fed Z.1 ≈ $169T), net short-term capital gains + **−$77.5B** against the −$76.8B PUF-anchored target, tips $52.9B, rent + $757.5B. ## Validation @@ -67,47 +95,28 @@ targets never seen by either side's refit. Lower is better. | metric | populace-us | enhanced CPS | | --- | --- | --- | -| training loss (2,965 targets) | **0.132** | 1.089 | -| held-out loss (739 unseen targets) | **0.032** | 0.317 | -| full-surface loss (3,704 targets) | **0.164** | 1.406 | - -Per individual target the incumbent still wins more often (2,484 of 3,704 to -our 1,168, 52 ties): populace wins big where it wins and loses narrowly where -it loses. Both facts are the story. - -Shipped-file properties: **0 parity gaps** (every PolicyEngine input the -enhanced CPS populates non-degenerately, this file populates), **95.55% of -3,704 calibration targets within 10%** (calibration loss 0.022), max household -weight 382,478 with **zero records above 500k** (the enhanced CPS ships 21, -max 1.05M). End-to-end through `Microsimulation`: 332.7M people, $93.2B SNAP, -$338.4B traditional 401(k) contributions, $163.6T net worth (incumbent: -$163.4T). +| training loss (2,965 targets) | **0.190** | 1.089 | +| held-out loss (739 unseen targets) | **0.038** | 0.317 | +| full-surface loss (3,704 targets) | **0.228** | 1.406 | + +Per individual target the incumbent still wins more often (2,613 of 3,704 to our 1,040, 51 ties): populace wins big where it wins and loses narrowly where it loses. Both facts are the story. ## Known gaps We publish the misses with the hits: -- **Short-term capital gains over-weight large losses**: the aggregate is - −$0.9T against the donor's weighted −$77B. The pool's records are faithful - at design weights (−$164B); calibration amplifies loss-heavy records to hit - the targets it can see, and net STCG is not on the target surface. A - net-STCG calibration target is on the roadmap. -- **Donor-imputed layers inherit the incumbent's model error.** Wealth, - mortgages, vehicles, premiums, and prior-year income are drawn from models - trained on the enhanced CPS, whose own values for those layers are - themselves survey-imputed (SCF/SIPP/ACS). -- **Aggregate household net income is $13.7T** vs the incumbent's $22.2T — - most of the gap is the STCG item above plus thinner tail capital income; - several program totals land closer to administrative actuals than the - incumbent (SNAP $93.2B, SSI near the ~$60B outlay). Results from the two - populations are not interchangeable. -- **Per-target wins still favor the incumbent** (see Validation): the - aggregate losses are far lower, but on a majority of individual targets the - enhanced CPS sits closer. +- **Net worth runs ~4% above Fed Z.1** ($176.5T vs ≈ $169T): the calibration + target ($160T) sits below Z.1 and the achieved total lands between them. +- **Investment interest expense is thin** ($5.1B against IRS SOI ≈ $24B): + the PUF-residual rule populates the layer conservatively; a dedicated SOI + calibration target is the roadmap item. +- **Per-target wins vs the incumbent**: see Validation — aggregate losses + are what the comparison gates on, but per-target patterns differ between + the two populations. Results are not interchangeable. The dashboard at [populace.dev/dashboard](https://populace.dev/dashboard) -shows the full per-family calibration fit, the worst-fit targets by name, and -the weight distribution. Methodology and evidence: -[populace.dev](https://populace.dev); loader and registry: +shows the full per-family calibration fit, the worst-fit targets by name, the +weight distribution, and a live strip while a build chain runs. Methodology +and evidence: [populace.dev](https://populace.dev); loader and registry: [github.com/PolicyEngine/populace](https://github.com/PolicyEngine/populace) (`packages/populace-data`). diff --git a/packages/populace-data/build/us/primary_source_impute.py b/packages/populace-data/build/us/primary_source_impute.py new file mode 100644 index 00000000..fd2101ad --- /dev/null +++ b/packages/populace-data/build/us/primary_source_impute.py @@ -0,0 +1,601 @@ +"""Primary-source imputation stages (v3, eCPS-free). + +Replaces ecps_donor_impute.py: every layer draws from its primary survey via +the usdata loaders in the worktree (Fed SCF for wealth, SIPP for tips, CPS-ORG +for hourly wage, MEPS-IC parameters for ESI premiums). The enhanced CPS +appears nowhere — it is only ever the benchmark in scoring. Each imputed +block is support-guarded to ITS OWN donor's realized per-record range. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +#: SCF-sourced wealth block (household grain in the pool; donor is the +#: summarized Fed SCF with CPS-comparable predictor names). +SCF_TARGETS = [ + "net_worth", + "scf_primary_residence_value", + "scf_retirement_assets", + "scf_business_equity", + "scf_mortgage_debt", + "scf_other_residential_real_estate", + "scf_nonresidential_real_estate_equity", + "scf_other_residential_debt", + "scf_other_financial_assets", + "scf_other_nonfinancial_assets", + "scf_other_managed_assets", + "scf_cash_value_life_insurance", + "scf_certificates_of_deposit", + "scf_savings_bonds", + "scf_credit_card_debt", + "scf_student_loan_debt", + "scf_vehicle_installment_debt", + "scf_other_installment_debt", + "scf_other_lines_of_credit", + "scf_other_debt", + "household_vehicles_owned", + "household_vehicles_value", + "auto_loan_balance", + "auto_loan_interest", + "bank_account_assets", + "bond_assets", + "stock_assets", +] + + +def _support_guard(values: np.ndarray, donor: np.ndarray, name: str, log) -> np.ndarray: + lo, hi = float(np.nanmin(donor)), float(np.nanmax(donor)) + clipped = np.clip(values, lo, hi) + n = int((clipped != values).sum()) + if n: + log(f" support-guard {name}: clipped {n} to donor range [{lo:,.0f}, {hi:,.0f}]") + return clipped + + +def add_scf_wealth(person: pd.DataFrame, hh: pd.DataFrame, seed: int, log) -> pd.DataFrame: + """Impute the wealth block onto households from SCF 2022 (usdata blueprint). + + Mirrors usdata cps.py's own SCF stage: SCF_2022 donor, `wgt` weights, the + same predictor list, target lists from the same helper functions; imputed + at household-head grain and attached to households, support-guarded to the + SCF's own realized ranges. + """ + from microimpute import Imputer + from policyengine_us_data.datasets.cps.cps import ( + add_scf_financial_asset_targets, + add_scf_household_asset_targets, + add_scf_net_worth_component_targets, + add_scf_net_worth_target, + ) + from policyengine_us_data.datasets.scf.scf import SCF_2022 + + scf_raw = SCF_2022().load_dataset() + scf = pd.DataFrame({k: scf_raw[k] for k in scf_raw.keys()}) + targets = list( + dict.fromkeys( + list(add_scf_net_worth_target(scf)) + + ["auto_loan_balance", "auto_loan_interest"] + + list(add_scf_financial_asset_targets(scf)) + + list(add_scf_household_asset_targets(scf)) + + list(add_scf_net_worth_component_targets(scf)) + ) + ) + PREDICTORS = [ + "age", "is_female", "cps_race", "is_married", + "own_children_in_household", "employment_income", + "interest_dividend_income", "social_security_pension_income", + ] + log(f" SCF 2022 donor: {len(scf):,} rows, {len(targets)} targets") + + num = lambda c: pd.to_numeric(person.get(c, 0), errors="coerce").fillna(0) # noqa: E731 + if "is_household_head" not in person.columns: + raise RuntimeError( + "head-carry requires is_household_head on the person frame " + "(derive it from ASEC P_SEQ == 1 before stage F2); refusing to " + "head-carry onto an all-False mask." + ) + head = person["is_household_head"].astype(bool) + pf = pd.DataFrame( + { + "hh": person["person_household_id"], + "head": head, + "age": num("A_AGE") if "A_AGE" in person.columns else num("age"), + "is_female": person.get("is_female", False), + "cps_race": num("cps_race"), + "is_married": person.get("A_MARITL", pd.Series(0, index=person.index)).isin([1, 2]).astype(float) if "A_MARITL" in person.columns else 0.0, + "own_children_in_household": num("own_children_in_household"), + "employment_income": num("employment_income"), + "interest_dividend_income": num("taxable_interest_income") + num("dividend_income") + num("qualified_dividend_income") + num("non_qualified_dividend_income"), + "social_security_pension_income": num("social_security") + num("taxable_pension_income"), + } + ) + heads = pf[pf["head"]].drop_duplicates("hh").set_index("hh") + # Households with no flagged head: use the eldest member. + missing = set(hh["household_id"]) - set(heads.index) + if missing: + eldest = ( + pf[pf["hh"].isin(missing)] + .sort_values("age", ascending=False) + .drop_duplicates("hh") + .set_index("hh") + ) + heads = pd.concat([heads, eldest]) + recv = heads.reindex(hh["household_id"]).fillna(0.0) + + donor_cols = [c for c in PREDICTORS if c in scf.columns] + targets = [t for t in targets if t in scf.columns] + donor = scf[donor_cols + targets + ["wgt"]].dropna() + fitted = Imputer(seed=seed, log_level="WARNING").fit( + donor, donor_cols, targets, weight_col="wgt" + ) + draws = fitted.predict(recv[donor_cols].copy().reset_index(drop=True)) + hh = hh.copy() + for t in targets: + vals = np.asarray(draws[t], dtype=np.float64) + hh[t] = _support_guard(vals, scf[t].to_numpy(dtype=np.float64), t, log) + + def hh_sum(cols): + present = [c for c in cols if c in hh.columns] + return sum(hh[c].to_numpy(dtype=np.float64) for c in present) if present else None + + # The usdata target helpers impute PE-shaped names under an scf_ prefix + # (scf_bank_account_assets etc.); person-entity assets head-carry from + # those. Loud: a missing source column is a build bug, not a skip. + PE_FROM_SCF = { + "bank_account_assets": "scf_bank_account_assets", + "stock_assets": "scf_stock_assets", + "bond_assets": "scf_bond_assets", + } + head_carry_to_person = {} + for pe_name, scf_name in PE_FROM_SCF.items(): + if scf_name not in hh.columns: + raise RuntimeError( + f"SCF stage expected imputed column {scf_name!r} for " + f"{pe_name!r}; imputed scf columns: " + f"{[c for c in hh.columns if c.startswith('scf_')][:8]}..." + ) + head_carry_to_person[pe_name] = hh[scf_name].to_numpy(dtype=np.float64) + # net worth = sum of imputed SCF components (usdata computes it the same way). + # Net worth is usdata's own imputed measure — scf_net_worth already nets + # assets against debts. Summing the scf_ component columns ON TOP of it + # double-counts (the +34% miss on nation/net_worth/total); the vehicle + # value joins later in the SIPP vehicle stage, mirroring usdata's + # net_worth_components assembly. + if "scf_net_worth" not in hh.columns: + raise RuntimeError( + "SCF stage expected imputed column 'scf_net_worth'; got " + f"{[c for c in hh.columns if c.startswith('scf_')][:8]}..." + ) + hh["net_worth"] = hh["scf_net_worth"].to_numpy(dtype=np.float64) + # person-entity assets: head-carry onto persons (export entity mover is + # person->tax_unit only; person columns export directly). + if "is_household_head" not in person.columns: + raise RuntimeError( + "head-carry requires is_household_head on the person frame " + "(derive it from ASEC P_SEQ == 1 before stage F2); refusing to " + "head-carry onto an all-False mask." + ) + headp = person["is_household_head"].astype(bool) + hmap = dict(zip(hh["household_id"].tolist(), range(len(hh)))) + pidx = person["person_household_id"].map(hmap) + for pe_name, vals in head_carry_to_person.items(): + person[pe_name] = np.where( + headp & pidx.notna(), np.asarray(vals)[pidx.fillna(0).astype(int)], 0.0 + ) + log(f" SCF wealth block: {len(targets)} scf vars + net_worth + {sorted(head_carry_to_person)} (weighted=wgt)") + return person, hh + + +def add_sipp_tips(person: pd.DataFrame, log) -> pd.DataFrame: + """Tips from the SIPP-trained model (usdata get_tip_model).""" + from policyengine_us_data.datasets.sipp import get_tip_model + + model = get_tip_model() + x = pd.DataFrame(index=person.index) + emp = pd.to_numeric(person["employment_income"], errors="coerce").fillna(0) + x["employment_income"] = emp + x["is_tipped_occupation"] = person.get("is_tipped_occupation", False) + x["age"] = pd.to_numeric(person.get("age", person.get("A_AGE", 0)), errors="coerce").fillna(0) + # usdata's call site builds pension/retirement/non-SSI aggregates first; + # provide every feature the model declares, defaulting absent ones to 0. + needed = list(getattr(model, "predictors", []) or []) + for c in needed: + if c not in x.columns: + src = person.get(c) + x[c] = pd.to_numeric(src, errors="coerce").fillna(0) if src is not None else 0.0 + try: + person = person.copy() + person["tip_income"] = np.asarray( + model.predict(X_test=x, mean_quantile=0.5).tip_income.values + ) + person.loc[~person.get("is_tipped_occupation", pd.Series(False, index=person.index)).astype(bool), "tip_income"] = 0.0 + log(f" SIPP tips: nz {(person['tip_income']>0).mean()*100:.1f}%") + except Exception as exc: + log(f" SIPP tips FAILED ({exc}); leaving zeros") + person["tip_income"] = 0.0 + return person + + +def add_org_wages(person: pd.DataFrame, hh: pd.DataFrame, year: int, log) -> pd.DataFrame: + """Hourly wage / hourly-pay status / overtime from CPS-ORG donors. + + usdata's add_org_labor_market_inputs operates on an h5-like mapping of + arrays; a plain dict satisfies its read/write protocol. + """ + from policyengine_us_data.datasets.cps.cps import add_org_labor_market_inputs + + hh_state = pd.to_numeric(hh.get("state_fips", 0), errors="coerce").fillna(0) + + class _ZeroFallback(dict): + """h5-like mapping: unknown reads return zeros (logged once).""" + + def __init__(self, n, *a, **k): + super().__init__(*a, **k) + self._n = n + self._missed = set() + + def __getitem__(self, key): + if key in self: + return super().__getitem__(key) + if key not in self._missed: + self._missed.add(key) + return np.zeros(self._n, dtype=np.float32) + + n_persons = len(person) + cps = _ZeroFallback(n_persons) + cps.update({ + "age": pd.to_numeric(person.get("age", person.get("A_AGE", 0)), errors="coerce").fillna(0).to_numpy(np.float32), + "household_id": hh["household_id"].to_numpy(np.int64), + "person_household_id": person["person_household_id"].to_numpy(np.int64), + "state_fips": hh_state.to_numpy(np.float32), + "employment_income": pd.to_numeric(person["employment_income"], errors="coerce").fillna(0).to_numpy(np.float32), + "is_female": person.get("is_female", pd.Series(False, index=person.index)).astype(bool).to_numpy(), + "cps_race": pd.to_numeric(person.get("cps_race", 0), errors="coerce").fillna(0).to_numpy(np.float32), + "weekly_hours_worked": pd.to_numeric(person.get("hours_worked_last_week", 0), errors="coerce").fillna(0).to_numpy(np.float32), + "hours_worked_last_week": pd.to_numeric(person.get("hours_worked_last_week", 0), errors="coerce").fillna(0).to_numpy(np.float32), + "weeks_worked": pd.to_numeric(person.get("weeks_worked", 0), errors="coerce").fillna(0).to_numpy(np.float32), + "is_hispanic": person.get("is_hispanic", pd.Series(False, index=person.index)).astype(bool).to_numpy(), + }) + # Occupation flags the ORG models read — pass real pool values when present. + for flag in ("has_never_worked", "is_computer_scientist", + "is_executive_administrative_professional", + "is_farmer_fisher", "is_military"): + if flag in person.columns: + cps[flag] = person[flag].astype(bool).to_numpy() + try: + add_org_labor_market_inputs(cps, year) + if cps._missed: + log(f" ORG zero-fallback keys: {sorted(cps._missed)}") + person = person.copy() + for out in ("hourly_wage", "is_paid_hourly", "is_union_member_or_covered", + "weekly_hours_worked_before_lsr", "fsla_overtime_premium"): + if out in cps: + person[out] = np.asarray(cps[out]) + log(" ORG labor-market inputs imputed") + except Exception as exc: + log(f" ORG stage FAILED ({exc}); leaving defaults") + return person + + +def add_meps_esi_premiums(person: pd.DataFrame, log) -> pd.DataFrame: + """ESI premiums from MEPS-IC plan-type parameters (usdata rule, verbatim).""" + from policyengine_us_data.datasets.cps.cps import ( + impute_employer_sponsored_insurance_premiums, + ) + + person = person.copy() + person["employer_sponsored_insurance_premiums"] = ( + impute_employer_sponsored_insurance_premiums(person) + ) + nz = float((person["employer_sponsored_insurance_premiums"] > 0).mean()) + log(f" MEPS ESI premiums: nz {nz*100:.1f}%") + return person + + +def add_prior_year_income(person: pd.DataFrame, asec_year: int, log) -> pd.DataFrame: + """Prior-year earnings via the consecutive-ASEC PERIDNUM join (usdata rule). + + Maps last year's WSAL_VAL/SEMP_VAL onto matched persons; sentinel values + {-1, -9999} mean unavailable. + """ + from microplex.data_sources.cps import load_cps_asec + + person = person.copy() + if "PERIDNUM" not in person.columns: + log(" prior-year: PERIDNUM missing from pool; skipping") + return person + prior = load_cps_asec( + year=asec_year - 1, + extra_person_columns=["PERIDNUM", "WSAL_VAL", "SEMP_VAL"], + ).persons.to_pandas() + prior = prior.drop_duplicates("PERIDNUM").set_index("PERIDNUM") + sentinels = {-1, -9999} + cur_ids = person["PERIDNUM"] + emp = cur_ids.map(prior["WSAL_VAL"]) if "WSAL_VAL" in prior.columns else pd.Series(np.nan, index=person.index) + se = cur_ids.map(prior["SEMP_VAL"]) if "SEMP_VAL" in prior.columns else pd.Series(np.nan, index=person.index) + matched = emp.notna() & se.notna() & ~emp.isin(sentinels) & ~se.isin(sentinels) + person["employment_income_last_year"] = pd.to_numeric(emp, errors="coerce").where(matched, 0.0).fillna(0.0) + person["self_employment_income_last_year"] = pd.to_numeric(se, errors="coerce").where(matched, 0.0).fillna(0.0) + person["previous_year_income_available"] = matched.astype(bool) + log(f" prior-year join: matched {matched.mean()*100:.1f}% of persons (ASEC {asec_year-1})") + return person + + +def add_mortgage_conversion(person: pd.DataFrame, hh: pd.DataFrame, year: int, log) -> pd.DataFrame: + """Structural mortgages from SCF hints + PUF deductible interest. + + Ports usdata's two-step conversion (extended_cps.py:1183-1190): SCF-donor + balance hints, then conversion of the PUF-imputed interest deduction into + tax-unit mortgage balances/interest/origination plus person-level + home_mortgage_interest and the investment_interest_expense residual. + Failures raise — no silent zero fallbacks (charter rule). + """ + from policyengine_us_data.utils.mortgage_interest import ( + convert_mortgage_interest_to_structural_inputs, + impute_tax_unit_mortgage_balance_hints, + ) + + person = person.copy() + tp = year + p_tu = person["person_tax_unit_id"].to_numpy() + tu_ids = np.sort(np.unique(p_tu)) + tu_index = {t: i for i, t in enumerate(tu_ids.tolist())} + p_tu_idx = np.fromiter((tu_index[t] for t in p_tu.tolist()), dtype=np.int64) + + def person_col(name, default=0.0): + return pd.to_numeric(person.get(name, default), errors="coerce").fillna(0.0).to_numpy(np.float32) + + # Tax-unit interest deduction from the head-carried person values. + itd_person = person_col("interest_deduction") + itd_tu = np.zeros(len(tu_ids), dtype=np.float32) + np.add.at(itd_tu, p_tu_idx, itd_person) + + class _PersonZeroFallback(dict): + """Missing person-grain reads return person-length zeros (logged).""" + + def __init__(self, n): + super().__init__() + self._n = n + self.missed: set[str] = set() + + def get(self, key, default=None): + if key in self: + return super().__getitem__(key) + self.missed.add(key) + return {tp: np.zeros(self._n, dtype=np.float32)} + + def __getitem__(self, key): + if key in self: + return super().__getitem__(key) + self.missed.add(key) + return {tp: np.zeros(self._n, dtype=np.float32)} + + data = _PersonZeroFallback(len(person)) + data.update({ + "interest_deduction": {tp: itd_tu}, + "deductible_mortgage_interest": {tp: np.maximum(itd_person, 0)}, + "person_tax_unit_id": {tp: p_tu}, + "tax_unit_id": {tp: tu_ids}, + "person_id": {tp: person["person_id"].to_numpy()}, + "age": {tp: person_col("A_AGE" if "A_AGE" in person.columns else "age")}, + "employment_income": {tp: person_col("employment_income")}, + "self_employment_income": {tp: person_col("self_employment_income")}, + "social_security": {tp: person_col("social_security")}, + "taxable_pension_income": {tp: person_col("taxable_pension_income")}, + "taxable_interest_income": {tp: person_col("taxable_interest_income")}, + "is_female": {tp: person.get("is_female", pd.Series(False, index=person.index)).astype(bool).to_numpy()}, + "cps_race": {tp: person_col("cps_race")}, + "is_tax_unit_head": {tp: (person.get("tax_unit_role_input", "") == "HEAD").to_numpy()}, + "is_tax_unit_spouse": {tp: (person.get("tax_unit_role_input", "") == "SPOUSE").to_numpy()}, + "person_household_id": {tp: person["person_household_id"].to_numpy()}, + "person_spm_unit_id": {tp: person["person_spm_unit_id"].to_numpy()}, + "household_id": {tp: hh["household_id"].to_numpy()}, + "spm_unit_id": {tp: np.sort(person["person_spm_unit_id"].unique())}, + "tenure_type": {tp: hh.get("tenure_type", pd.Series("NONE", index=hh.index)).fillna("NONE").astype(str).to_numpy()}, + }) + # filing_status per tax unit: JOINT when a spouse is present, else SINGLE + # (the converter uses it only for the mortgage debt-cap split). + has_spouse_tu = np.zeros(len(tu_ids), dtype=bool) + np.add.at( + has_spouse_tu, + p_tu_idx, + (person.get("tax_unit_role_input", "") == "SPOUSE").to_numpy(), + ) + data["filing_status"] = {tp: np.where(has_spouse_tu, "JOINT", "SINGLE")} + # spm_unit_tenure_type: the SPM unit inherits its household's tenure. + hh_tenure = dict(zip(hh["household_id"], data["tenure_type"][tp])) + spm_ids = data["spm_unit_id"][tp] + spm_hh = ( + pd.DataFrame({"spm": person["person_spm_unit_id"], "hh": person["person_household_id"]}) + .drop_duplicates("spm").set_index("spm")["hh"] + ) + data["spm_unit_tenure_type"] = {tp: np.array([ + hh_tenure.get(spm_hh.get(sid), "NONE") for sid in spm_ids.tolist() + ])} + before = set(data) + data = impute_tax_unit_mortgage_balance_hints(data, tp) + data = convert_mortgage_interest_to_structural_inputs(data, tp) + + if data.missed: + log(f" mortgage converter zero-fallback keys: {sorted(data.missed)}") + tu_outputs = [] + for key in set(data) - before | {"interest_deduction"}: + arr = np.asarray(data[key][tp]) + if key.startswith("imputed_"): + continue + if len(arr) == len(person): + person[key] = arr + elif len(arr) == len(tu_ids): + # Head-carry tax-unit outputs onto persons; the export-time entity + # mover places them on the tax-unit table. + head = (person.get("tax_unit_role_input", "") == "HEAD").to_numpy() + vals = arr[p_tu_idx] + person[key] = np.where(head, vals, 0.0) + tu_outputs.append(key) + else: + raise ValueError(f"mortgage conversion output {key!r} has odd length {len(arr)}") + log(f" mortgage conversion: person outputs + head-carried {sorted(tu_outputs)}") + inv = person.get("investment_interest_expense") + if inv is not None: + log(f" investment_interest_expense: nz {(pd.to_numeric(inv, errors='coerce').fillna(0)>0).mean()*100:.1f}%") + return person + + +def add_acs_rent(person: pd.DataFrame, hh: pd.DataFrame, seed: int, log): + """Rent + vehicle ownership from the Census ACS 2022 artifact (usdata + storage), imputed at household grain with ACS household weights; + pre_subsidy_rent is person-entity and head-carried.""" + import h5py + from microimpute import Imputer + + acs_path = ( + "/Users/maxghenis/.claude-worktrees/usdata-populace/" + "policyengine_us_data/storage/acs_2022.h5" + ) + with h5py.File(acs_path) as f: + def col(name): + v = f[name][:] + return v + # rent is stored at person grain in the ACS artifact (head-carried); + # everything else here is household grain — aggregate before framing. + d_pers = pd.DataFrame({ + "person_household_id": col("person_household_id"), + "is_household_head": col("is_household_head").astype(bool), + "employment_income": col("employment_income"), + "rent": col("rent"), + }) + d_hh = pd.DataFrame({ + "household_id": col("household_id"), + "household_weight": col("household_weight"), + "state_fips": col("household_state_fips"), + "household_vehicles_owned": col("household_vehicles_owned"), + }) + g = d_pers.groupby("person_household_id") + d_hh = d_hh.merge( + pd.DataFrame({ + "hh_employment_income": g["employment_income"].sum(), + "hh_size": g.size(), + "rent": g["rent"].sum(), + }), + left_on="household_id", right_index=True, how="left", + ).fillna({"hh_employment_income": 0, "hh_size": 1, "rent": 0}) + d_hh = d_hh[d_hh["household_weight"] > 0] + + PRED = ["state_fips", "hh_employment_income", "hh_size"] + fitted = Imputer(seed=seed, log_level="WARNING").fit( + d_hh.dropna(subset=PRED + ["rent"]), + PRED, ["rent"], + weight_col="household_weight", + ) + pg = person.groupby(person["person_household_id"]) + recv = pd.DataFrame({ + "hh_employment_income": pd.to_numeric(person["employment_income"], errors="coerce").fillna(0).groupby(person["person_household_id"]).sum(), + "hh_size": pg.size(), + }) + recv = recv.reindex(hh["household_id"]).fillna(0) + recv["state_fips"] = pd.to_numeric(hh.get("state_fips", 0), errors="coerce").fillna(0).to_numpy() + draws = fitted.predict(recv[PRED].reset_index(drop=True)) + + hh = hh.copy() + rent = _support_guard(np.asarray(draws["rent"], dtype=np.float64), d_hh["rent"].to_numpy(np.float64), "rent", log) + # Rent applies to renter households only (tenure from the CPS H_TENURE map). + tenure = hh.get("tenure_type", pd.Series("NONE", index=hh.index)).astype(str) + rent = np.where(tenure.str.upper() == "RENTED", rent, 0.0) + # pre_subsidy_rent is person-entity: head-carry the household rent. + if "is_household_head" not in person.columns: + raise RuntimeError( + "head-carry requires is_household_head on the person frame " + "(derive it from ASEC P_SEQ == 1 before stage F2); refusing to " + "head-carry onto an all-False mask." + ) + headp = person["is_household_head"].astype(bool) + hmap = dict(zip(hh["household_id"].tolist(), range(len(hh)))) + pidx = person["person_household_id"].map(hmap) + person = person.copy() + person["pre_subsidy_rent"] = np.where( + headp & pidx.notna(), rent[pidx.fillna(0).astype(int)], 0.0 + ) + log(f" ACS rent: renter-hh rent nz {(rent>0).mean()*100:.1f}%") + return person, hh + + +def add_vehicle_assets(person: pd.DataFrame, hh: pd.DataFrame, log): + """Household vehicles (count + value) from the SIPP-trained QRF donor + (usdata get_vehicle_model + receiver builder), mirroring usdata's + auto-loan/vehicle imputation. Writes household grain.""" + from policyengine_us_data.datasets.sipp import get_vehicle_model + from policyengine_us_data.utils.asset_imputation import ( + build_household_vehicle_receiver, + ) + + model = get_vehicle_model() + # Build the receiver from a controlled column set: the person frame may + # already carry a household_id-named column, and a duplicate name breaks + # the builder's groupby. + receiver_cols = [ + c + for c in ( + "employment_income", + "interest_income", + "dividend_income", + "interest_dividend_income", + "rental_income", + "age", + "is_female", + "is_married", + "is_household_head", + ) + if c in person.columns + ] + receiver_person = person[receiver_cols].copy() + receiver_person["household_id"] = person["person_household_id"].to_numpy() + tenure = hh.get("tenure_type") + receiver = build_household_vehicle_receiver( + receiver_person, + tenure_type=(np.asarray(tenure) if tenure is not None else None), + ) + pred = model.predict(X_test=receiver, mean_quantile=0.5) + owned = np.clip( + np.rint(np.asarray(pred["household_vehicles_owned"], dtype=np.float64)), + 0, + None, + ) + value = np.clip( + np.asarray(pred["household_vehicles_value"], dtype=np.float64), 0, None + ) + # receiver rows are one per household in hh order (builder groups by + # household_id of the persons); align defensively by id. + rid = np.asarray(receiver["household_id"]) if "household_id" in receiver else None + hh = hh.copy() + if rid is not None: + owned_by = dict(zip(rid.tolist(), owned.tolist())) + value_by = dict(zip(rid.tolist(), value.tolist())) + hh["household_vehicles_owned"] = ( + hh["household_id"].map(owned_by).fillna(0.0).astype(float) + ) + hh["household_vehicles_value"] = ( + hh["household_id"].map(value_by).fillna(0.0).astype(float) + ) + else: + if len(owned) != len(hh): + raise RuntimeError( + f"vehicle receiver rows ({len(owned)}) != households " + f"({len(hh)}) and no household_id to align by." + ) + hh["household_vehicles_owned"] = owned + hh["household_vehicles_value"] = value + # usdata folds the vehicle value into net worth (cps.py net_worth + # components assembly); mirror that here, where the value is imputed. + if "net_worth" in hh.columns: + hh["net_worth"] = ( + hh["net_worth"].to_numpy(dtype=np.float64) + + hh["household_vehicles_value"].to_numpy(dtype=np.float64) + ) + log( + f" SIPP vehicles: owned nz {(hh['household_vehicles_owned']>0).mean()*100:.1f}%, " + f"value nz {(hh['household_vehicles_value']>0).mean()*100:.1f}%, " + f"folded into net_worth" + ) + return person, hh diff --git a/packages/populace-data/build/us/release_manifest.json b/packages/populace-data/build/us/release_manifest.json new file mode 100644 index 00000000..cd824232 --- /dev/null +++ b/packages/populace-data/build/us/release_manifest.json @@ -0,0 +1,58 @@ +{ + "build_id": "populace-us-2024-9f1260b-20260611", + "builder": "populace", + "build_sha": "9f1260b", + "build_date": "2026-06-11", + "dataset": { + "filename": "populace_us_2024.h5", + "sha256": "dc75c0d4fdedd57946db84a8d838dbc5b61a284365c3ce6eb6508b8e81111a4b" + }, + "calibration": { + "filename": "populace_us_2024_calibration.npz", + "sha256": "a3da2f59085c45f0e16b06337818e3513c2635911dc0d16fa7deb5006263c12a" + }, + "construction": "eCPS-free: every layer from primary sources (CPS ASEC, IRS PUF 2015 uprated, Fed SCF 2022, Census SIPP, CPS-ORG, MEPS-IC parameters, Census ACS 2022); enhanced CPS used only as the scoring benchmark", + "gates": { + "parity_gaps": 0, + "exported_nonzero": { + "passed": true, + "stored_columns": 309 + }, + "calibration": { + "within_10pct_share": 0.9509, + "loss": 0.022, + "max_weight": 297651, + "weights_above_500k": 0, + "max_weight_ratio": 50 + }, + "smoke": { + "people_m": 332.3, + "snap_b": 97.8, + "net_worth_t": 176.5, + "net_stcg_b": -77.5, + "tips_b": 52.9, + "pre_subsidy_rent_b": 757.5, + "investment_interest_expense_b": 5.1 + } + }, + "score_vs_enhanced_cps": { + "protocol": "matched 41,314 households, symmetric refit, 739-target holdout (seed 20260529)", + "train_loss": { + "populace": 0.18957, + "enhanced_cps": 1.08879 + }, + "holdout_loss": { + "populace": 0.03837, + "enhanced_cps": 0.3167 + }, + "full_loss": { + "populace": 0.22794, + "enhanced_cps": 1.40549 + }, + "per_target_wins": { + "populace": 1040, + "enhanced_cps": 2613, + "ties": 51 + } + } +} \ No newline at end of file diff --git a/packages/populace-data/build/us/run_chain.sh b/packages/populace-data/build/us/run_chain.sh new file mode 100755 index 00000000..cff852e2 --- /dev/null +++ b/packages/populace-data/build/us/run_chain.sh @@ -0,0 +1,30 @@ +#!/bin/zsh +# The v3 build chain, end to end. Each step's venv is the one that owns its +# dependencies (worktree .venv: engine + usdata; populace venv: calibrate + +# gates). Fails fast; the log narrates which step died. +set -e +WT=~/.claude-worktrees/microplex-spec-build +# One chain = one run id = one fresh log; never inherit a prior chain's. +rm -f /tmp/populace_run_id +echo $$ > /tmp/populace_chain.pid +trap 'rm -f /tmp/populace_chain.pid' EXIT +BUILD_PY=$WT/.venv/bin/python +POP_PY=/tmp/populace-build-venv/bin/python + +echo "=== CHAIN step 1: full build ===" +$BUILD_PY -u $WT/scripts/build_us_candidate.py --mode full \ + --usdata-repo ~/.claude-worktrees/usdata-populace + +echo "=== CHAIN step 2: extract target surface ===" +$BUILD_PY -u $WT/scripts/extract_target_surface.py + +echo "=== CHAIN step 3: calibrate + artifact ===" +$POP_PY -u $WT/scripts/build_dataset.py + +echo "=== CHAIN step 4: enrich (sim-dependent layers) ===" +$BUILD_PY -u $WT/scripts/enrich_artifact.py + +echo "=== CHAIN step 5: gates (exported-nonzero + parity + smoke) ===" +$POP_PY -u $WT/scripts/check_parity.py + +echo "=== CHAIN COMPLETE: all gates green ==="