diff --git a/.gitignore b/.gitignore index a815b5ec..cd4e0f41 100644 --- a/.gitignore +++ b/.gitignore @@ -46,3 +46,6 @@ docs/book/_build/* oguk/data/ uv.lock test_api.py + +# ogcore demographics drops this in the CWD +un_api_token.txt diff --git a/oguk/api.py b/oguk/api.py index 0ec2e71c..4b84491a 100644 --- a/oguk/api.py +++ b/oguk/api.py @@ -288,13 +288,91 @@ class Config: year: int +def _resolve_year_dataset(datasets: dict, year: int): + """Pick the per-year dataset from an ensure_datasets result, across + policyengine-uk vintages. + + policyengine-uk >= 2.89 keys per-year datasets as + ``{stem}_{year}`` with a ``populace_uk_*`` default stem; 2.88 used + ``enhanced_frs_2023_24_{year}``. Matching on the year suffix instead of + hardcoding the stem keeps OG-UK working across the rename (and the next + one); when a vintage ships several stems for the same year (2.88 + returned both frs_* and enhanced_frs_*), the calibrated stem is + preferred. See PSLmodels/OG-UK#68. + """ + matches = sorted(k for k in datasets if k.endswith(f"_{year}")) + if not matches: + raise KeyError( + f"no dataset for year {year} in ensure_datasets result; " + f"available keys: {sorted(datasets)}" + ) + if len(matches) == 1: + return datasets[matches[0]] + # Multiple stems for the year: policyengine-uk 2.88's default returned + # BOTH frs_2023_24_* and enhanced_frs_2023_24_* — prefer the calibrated + # stems in order, and only refuse when preference cannot break the tie. + for stem in ("populace_uk", "enhanced_frs"): + preferred = [k for k in matches if k.startswith(stem)] + if len(preferred) == 1: + return datasets[preferred[0]] + if len(preferred) > 1: + break + raise KeyError( + f"ambiguous datasets for year {year}: {matches}; pass a single " + "dataset stem to ensure_datasets" + ) + + +def _perturb_first_populated(s, candidates, delta): + """Add ``delta`` to the first POPULATED input among ``candidates``. + + Robustness across policyengine-uk data layouts: employment income + may be carried by employment_income_before_lsr (populated at + simulation construction, with employment_income derived from it) + or by employment_income directly. Perturbing whichever holder + actually holds values — and raising when none does — guarantees + the finite difference is never silently flat. + """ + for name in candidates: + holder = s.get_holder(name) + periods = list(holder.get_known_periods()) + if not periods: + continue + for period in periods: + values = holder.get_array(period) + holder.delete_arrays(period) + s.set_input(name, period, values + delta) + return s + raise RuntimeError( + f"no populated input among {candidates}; cannot apply the " + "MTR perturbation — refusing to return flat marginal rates" + ) + + +def _build_perturbation_modifier(reform_modifier, year, candidates): + """Compose the reform (first) with a GBP 1 per-adult income perturbation. + + Ordering is load-bearing: the perturbation must land on the REFORMED + world so the finite difference measures marginal rates under the + policy being scored, not under the baseline. + """ + + def _modify(s): + if reform_modifier is not None: + s = reform_modifier(s) or s + adult = s.calculate("is_adult", year) + return _perturb_first_populated(s, candidates, adult) + + return _modify + + def _get_micro_data(year: int, policy: Policy | None, data_folder: str) -> _MicroData: """Extract microdata from PolicyEngine-UK (internal).""" from policyengine.core import Policy, Simulation from policyengine.tax_benefit_models.uk import ensure_datasets, uk_latest datasets = ensure_datasets(data_folder=data_folder, years=[year]) - dataset = datasets[f"enhanced_frs_2023_24_{year}"] + dataset = _resolve_year_dataset(datasets, year) sim = Simulation( dataset=dataset, tax_benefit_model_version=uk_latest, policy=policy @@ -326,61 +404,125 @@ def _hbai_net_per_person(microsim): baseline_net = _hbai_net_per_person(sim) - # Build perturbation policies. PolicyEngine silently drops parameter_values - # when a simulation_modifier is present, so we must apply any reform - # parameter changes *inside* the modifier using the TBS parameter tree. + # Household context for the finite-difference MTRs below. The GBP 1 + # perturbation is applied to EVERY adult at once (one extra simulation, + # not one per person), so person-level rates must be recovered from + # person-level tax deltas plus a per-adult share of the household-level + # residual (benefit withdrawal). Averaging the whole household delta is + # NOT enough: it assigns each adult the household-average rate (a + # GBP 30k / GBP 0 couple would both read ~0.14 instead of 0.28 / 0) — + # and the previous code did not divide at all, so 1 - delta went + # negative for every multi-adult household and np.clip floored it to + # ZERO (measured: 80% of GBP 25-35k earners at a 0.000 labour MTR). + is_adult_arr = (age >= 18).astype(float) + hh_ids = person["household_id"].values + adults_per_hh = ( + pd.DataFrame({"hh": hh_ids, "adult": is_adult_arr}) + .groupby("hh")["adult"] + .transform("sum") + .values + ) + adults_per_hh = np.maximum(adults_per_hh, 1.0) + + def _person_tax(s): + p = s.output_dataset.data.person + n = len(p) + return ( + p.get("income_tax", np.zeros(n)).values + + p.get("national_insurance", np.zeros(n)).values + ) + + base_tax = _person_tax(sim) + + def _person_mtr(pert_sim): + """Per-person MTR from a simultaneously perturbed simulation. + + Identity per household (GBP 1 to each of its n adults): + + dnet_hh = n - sum_adults(dtax_i) + resid_hh + + where resid_hh is the benefit response net of any NON-perturbed + members' tax response (e.g. a tax-paying child under an HICBC + shift). + + This is an explicit INCIDENCE APPROXIMATION, not an exact + per-person finite difference (which would need one simulation per + person). dtax_i is person i's own tax delta under the simultaneous + perturbation; the household residual — benefit withdrawal, + non-perturbed members' tax responses, and any other household-level + effect — is shared equally across the perturbed adults: + + mtr_i = dtax_i - resid_hh / n + + Where taxes couple ACROSS perturbed adults (marriage allowance + transfers, HICBC), the coupled component is misattributed between + the partners. Before clipping, the adult rates sum to the exact + household withdrawal; partner-level and conditional distributions + (quantiles, variance) remain approximate, and the clip below can + move even the household sum. Genuine negative marginal + rates (benefit phase-ins) are clipped to 0 by the long-standing + [0, 1] clip. Policies carrying a custom simulation_modifier are + not composed into the perturbation runs (pre-existing behaviour); + parametric reforms are, via the engine helper below. + """ + dtax = _person_tax(pert_sim) - base_tax + dnet_hh = _hbai_net_per_person(pert_sim) - baseline_net + # Only the PERTURBED (adult) members' own tax deltas enter the + # identity term; a non-perturbed member whose tax moved (e.g. a + # tax-paying child under an HICBC shift) is a cross-person response + # and belongs in the shared residual, not in anyone's own rate. + sum_adult_dtax_hh = ( + pd.DataFrame({"hh": hh_ids, "dtax": dtax * is_adult_arr}) + .groupby("hh")["dtax"] + .transform("sum") + .values + ) + resid_hh = dnet_hh - (adults_per_hh - sum_adult_dtax_hh) + return np.clip(dtax - resid_hh / adults_per_hh, 0, 1) + + # Reform parameter replay: the engine drops parameter_values when a + # simulation_modifier is present, so the reform must be applied inside + # the modifier — via the engine's OWN helper, which honours each + # ParameterValue's full dating. (A hand-rolled replay here applied a + # single year only, so on multi-year windows the baseline kept the + # reform while the perturbation runs lost it.) + from policyengine.utils.parametric_reforms import ( + simulation_modifier_from_parameter_values, + ) + reform_param_values = policy.parameter_values if policy else [] + reform_modifier = ( + simulation_modifier_from_parameter_values(reform_param_values) + if reform_param_values + else None + ) - def _apply_reform_params(microsim): - """Replay reform parameter_values on the internal TBS.""" - for pv in reform_param_values: - param_path = ( - pv.parameter.name - ) # e.g. "gov.hmrc.income_tax.rates.uk[0].rate" - node = microsim.tax_benefit_system.parameters - for part in param_path.split("."): - # Handle bracket indexing like "uk[0]" - if "[" in part: - name, idx = part.split("[") - idx = int(idx.rstrip("]")) - node = getattr(node, name).brackets[idx] - else: - node = getattr(node, part) - start = pv.start_date - period = f"year:{start.year}:1" - node.update(period=period, value=pv.value) - - # Labour MTR: perturb employment income by £1 - def add_labor(s): - _apply_reform_params(s) - emp = s.calculate("employment_income", year) - adult = s.calculate("is_adult", year) - s.set_input("employment_income", year, emp + adult) - return s + # Labour MTR: perturb employment income by GBP 1 per adult + add_labor = _build_perturbation_modifier( + reform_modifier, year, ("employment_income_before_lsr", "employment_income") + ) - labor_pol = Policy(name="labor_perturb", simulation_modifier=add_labor) + # Distinct policy names per run: purely for traceability in logs and + # cached-artifact listings (simulation ids are already unique). + import uuid as _uuid + + run_tag = _uuid.uuid4().hex[:10] + labor_pol = Policy(name=f"labor_perturb_{run_tag}", simulation_modifier=add_labor) labor_sim = Simulation( dataset=dataset, tax_benefit_model_version=uk_latest, policy=labor_pol ) labor_sim.ensure() - labor_net = _hbai_net_per_person(labor_sim) - mtr_labor = np.clip(1 - (labor_net - baseline_net), 0, 1) + mtr_labor = _person_mtr(labor_sim) - # Capital MTR: perturb dividend income by £1 - def add_cap(s): - _apply_reform_params(s) - div = s.calculate("dividend_income", year) - adult = s.calculate("is_adult", year) - s.set_input("dividend_income", year, div + adult) - return s + # Capital MTR: perturb dividend income by GBP 1 per adult + add_cap = _build_perturbation_modifier(reform_modifier, year, ("dividend_income",)) - cap_pol = Policy(name="cap_perturb", simulation_modifier=add_cap) + cap_pol = Policy(name=f"cap_perturb_{run_tag}", simulation_modifier=add_cap) cap_sim = Simulation( dataset=dataset, tax_benefit_model_version=uk_latest, policy=cap_pol ) cap_sim.ensure() - cap_net = _hbai_net_per_person(cap_sim) - mtr_capital = np.clip(1 - (cap_net - baseline_net), 0, 1) + mtr_capital = _person_mtr(cap_sim) market_inc = labor_inc + cap_inc tax = ( diff --git a/oguk/tests/test_get_micro_data.py b/oguk/tests/test_get_micro_data.py index 788f2dd9..16476bfa 100644 --- a/oguk/tests/test_get_micro_data.py +++ b/oguk/tests/test_get_micro_data.py @@ -1,14 +1,60 @@ """Tests for OG-UK calibration API.""" +import os from datetime import datetime +import pytest from policyengine.core import ParameterValue, Policy from policyengine.tax_benefit_models.uk import uk_latest from oguk import CalibrationResult, calibrate -def test_baseline_calibration(): +# The UK microdata lives in a private Hugging Face repo. Same-repo CI has the +# token secret; FORK pull requests do not (GitHub withholds secrets), and +# contributors may not have access either — so data-dependent tests skip +# cleanly without a token instead of failing on a 401. +def _hf_token_present() -> bool: + """True when any Hugging Face credential is available — env vars or the + hub's stored login (huggingface-cli login writes a token file that the + hub uses regardless of the environment).""" + if os.environ.get("HUGGING_FACE_TOKEN") or os.environ.get("HF_TOKEN"): + return True + try: + from huggingface_hub import get_token + + return bool(get_token()) + except Exception: + return False + + +requires_uk_microdata = pytest.mark.skipif( + not _hf_token_present(), + reason="needs a Hugging Face token with access to the private UK microdata", +) + + +@pytest.fixture() +def _headless_un_token(): + """ogcore.demographics.get_un_data prompts on stdin for a UN API token + unless un_api_token.txt exists in the CWD; under pytest the prompt + raises (stdin is captured, and ogcore catches only EOFError). Pre-seed + an empty token: the UN API then returns 401 and ogcore falls back to + its public GitHub population-data mirror, which needs no auth. The + file is removed only if this fixture created it.""" + import pathlib + + path = pathlib.Path("un_api_token.txt") + created = not path.exists() + if created: + path.write_text("") + yield + if created: + path.unlink(missing_ok=True) + + +@requires_uk_microdata +def test_baseline_calibration(_headless_un_token): """Test baseline calibration produces valid results.""" result = calibrate(start_year=2026, years=1) @@ -18,7 +64,8 @@ def test_baseline_calibration(): assert len(result.omega_SS) > 0 -def test_reform_calibration(): +@requires_uk_microdata +def test_reform_calibration(_headless_un_token): """Test calibration with a policy reform.""" pa_param = uk_latest.get_parameter( "gov.hmrc.income_tax.allowances.personal_allowance.amount" @@ -34,13 +81,16 @@ def test_reform_calibration(): ], ) - result = calibrate(start_year=2026, years=1, policy=reform) + result = calibrate(start_year=2026, years=2, policy=reform) assert isinstance(result, CalibrationResult) assert result.mean_income > 0 + # multi-year: one fitted ETR set per year, both under the reform + assert len(result.etr_params) == 2 -def test_demographic_outputs(): +@requires_uk_microdata +def test_demographic_outputs(_headless_un_token): """Test demographic parameters are valid.""" result = calibrate(start_year=2026, years=1) @@ -53,3 +103,227 @@ def test_demographic_outputs(): # Population shares should sum to 1 assert abs(result.omega_SS.sum() - 1.0) < 0.01 + + +# --- dataset-key resolution across policyengine-uk vintages (issue #68) --- + + +def test_resolve_year_dataset_populace_keys(): + from oguk.api import _resolve_year_dataset + + ds = {"populace_uk_2023_2026": "a", "populace_uk_2023_2027": "b"} + assert _resolve_year_dataset(ds, 2026) == "a" + assert _resolve_year_dataset(ds, 2027) == "b" + + +def test_resolve_year_dataset_legacy_enhanced_frs_keys(): + from oguk.api import _resolve_year_dataset + + ds = {"enhanced_frs_2023_24_2026": "x"} + assert _resolve_year_dataset(ds, 2026) == "x" + + +def test_resolve_year_dataset_missing_year_lists_available(): + import pytest + + from oguk.api import _resolve_year_dataset + + with pytest.raises(KeyError, match="available keys"): + _resolve_year_dataset({"populace_uk_2023_2027": "b"}, 2026) + + +def test_resolve_year_dataset_prefers_calibrated_stems(): + """2.88's default ensure_datasets returned BOTH frs_* and + enhanced_frs_* for each year — the calibrated stem must win, and + populace_uk_* outranks both.""" + from oguk.api import _resolve_year_dataset + + legacy_pair = {"frs_2023_24_2026": "raw", "enhanced_frs_2023_24_2026": "x"} + assert _resolve_year_dataset(legacy_pair, 2026) == "x" + mixed = {"populace_uk_2023_2026": "a", "enhanced_frs_2023_24_2026": "x"} + assert _resolve_year_dataset(mixed, 2026) == "a" + + +def test_resolve_year_dataset_ambiguous_within_stem_refused(): + from oguk.api import _resolve_year_dataset + + with pytest.raises(KeyError, match="ambiguous"): + _resolve_year_dataset( + {"populace_uk_2023_2026": "a", "populace_uk_2024_2026": "b"}, 2026 + ) + + +def _extract_or_skip(year, policy=None): + """Run the microdata extraction, skipping (not failing) when the private + UK dataset is unreachable from this environment. Attempt-based, because + availability is more than env vars: a huggingface-cli login or a warm + hub cache also works.""" + import tempfile + + from oguk.api import _get_micro_data + + try: + with tempfile.TemporaryDirectory() as tmp: + return _get_micro_data(year, policy, tmp) + except Exception as e: # noqa: BLE001 — availability gate + marker = f"{type(e).__name__}: {e}" + if any( + s in marker + for s in ( + "401", + "Unauthorized", + "RepositoryNotFound", + "GatedRepo", + "LocalEntryNotFound", + ) + ): + pytest.skip(f"UK microdata unavailable here: {marker[:200]}") + raise + + +def test_person_level_mtrs_discriminate_household_structure(): + """The two historical failure modes, pinned separately: + + (1) no household normalisation -> multi-adult households clip to a + 0.000 MTR (80% of mid-band earners measured); + (2) household-AVERAGING -> every adult gets the household mean (a + GBP 30k / GBP 0 couple both read ~0.14). + + Person-level rates must show mid-band earners near the statutory + ~28% in BOTH single- and multi-adult households, while zero-earning + partners in those same households stay near zero. + """ + import numpy as np + + md = _extract_or_skip(2026) + # _MicroData has no household ids, so rebuild bands from income alone: + mid = (md.labor_income > 25_000) & (md.labor_income < 35_000) + med_mid = float(np.median(md.mtr_labor[mid])) + assert 0.20 < med_mid < 0.45, f"mid-band median {med_mid:.3f}" + # THE household-structure discriminator is the mid-band LOWER TAIL: + # under household-averaging, a GBP 30k earner married to a non-earner + # reads ~0.14, dragging q10 to ~0.14 (measured); person-level rates + # put every mid-band earner near the statutory ~0.28 (q10 measured + # 0.277). Under the original no-normalisation clipping, the median + # itself was 0.000. + q10_mid = float(np.percentile(md.mtr_labor[mid], 10)) + assert q10_mid > 0.20, ( + f"mid-band q10 {q10_mid:.3f} — household averaging is back " + "(earners in multi-adult households diluted toward the mean)" + ) + # Low-EARNINGS adults are not low-INCOME (pensioners' marginal rate on + # GBP 1 of earnings is legitimately ~20%), so no near-zero assertion + # there — but averaging would also DILUTE their rates toward household + # means; sanity-bound the population mean instead. + assert 0.15 < float(md.mtr_labor.mean()) < 0.40 + # Capital: the GBP 1 dividend perturbation sits inside the GBP 500 + # dividend allowance for anyone without existing dividends, so a zero + # MEDIAN among broad capital-income holders (mostly pension income) is + # correct; require instead that a real mass faces positive dividend + # marginal rates (8.75%/33.75% bands). + divs = md.capital_income > 2_000 + assert divs.sum() > 1_000, ( + f"only {int(divs.sum())} capital-income holders in the sample — " + "the capital assertions below would be vacuous" + ) + positive_share = float(np.mean(md.mtr_capital[divs] > 0.05)) + assert positive_share > 0.10, ( + f"only {positive_share:.0%} of capital-income holders face a " + "positive marginal dividend rate — perturbation likely flat" + ) + assert 0.02 < float(md.mtr_capital.mean()) < 0.30 + + +def test_labor_mtr_bites_for_midband_earners(): + """Regression for silently flat labour MTRs. The originally observed + failure (~80% of £25-35k earners at a 0.000 MTR) was clipping caused + by comparing household-level net-income deltas against per-person tax + deltas without household normalisation — NOT an input-layout change; + the perturbation reaches the tax pipeline on both 2.88 and 2.89 + layouts via _perturb_first_populated. A mid-band earner's marginal + rate is ~28% (basic-rate income tax + employee NI); the median across + the band must land near it.""" + import numpy as np + + md = _extract_or_skip(2026) + mask = (md.labor_income > 25_000) & (md.labor_income < 35_000) + median_mtr = float(np.median(md.mtr_labor[mask])) + assert 0.20 < median_mtr < 0.45, ( + f"mid-band labour MTR median {median_mtr:.3f} — the perturbation " + "is not reaching the tax pipeline" + ) + + +def test_reform_modifier_composes_before_perturbation(): + """Parametric reforms are applied INSIDE the perturbation modifier, + before the income perturbation — pinned with fakes so the ordering + cannot silently invert (perturbing the baseline world instead of the + reformed one).""" + import numpy as np + + from oguk.api import _build_perturbation_modifier + + calls = [] + + class _FakeHolder: + def __init__(self): + self.arrays = {2026: np.array([1.0, 2.0])} + + def get_known_periods(self): + return list(self.arrays) + + def get_array(self, period): + return self.arrays[period] + + def delete_arrays(self, period): + del self.arrays[period] + + class _FakeSim: + def __init__(self): + self.holder = _FakeHolder() + self.values = None + + def get_holder(self, name): + return self.holder + + def calculate(self, name, year): + calls.append(f"calc:{name}") + return np.array([1.0, 0.0]) # adult, child + + def set_input(self, name, period, values): + calls.append(f"perturb:{name}") + self.values = values + + def fake_reform(s): + calls.append("reform") + return s + + modifier = _build_perturbation_modifier(fake_reform, 2026, ("employment_income",)) + sim = _FakeSim() + modifier(sim) + assert calls[0] == "reform", calls + assert "perturb:employment_income" in calls, calls + # GBP 1 added only for the adult member + assert list(sim.values) == [2.0, 2.0] + + +def test_perturbation_refuses_silently_flat(): + """No populated candidate input -> hard error, never flat MTRs.""" + import numpy as np + + from oguk.api import _build_perturbation_modifier + + class _EmptyHolder: + def get_known_periods(self): + return [] + + class _FakeSim: + def get_holder(self, name): + return _EmptyHolder() + + def calculate(self, name, year): + return np.array([1.0]) + + modifier = _build_perturbation_modifier(None, 2026, ("employment_income",)) + with pytest.raises(RuntimeError, match="refusing"): + modifier(_FakeSim()) diff --git a/pyproject.toml b/pyproject.toml index 088991bd..438e34de 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,9 +28,13 @@ dependencies = [ "pandas>=2.0", "matplotlib>=3.7", "pydantic>=2.0", - "ogcore>=0.15.6", - "policyengine>=4.3", - "policyengine-uk==2.88.0", + # ogcore 0.18 makes income_percentiles effectively required in + # get_pop_objs (expand_pop_obj_J asserts before its homogeneous-case + # branch), breaking calibrate(); cap until oguk migrates to the + # income-heterogeneity API. + "ogcore>=0.15.6,<0.18", + "policyengine>=4.22", + "policyengine-uk>=2.89.2", "requests>=2.31", "rich>=13.0", "openpyxl>=3.1",