From c4cddcdc9746dd2aa74c790ad70685ecb10833eb Mon Sep 17 00:00:00 2001 From: Maria Juaristi Date: Tue, 20 Oct 2026 14:00:41 +0100 Subject: [PATCH 1/4] Add the #1326 tool: 2023-24 Scottish property income by NSND band A project tool for one external request, not for main. It reads a released UK national file and its build's national target registry and writes private outputs outside the repository: - a one-year back-cast that inverts policyengine-uk's own forward uprating, with a guard against the engine's silent zeros for earlier periods; - a bounded reweight to the SPI 2023-24 facts the build binds (UK property income counts and amounts by band; Scottish taxpayers and total income by band), holding each region's households; - Release 1 (the ten requested bands) and Release 2 (four bins aligned to the 2023-24 Scottish rate thresholds), with distinct records, effective records and source-bootstrap intervals per band, and restricted and public versions with hierarchical suppression; - an illustrative 1p on the Scottish rates charged on property income, with and without the finance-cost reduction following. Co-Authored-By: Claude Opus 5.5 --- ...test_uk_scotland_property_income_engine.py | 135 ++++++ .../uk/test_uk_scotland_property_income.py | 283 ++++++++++++ tools/uk_scotland_property_income/__init__.py | 17 + tools/uk_scotland_property_income/__main__.py | 424 ++++++++++++++++++ tools/uk_scotland_property_income/backcast.py | 86 ++++ tools/uk_scotland_property_income/bands.py | 79 ++++ tools/uk_scotland_property_income/measures.py | 106 +++++ tools/uk_scotland_property_income/rake.py | 110 +++++ tools/uk_scotland_property_income/reform.py | 82 ++++ tools/uk_scotland_property_income/registry.py | 101 +++++ tools/uk_scotland_property_income/release.py | 363 +++++++++++++++ 11 files changed, 1786 insertions(+) create mode 100644 packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py create mode 100644 packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py create mode 100644 tools/uk_scotland_property_income/__init__.py create mode 100644 tools/uk_scotland_property_income/__main__.py create mode 100644 tools/uk_scotland_property_income/backcast.py create mode 100644 tools/uk_scotland_property_income/bands.py create mode 100644 tools/uk_scotland_property_income/measures.py create mode 100644 tools/uk_scotland_property_income/rake.py create mode 100644 tools/uk_scotland_property_income/reform.py create mode 100644 tools/uk_scotland_property_income/registry.py create mode 100644 tools/uk_scotland_property_income/release.py diff --git a/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py b/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py new file mode 100644 index 000000000..acfa33a1e --- /dev/null +++ b/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py @@ -0,0 +1,135 @@ +"""The microcosm#1326 tool against policyengine-uk: the back-cast inverts the +engine's own forward step, and the Scottish property-rate reform charges +exactly ``delta`` on the property slice.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest +from policyengine_uk import CountryTaxBenefitSystem, Microsimulation +from policyengine_uk.data import UKSingleYearDataset +from policyengine_uk.data.economic_assumptions import extend_single_year_dataset + +from tools.uk_scotland_property_income.backcast import ( + backcast_one_year, + refuse_earlier_period, + with_household_weights, +) +from tools.uk_scotland_property_income.reform import scottish_property_rate_reform + +REGIONS = ("SCOTLAND", "SCOTLAND", "LONDON") + + +def _dataset(year: int) -> UKSingleYearDataset: + ids = np.arange(len(REGIONS), dtype=np.int64) + return UKSingleYearDataset( + person=pd.DataFrame( + { + "person_id": ids, + "person_benunit_id": ids, + "person_household_id": ids, + "age": [45, 52, 45], + "employment_income": [30_000.0, 48_000.0, 30_000.0], + "property_income": [10_000.0, 14_000.0, 10_000.0], + "property_finance_costs": [0.0, 3_000.0, 0.0], + "savings_interest_income": [200.0, 0.0, 200.0], + } + ), + benunit=pd.DataFrame({"benunit_id": ids}), + household=pd.DataFrame( + { + "household_id": ids, + "household_weight": [1_000.0, 800.0, 1_000.0], + "region": list(REGIONS), + "council_tax": [1_500.0, 1_800.0, 1_500.0], + "rent": 0.0, + "tenure_type": "OWNED_OUTRIGHT", + } + ), + fiscal_year=year, + ) + + +def test_the_backcast_inverts_the_engines_forward_step(): + parameters = CountryTaxBenefitSystem().parameters + released = _dataset(2024) + earlier, divided = backcast_one_year(released, parameters) + assert earlier.time_period == "2023" + assert { + "employment_income", + "property_income", + "property_finance_costs", + "household_weight", + } <= set(divided) + assert ( + earlier.person["employment_income"].iloc[0] + < released.person["employment_income"].iloc[0] + ) + + forward = extend_single_year_dataset(earlier, parameters, end_year=2024).datasets[ + 2024 + ] + for table in ("person", "household"): + for column in getattr(released, table).columns: + if column in divided: + np.testing.assert_allclose( + getattr(forward, table)[column].to_numpy(dtype=float), + getattr(released, table)[column].to_numpy(dtype=float), + rtol=1e-12, + err_msg=column, + ) + + +def test_an_earlier_period_is_refused(): + with pytest.raises(ValueError, match="Back-cast it first"): + refuse_earlier_period(_dataset(2024), 2023) + refuse_earlier_period(_dataset(2024), 2024) + + +def _income_tax(dataset: UKSingleYearDataset, reform=None) -> np.ndarray: + simulation = Microsimulation( + dataset=with_household_weights(dataset, dataset.household["household_weight"]), + reform=reform, + ) + return np.asarray(simulation.calculate("income_tax", 2023), dtype=float), simulation + + +def test_the_reform_adds_delta_on_the_scottish_property_slice_only(): + dataset = _dataset(2023) + baseline, simulation = _income_tax(dataset) + property_slice = np.asarray( + simulation.calculate("taxable_property_income_after_allowances", 2023), + dtype=float, + ) + reformed, _ = _income_tax(dataset, scottish_property_rate_reform(0.01)) + scottish = np.array([region == "SCOTLAND" for region in REGIONS]) + np.testing.assert_allclose( + reformed - baseline, np.where(scottish, 0.01 * property_slice, 0.0), atol=1e-6 + ) + assert (reformed - baseline)[scottish].min() > 0 + + +def test_the_relief_variant_matches_the_engine_at_zero_and_follows_the_rate(): + dataset = _dataset(2023) + baseline, simulation = _income_tax(dataset) + unchanged, _ = _income_tax( + dataset, scottish_property_rate_reform(0.0, relief_follows=True) + ) + np.testing.assert_allclose(unchanged, baseline, atol=1e-6) + + property_slice = np.asarray( + simulation.calculate("taxable_property_income_after_allowances", 2023), + dtype=float, + ) + relieved = np.asarray( + simulation.calculate("property_finance_costs_relieved", 2023), dtype=float + ) + assert relieved[1] > 0 + following, _ = _income_tax( + dataset, scottish_property_rate_reform(0.01, relief_follows=True) + ) + expected = np.array( + [0.01 * property_slice[0], 0.01 * (property_slice[1] - relieved[1]), 0.0] + ) + np.testing.assert_allclose(following - baseline, expected, atol=1e-6) diff --git a/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py b/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py new file mode 100644 index 000000000..edc3c1836 --- /dev/null +++ b/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py @@ -0,0 +1,283 @@ +"""The microcosm#1326 tool's engine-free pieces: bands, published values, the +reweight, the record audit and the disclosure statuses.""" + +from __future__ import annotations + +import json + +import numpy as np +import pandas as pd +import pytest +from scipy import sparse + +from tools.uk_scotland_property_income import measures, rake, registry, release +from tools.uk_scotland_property_income.bands import ( + INTENTIONAL_EDGES, + REQUESTED_EDGES, + SCOTTISH_RATE_STARTS_2023, + assign_bands, + band_labels, + coarse_index_of_fine, +) + +AMOUNTS = ("taxable_property_income",) + + +def test_bands_are_closed_below_and_open_above(): + values = [12_569.99, 12_570, 14_999.99, 15_000, 89_999, 90_000, 2e6] + assert assign_bands(values, REQUESTED_EDGES).tolist() == [-1, 0, 0, 1, 8, 9, 9] + assert band_labels(REQUESTED_EDGES)[0] == "PA–£14,999" + assert band_labels(REQUESTED_EDGES)[-1] == "£90,000+" + assert band_labels(INTENTIONAL_EDGES) == [ + "PA–£24,999", + "£25,000–£39,999", + "£40,000–£89,999", + "£90,000+", + ] + + +def test_release_2_bins_are_unions_of_the_requested_bands(): + assert coarse_index_of_fine(REQUESTED_EDGES, INTENTIONAL_EDGES).tolist() == [ + 0, + 0, + 0, + 1, + 1, + 1, + 2, + 2, + 2, + 3, + ] + with pytest.raises(ValueError, match="not unions"): + coarse_index_of_fine(REQUESTED_EDGES, (12_570, 26_000)) + + +def test_each_bin_edge_is_the_requested_edge_nearest_a_rate_start(): + starts = [SCOTTISH_RATE_STARTS_2023[name] for name in ("intermediate", "higher")] + for edge, start in zip(INTENTIONAL_EDGES[1:3], starts, strict=True): + assert edge == min( + REQUESTED_EDGES, key=lambda candidate: abs(candidate - start) + ) + + +def _spec(name, value, **metadata): + return { + "name": name, + "value": value, + "metadata": {"uprating_from_period": "2023", **metadata}, + } + + +def test_published_value_divides_out_the_uprating_and_the_cross_grain_control(): + plain = _spec( + "hmrc/property_income_income_band_20_000_to_30_000", + 3_831.0, + uprating_factor="1.0025", + ) + assert registry.published_value(plain) == pytest.approx(3_831.0 / 1.0025) + controlled = _spec( + "hmrc.spi_region.taxpayers_by_region_20000_30000@S92000003", + 849_390.0, + uprating_factor="1.00245", + cross_grain_value_before="843061.0", + ) + assert registry.published_value(controlled) == pytest.approx(843_061.0 / 1.00245) + with pytest.raises(ValueError, match="not a 2023-24 fact"): + registry.published_value( + {**plain, "metadata": {"uprating_from_period": "2024"}} + ) + + +def test_spi_facts_read_bands_measures_and_regions_from_the_registry(): + predicate = json.dumps( + {"variable": "region", "operator": "==", "value": "SCOTLAND"} + ) + specs = [ + _spec( + "hmrc/property_income_income_band_1_000_000_to_inf", + 4.0, + uprating_factor="2", + ), + _spec( + "hmrc/property_income_count_income_band_12_570_to_15_000", + 9.0, + uprating_factor="3", + ), + _spec( + "hmrc.spi_region.income_tax_by_region_200000_plus@S92000003", + 10.0, + uprating_factor="2", + geography_predicate=predicate, + ), + _spec("obr.income_tax", 1.0, uprating_factor="1"), + ] + facts = registry.spi_2023_facts({"specs": specs}) + assert [(f.measure, f.lower, f.upper, f.region, f.value) for f in facts] == [ + ("property_amount", 1e6, np.inf, None, 2.0), + ("property_count", 12_570, 15_000, None, 3.0), + ("income_tax", 200_000, np.inf, "SCOTLAND", 5.0), + ] + + +def _persons(): + return pd.DataFrame( + { + "household_row": [0, 0, 1, 2], + "region": ["SCOTLAND", "SCOTLAND", "SCOTLAND", "LONDON"], + "income_tax": [1_000.0, 0.0, 2_000.0, 3_000.0], + "property_income": [5_000.0, 7_000.0, 0.0, 9_000.0], + "total_income_after_loss_relief": [25_000.0, 22_000.0, 28_000.0, 26_000.0], + } + ) + + +def _band_fact(measure, region=None): + return registry.SpiFact("x", measure, 20_000, 30_000, region, 1.0) + + +def _contribution(persons, measure, region=None): + return measures.contribution(persons, _band_fact(measure, region)).tolist() + + +def test_contributions_follow_the_build_bindings(): + persons = _persons() + fact = _band_fact + # Property amounts and counts are taxpayers' only. + assert _contribution(persons, "property_amount") == [5_000, 0, 0, 9_000] + assert _contribution(persons, "property_count") == [1, 0, 0, 1] + # Regional total income has no taxpayer filter; the region predicate applies. + assert _contribution(persons, "total_income", "SCOTLAND") == [ + 25_000, + 22_000, + 28_000, + 0, + ] + assert _contribution(persons, "taxpayers", "SCOTLAND") == [1, 0, 1, 0] + matrix = measures.household_matrix(persons, [fact("property_amount")], 3) + assert matrix.toarray().ravel().tolist() == [5_000, 0, 9_000] + + +def test_total_income_reads_trading_profit_after_losses_brought_forward(): + after = measures.total_income_after_loss_relief( + [30_000, 30_000], [10_000, 10_000], [4_000, 15_000] + ) + assert after.tolist() == [26_000, 20_000] + + +def test_rake_meets_reachable_facts_within_its_bounds(): + rng = np.random.default_rng(1) + matrix = sparse.csr_matrix( + rng.uniform(0, 1, size=(400, 5)) * (rng.uniform(size=(400, 5)) < 0.4) + ) + initial = rng.uniform(50, 150, size=400) + true_ratio = rng.uniform(0.7, 1.4, size=400) + targets = matrix.T @ (initial * true_ratio) + result = rake.rake(matrix, initial, targets, lower=0.5, upper=2.0) + assert np.abs(result.relative_error_after).max() < 1e-5 + assert result.ratio.min() >= 0.5 and result.ratio.max() <= 2.0 + + +def test_rake_stops_at_its_bounds_when_a_fact_is_out_of_reach(): + matrix = sparse.csr_matrix(np.ones((10, 1))) + result = rake.rake(matrix, np.ones(10), [100.0], lower=0.5, upper=2.0) + assert np.allclose(result.ratio, 2.0) + assert result.after[0] == pytest.approx(20.0) + + +def test_source_keys_count_clones_and_copies_once(): + table = pd.DataFrame( + { + "person_support_channel": [b"frs", b"frs", b"spi", b"spi", b"spi"], + "person_source_id": [1001, 1001, 2001, 2001, 2001], + **{name: [0.0] * 5 for name in release.SOURCE_INCOME_COLUMNS}, + } + ) + table.loc[[2, 3], "property_income"] = 12_000.0 + table.loc[4, "property_income"] = 13_000.0 + keys = release.source_keys(table) + assert keys[0] == keys[1] == "frs:1001" + # SPI rows share their FRS host's id, so their incomes tell them apart. + assert keys[2] == keys[3] != keys[4] + + +def _cells(sources_per_band, rows_per_source=3, band_values=(13_000, 16_000, 21_000)): + rows = [] + for band_value, n_sources in zip(band_values, sources_per_band, strict=True): + for source in range(n_sources): + for _clone in range(rows_per_source): + rows.append( + { + "nsnd": band_value, + "source": f"{band_value}:{source}", + "w": 10.0, + "taxable_property_income": 1_000.0 + source, + } + ) + return pd.DataFrame(rows) + + +def test_audit_counts_sources_not_rows(): + edges = (12_570, 15_000, 20_000) + table = release.audit(_cells((4, 12, 20)), edges, AMOUNTS) + assert table["rows"].tolist() == [12, 36, 60] + assert table["sources"].tolist() == [4, 12, 20] + # Equal weight per source: the effective sample size is the source count. + assert table["source_ess"].tolist() == pytest.approx([4, 12, 20]) + + +def test_bootstrap_moves_clones_together(): + cells = _cells((1, 0, 0), rows_per_source=5) + intervals = release.bootstrap( + cells, (12_570, 15_000, 20_000), AMOUNTS, replicates=50 + ) + # A single source is drawn once in every replicate, so its band never varies. + assert intervals.loc[0, "count_low"] == intervals.loc[0, "count_high"] == 50.0 + + +FEW, SECOND, PUBLISHED = release.FEWER_THAN_FLOOR, release.SECONDARY, release.PUBLISHED +FINE_EDGES = (12_570, 15_000, 20_000) + + +def test_a_lone_withheld_band_is_protected_inside_its_bin(): + statuses = release.disclosure_statuses( + _cells((4, 12, 20)), FINE_EDGES, (12_570,), AMOUNTS + ) + assert statuses == ([FEW, SECOND, PUBLISHED], [PUBLISHED]) + + +def test_a_published_total_protects_a_lone_withheld_bin(): + cells = _cells((4, 12, 20)) + # The first bin is the thin band alone. With the total published, it + # equals the total less the other bin, so that bin is withheld too, and + # one of its bands, so it can't be rebuilt from them. + assert release.disclosure_statuses( + cells, FINE_EDGES, (12_570, 15_000), AMOUNTS + ) == ([FEW, SECOND, PUBLISHED], [FEW, SECOND]) + assert release.disclosure_statuses( + cells, FINE_EDGES, (12_570, 15_000), AMOUNTS, total_published=False + ) == ([FEW, PUBLISHED, PUBLISHED], [FEW, PUBLISHED]) + + +def test_one_dominant_record_withholds_its_band(): + cells = _cells((0, 12, 0)) + cells.loc[cells["source"] == "16000:0", "w"] = 1_000.0 + fine_status, _ = release.disclosure_statuses(cells, FINE_EDGES, (12_570,), AMOUNTS) + assert fine_status[1] == release.DOMINATED + + +def test_public_tables_blank_withheld_bands_and_floor_record_counts(): + edges = (12_570, 15_000, 20_000) + cells = _cells((4, 12, 20)) + table = release.audit(cells, edges, AMOUNTS) + intervals = release.bootstrap(cells, edges, AMOUNTS, replicates=20) + statuses = [release.FEWER_THAN_FLOOR, release.SECONDARY, release.PUBLISHED] + public = release.render(table, intervals, statuses, AMOUNTS, mode="public") + assert public["distinct_records"].tolist() == ["<10", "12", "20"] + assert pd.isna(public.loc[0, "count"]) and pd.isna( + public.loc[1, "taxable_property_income"] + ) + assert public.loc[2, "count"] == 600.0 + restricted = release.render(table, intervals, statuses, AMOUNTS, mode="restricted") + assert restricted["distinct_records"].tolist() == ["4", "12", "20"] + assert restricted.loc[0, "count"] == 100.0 diff --git a/tools/uk_scotland_property_income/__init__.py b/tools/uk_scotland_property_income/__init__.py new file mode 100644 index 000000000..d086a4947 --- /dev/null +++ b/tools/uk_scotland_property_income/__init__.py @@ -0,0 +1,17 @@ +"""2023-24 property income of Scottish taxpayers by NSND band (microcosm#1326). + +A project tool for one external request, kept off main. It reads a released UK +national file and the national target registry of the build that produced it, +and writes private outputs outside the repository: + +1. a one-year back-cast of the 2024-25 file with policyengine-uk's own growth + indices (``backcast``); +2. a bounded reweight to the SPI 2023-24 facts the build already binds + (``registry``, ``measures``, ``rake``); +3. the two release tables with the records behind every band, and their + disclosure statuses (``bands``, ``release``); +4. an illustrative 1p on the Scottish rates charged on property income + (``reform``). + +Run it with ``python -m tools.uk_scotland_property_income --help``. +""" diff --git a/tools/uk_scotland_property_income/__main__.py b/tools/uk_scotland_property_income/__main__.py new file mode 100644 index 000000000..2b11a622b --- /dev/null +++ b/tools/uk_scotland_property_income/__main__.py @@ -0,0 +1,424 @@ +"""Run the microcosm#1326 analysis on one released UK national file. + + python -m tools.uk_scotland_property_income \\ + --h5 /microcosm_uk_2024_25.h5 \\ + --registry /national_target_registry.json \\ + --out + +The output directory holds three kinds of file: + +- ``private/``: the per-band audit (rows, exact sources, shares), the reweight + report and the drift on every SPI 2023-24 fact. Licensed working files; + never shared. +- ``restricted/``: both releases in full, with exact record counts. Only for + recipients whose licence covers the FRS and the SPI tape. +- ``public/``: both releases with bands under the floor withheld and flagged, + and secondary suppression across the two. + +Everything is derived from UK Data Service licensed microdata, so the tool +refuses an output directory inside a git working tree. +""" + +from __future__ import annotations + +import argparse +import hashlib +import importlib.metadata +import json +import subprocess +from pathlib import Path +from typing import Any + +import numpy as np +import pandas as pd +from scipy import sparse + +from tools.uk_scotland_property_income import ( + backcast, + measures, + rake, + registry, + release, +) +from tools.uk_scotland_property_income.bands import ( + INTENTIONAL_EDGES, + PERSONAL_ALLOWANCE_2023, + REQUESTED_EDGES, + SCOTTISH_RATE_STARTS_2023, +) + +REPO_ROOT = Path(__file__).resolve().parents[2] +AMOUNTS = ("taxable_property_income", "property_income") +#: Pounds of float32 rounding allowed in a non-Scottish taxpayer's tax. +OUTSIDE_SCOTLAND_TOLERANCE = 0.05 +ENGINE_VARIABLES = ( + "person_id", + "pays_scottish_income_tax", + "income_tax", + "property_income", + "taxable_property_income", + "taxable_property_income_after_allowances", + "net_income", + "taxable_savings_interest_income", + "taxable_dividend_income", + "total_income", + "self_employment_income", + "loss_relief", +) + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for block in iter(lambda: handle.read(1 << 20), b""): + digest.update(block) + return digest.hexdigest() + + +def _refuse_working_tree(out: Path) -> None: + probe = out.resolve() + while not probe.exists(): + probe = probe.parent + inside = subprocess.run( + ["git", "-C", str(probe), "rev-parse", "--is-inside-work-tree"], + capture_output=True, + text=True, + ) + if inside.returncode == 0 and inside.stdout.strip() == "true": + raise SystemExit( + f"{out} is inside a git working tree; write licensed outputs elsewhere." + ) + + +def _git_state() -> dict[str, Any]: + def git(*args: str) -> str: + return subprocess.run( + ["git", "-C", str(REPO_ROOT), *args], capture_output=True, text=True + ).stdout.strip() + + return { + "commit": git("rev-parse", "HEAD"), + "dirty": bool( + git("status", "--porcelain", "--", "tools/uk_scotland_property_income") + ), + } + + +def _text(values: pd.Series) -> np.ndarray: + return values.map( + lambda value: value.decode() if isinstance(value, bytes) else str(value) + ).to_numpy() + + +def person_frame(simulation: Any, dataset: Any, year: int) -> pd.DataFrame: + """The person columns every step reads, at ``year``.""" + + backcast.refuse_earlier_period(dataset, year) + values = { + name: np.asarray(simulation.calculate(name, year), dtype=float) + for name in ENGINE_VARIABLES + } + person = dataset.person + if not np.array_equal( + values["person_id"], person["person_id"].to_numpy(dtype=float) + ): + raise ValueError("The engine's person order differs from the dataset's.") + household = dataset.household + position = pd.Series( + np.arange(len(household)), index=household["household_id"].to_numpy() + ) + household_row = position.loc[person["person_household_id"].to_numpy()].to_numpy() + frame = pd.DataFrame( + {name: values[name] for name in ENGINE_VARIABLES if name != "person_id"} + ) + frame.insert(0, "household_row", household_row) + frame.insert(1, "region", _text(household["region"])[household_row]) + frame["scottish"] = frame.pop("pays_scottish_income_tax") > 0 + frame["nsnd"] = ( + frame["net_income"] + - frame["taxable_savings_interest_income"] + - frame["taxable_dividend_income"] + ) + frame["total_income_after_loss_relief"] = measures.total_income_after_loss_relief( + frame["total_income"], frame["self_employment_income"], frame["loss_relief"] + ) + return frame + + +def _write(frame: pd.DataFrame, path: Path) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + frame.to_csv(path, index=False) + + +def run(args: argparse.Namespace) -> dict[str, Any]: + from policyengine_uk import CountryTaxBenefitSystem, Microsimulation + from policyengine_uk.data import UKSingleYearDataset + + from tools.uk_scotland_property_income.reform import scottish_property_rate_reform + + out = Path(args.out) + _refuse_working_tree(out) + h5, registry_path = Path(args.h5), Path(args.registry) + h5_sha = _sha256(h5) + if args.expect_sha256 and not h5_sha.startswith(args.expect_sha256): + raise SystemExit(f"{h5} has sha256 {h5_sha}, not {args.expect_sha256}.") + + released = UKSingleYearDataset(file_path=str(h5)) + if int(released.time_period) != args.year + 1: + raise SystemExit( + f"{h5} is for {released.time_period}; this back-casts one year to {args.year}." + ) + sources = release.source_keys(released.person) + parameters = CountryTaxBenefitSystem().parameters + earlier, divided = backcast.backcast_one_year(released, parameters) + w0 = earlier.household["household_weight"].to_numpy(dtype=float) + persons = person_frame( + Microsimulation(dataset=backcast.with_household_weights(earlier, w0)), + earlier, + args.year, + ) + persons["source"] = sources + + # A4: reweight to the SPI 2023-24 facts, holding each region's households. + facts = registry.spi_2023_facts(registry.load_registry(registry_path)) + raked_facts = [fact for fact in facts if fact.region is None] + [ + fact + for fact in facts + if fact.region == "SCOTLAND" and fact.measure in {"taxpayers", "total_income"} + ] + household_regions = _text(earlier.household["region"]) + regions = sorted(set(household_regions)) + n_households = len(w0) + matrix = sparse.hstack( + [ + measures.household_matrix(persons, raked_facts, n_households), + measures.region_matrix(regions, household_regions), + ], + format="csr", + ) + targets = np.concatenate( + [ + [fact.value for fact in raked_facts], + [w0[household_regions == region].sum() for region in regions], + ] + ) + raked = rake.rake( + matrix, w0, targets, lower=args.ratio_bounds[0], upper=args.ratio_bounds[1] + ) + w1 = raked.weights + all_matrix = measures.household_matrix(persons, facts, n_households) + raked_names = {fact.name for fact in raked_facts} + drift = pd.DataFrame( + { + "name": [fact.name for fact in facts], + "measure": [fact.measure for fact in facts], + "region": [fact.region or "UK" for fact in facts], + "lower": [fact.lower for fact in facts], + "upper": [fact.upper for fact in facts], + "published_2023_24": [fact.value for fact in facts], + "backcast": all_matrix.T @ w0, + "reweighted": all_matrix.T @ w1, + "reweighted_to": [fact.name in raked_names for fact in facts], + } + ) + drift["backcast_error"] = drift["backcast"] / drift["published_2023_24"] - 1 + drift["reweighted_error"] = drift["reweighted"] / drift["published_2023_24"] - 1 + + # A6: the two releases. + persons["w"] = w1[persons["household_row"].to_numpy()] + cells = persons[ + persons["scottish"] + & (persons["income_tax"] > 0) + & (persons["property_income"] > 0) + & (persons["nsnd"] >= PERSONAL_ALLOWANCE_2023) + ] + tables = {} + for key, edges in ( + ("release_1", REQUESTED_EDGES), + ("release_2", INTENTIONAL_EDGES), + ): + tables[key] = ( + release.audit(cells, edges, AMOUNTS), + release.bootstrap( + cells, edges, AMOUNTS, replicates=args.replicates, seed=args.seed + ), + ) + statuses = dict( + zip( + ("release_1", "release_2"), + release.disclosure_statuses( + cells, REQUESTED_EDGES, INTENTIONAL_EDGES, AMOUNTS + ), + strict=True, + ) + ) + for key, (table, intervals) in tables.items(): + audit = table.merge(intervals, on="band").assign(status=statuses[key]) + _write(audit, out / "private" / f"audit_{key}.csv") + for mode in ("restricted", "public"): + _write( + release.render(table, intervals, statuses[key], AMOUNTS, mode=mode), + out / mode / f"{key}.csv", + ) + + # B: the illustrative 1p, with and without the finance-cost reduction following. + scottish = persons["scottish"].to_numpy() + costing = {} + for variant, relief_follows in (("rates_only", False), ("relief_follows", True)): + reformed = Microsimulation( + dataset=backcast.with_household_weights(earlier, w0), + reform=scottish_property_rate_reform( + args.delta, relief_follows=relief_follows + ), + ) + change = ( + np.asarray(reformed.calculate("income_tax", args.year), dtype=float) + - persons["income_tax"].to_numpy() + ) + # The engine stores float32, and the relief variant's rate is a float64 + # array, so tax outside Scotland can move by a rounding penny. + if np.abs(change[~scottish]).max(initial=0.0) > OUTSIDE_SCOTLAND_TOLERANCE: + raise ValueError("The Scottish reform changed tax outside Scotland.") + change = np.where(scottish, change, 0.0) + weighted = persons["w"].to_numpy() * change + summary = release.bootstrap_total( + weighted[scottish], + persons.loc[scottish, "source"], + replicates=args.replicates, + seed=args.seed, + ) + affected = scottish & (change > 0) + summary["taxpayers_affected"] = float(persons.loc[affected, "w"].sum()) + summary["sources_affected"] = int(persons.loc[affected, "source"].nunique()) + costing[variant] = summary + _write_costing(costing, out) + + manifest = { + "issue": "PolicyEngine/microcosm#1326", + "year": args.year, + "inputs": { + "h5": {"path": str(h5.resolve()), "sha256": h5_sha}, + "registry": { + "path": str(registry_path.resolve()), + "sha256": _sha256(registry_path), + }, + }, + "engine": { + name: importlib.metadata.version(name) + for name in ("policyengine-uk", "policyengine-core") + }, + "tool": _git_state(), + "backcast": {"divided": divided, "not_inverted": list(backcast.NOT_INVERTED)}, + "reweight": { + "facts": [fact.name for fact in raked_facts], + "region_controls": regions, + "ratio_bounds": list(args.ratio_bounds), + "converged": raked.converged, + "message": raked.message, + "max_abs_error_before": float(np.abs(raked.relative_error_before).max()), + "max_abs_error_after": float(np.abs(raked.relative_error_after).max()), + "ratio_quantiles": dict( + zip( + ("min", "p01", "median", "p99", "max"), + np.quantile(raked.ratio, [0, 0.01, 0.5, 0.99, 1]).tolist(), + strict=True, + ) + ), + "share_at_bound": float( + np.mean( + (raked.ratio <= args.ratio_bounds[0] * 1.0001) + | (raked.ratio >= args.ratio_bounds[1] * 0.9999) + ) + ), + "effective_sample_size_ratio": raked.effective_sample_size_ratio, + }, + "bands": { + "release_1": list(REQUESTED_EDGES), + "release_2": list(INTENTIONAL_EDGES), + "scottish_rate_starts_2023": SCOTTISH_RATE_STARTS_2023, + }, + "bootstrap": { + "replicates": args.replicates, + "seed": args.seed, + "weights": "fixed", + }, + "delta": args.delta, + } + (out / "private").mkdir(parents=True, exist_ok=True) + drift.to_csv(out / "private" / "drift.csv", index=False) + (out / "manifest.json").write_text(json.dumps(manifest, indent=2) + "\n") + return manifest + + +def _write_costing(costing: dict[str, dict[str, float]], out: Path) -> None: + (out / "private").mkdir(parents=True, exist_ok=True) + (out / "private" / "costing.json").write_text(json.dumps(costing, indent=2) + "\n") + rows = [] + for variant, summary in costing.items(): + published = ( + summary["sources"] >= release.MIN_SOURCES + and summary["top_source_share"] <= release.DOMINANCE_SHARE + ) + rows.append( + { + "variant": variant, + "change_in_scottish_income_tax": round(summary["estimate"] / 1e5) * 1e5, + "low": round(summary["low"] / 1e5) * 1e5, + "high": round(summary["high"] / 1e5) * 1e5, + "cv": round(summary["cv"], 2), + "taxpayers_affected": round(summary["taxpayers_affected"] / 100) * 100, + "status": release.PUBLISHED if published else release.DOMINATED, + } + ) + table = pd.DataFrame(rows) + _write(table, out / "restricted" / "costing.csv") + public = table.copy() + withheld = public["status"] != release.PUBLISHED + public.loc[ + withheld, + ["change_in_scottish_income_tax", "low", "high", "cv", "taxpayers_affected"], + ] = None + _write(public, out / "public" / "costing.csv") + + +def parse_args(argv: list[str] | None = None) -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + parser.add_argument( + "--h5", required=True, help="The released 2024-25 national file." + ) + parser.add_argument( + "--registry", required=True, help="The build's national_target_registry.json." + ) + parser.add_argument( + "--out", required=True, help="A private directory outside any git working tree." + ) + parser.add_argument( + "--year", + type=int, + default=2023, + help="The year to back-cast to (default 2023).", + ) + parser.add_argument( + "--expect-sha256", help="Refuse an H5 whose sha256 doesn't start with this." + ) + parser.add_argument("--ratio-bounds", type=float, nargs=2, default=(0.25, 4.0)) + parser.add_argument("--replicates", type=int, default=1_000) + parser.add_argument("--seed", type=int, default=0) + parser.add_argument( + "--delta", + type=float, + default=0.01, + help="The change in each Scottish rate (default 1p).", + ) + return parser.parse_args(argv) + + +if __name__ == "__main__": + manifest = run(parse_args()) + reweight = manifest["reweight"] + print( + f"done: reweight converged={reweight['converged']}, " + f"max |error| {reweight['max_abs_error_before']:.3f} -> {reweight['max_abs_error_after']:.4f}" + ) diff --git a/tools/uk_scotland_property_income/backcast.py b/tools/uk_scotland_property_income/backcast.py new file mode 100644 index 000000000..ef0491195 --- /dev/null +++ b/tools/uk_scotland_property_income/backcast.py @@ -0,0 +1,86 @@ +"""One-year back-cast of a UK dataset with policyengine-uk's own indices (microcosm#1326 A2, A3). + +policyengine-uk projects a dataset forward only. Asked for a year before the +dataset's own, it returns zero incomes, zero tax and almost no weight rather +than refusing. So the only way to 2023-24 from a 2024-25 file is to invert +the engine's own forward step: divide each input by the growth +``apply_single_year_uprating`` would multiply it by. Household weights move +with ONS population growth, as they do going forward. + +The forward step's custom pieces are not inverted. Council tax and rent move +by region, and student-loan plans are re-labelled; none of them enters income +tax, and they stay at the file's own year. +""" + +from __future__ import annotations + +from typing import Any + +from policyengine_uk.data import UKSingleYearDataset +from policyengine_uk.data.economic_assumptions import ( + LEGACY_PROJECTION_END_YEAR, + _uprating_index_map, +) + +#: Inputs the forward step moves by its own rules, which this back-cast leaves +#: at the file's year. +NOT_INVERTED = ("council_tax", "rent", "student_loan_plan") + + +def refuse_earlier_period(dataset: UKSingleYearDataset, year: int) -> None: + """Stop before asking the engine for a year it would silently zero.""" + + if int(year) < int(dataset.time_period): + raise ValueError( + f"The dataset is for {dataset.time_period}; policyengine-uk returns zero " + f"incomes and weights for {year} instead of refusing. Back-cast it first." + ) + + +def backcast_one_year( + dataset: UKSingleYearDataset, parameters: Any +) -> tuple[UKSingleYearDataset, dict[str, dict[str, Any]]]: + """The dataset one year earlier, and the growth divided out of each input.""" + + year = int(dataset.time_period) + if year > LEGACY_PROJECTION_END_YEAR: + raise ValueError( + f"Inputs outside uprating_indices.yaml also move after {LEGACY_PROJECTION_END_YEAR}; " + "this back-cast inverts only the mapped ones." + ) + tables = { + "person": dataset.person.copy(), + "benunit": dataset.benunit.copy(), + "household": dataset.household.copy(), + } + divided: dict[str, dict[str, Any]] = {} + for index_name, variables in _uprating_index_map().items(): + # The forward step reads the growth at the later year's start. + growth = float(parameters.get_child(index_name)(str(year))) + for variable in variables: + for table in tables.values(): + if variable in table.columns: + table[variable] = table[variable] / (1.0 + growth) + divided[variable] = {"index": index_name, "growth": growth} + earlier = UKSingleYearDataset( + person=tables["person"], + benunit=tables["benunit"], + household=tables["household"], + fiscal_year=year - 1, + ) + return earlier, divided + + +def with_household_weights( + dataset: UKSingleYearDataset, weights: Any +) -> UKSingleYearDataset: + """A copy of the dataset with new household weights.""" + + household = dataset.household.copy() + household["household_weight"] = weights + return UKSingleYearDataset( + person=dataset.person.copy(), + benunit=dataset.benunit.copy(), + household=household, + fiscal_year=int(dataset.time_period), + ) diff --git a/tools/uk_scotland_property_income/bands.py b/tools/uk_scotland_property_income/bands.py new file mode 100644 index 000000000..434f7d4ca --- /dev/null +++ b/tools/uk_scotland_property_income/bands.py @@ -0,0 +1,79 @@ +"""The NSND bands of the two releases (microcosm#1326 A6). + +Bands are on individual non-savings, non-dividend income before the personal +allowance, closed below and open above. Release 1 keeps the ten bands the +request names. Release 2's bins are unions of them, each edge the requested +edge nearest a 2023-24 Scottish rate threshold, so the outturn can be summed to +match and each bin mostly faces one marginal rate. +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np + +PERSONAL_ALLOWANCE_2023 = 12_570 + +#: Lower edges of the ten requested bands; the last is open. +REQUESTED_EDGES: tuple[int, ...] = ( + 12_570, + 15_000, + 20_000, + 25_000, + 30_000, + 35_000, + 40_000, + 50_000, + 65_000, + 90_000, +) + +#: Lower edges of the Release 2 bins: starter and basic, intermediate, higher, +#: and higher with top. +INTENTIONAL_EDGES: tuple[int, ...] = (12_570, 25_000, 40_000, 90_000) + +#: Where each 2023-24 Scottish rate starts, in NSND income before the +#: allowance: the allowance plus the band's start in taxable income. The top +#: rate starts at 125,140, where the allowance has tapered away. +SCOTTISH_RATE_STARTS_2023: dict[str, int] = { + "starter": 12_571, + "basic": 14_733, + "intermediate": 25_689, + "higher": 43_663, + "top": 125_140, +} + + +def band_labels(edges: Sequence[int]) -> list[str]: + """``PA–£14,999`` style labels, the last one open.""" + + labels = [] + uppers = (*edges[1:], None) + for lower, upper in zip(edges, uppers, strict=True): + start = "PA" if lower == PERSONAL_ALLOWANCE_2023 else f"£{lower:,}" + labels.append(f"{start}+" if upper is None else f"{start}–£{upper - 1:,}") + return labels + + +def assign_bands(values: object, edges: Sequence[int]) -> np.ndarray: + """Each value's band index, or -1 below the first edge.""" + + boundaries = np.asarray(edges, dtype=float) + return ( + np.searchsorted(boundaries, np.asarray(values, dtype=float), side="right") - 1 + ) + + +def coarse_index_of_fine(fine: Sequence[int], coarse: Sequence[int]) -> np.ndarray: + """The coarse band holding each fine band. + + Refuses unless every coarse band is a union of fine bands, which is what + lets the outturn be summed to the coarse bins. + """ + + if coarse[0] != fine[0] or not set(coarse) <= set(fine): + raise ValueError( + f"Bins {list(coarse)} are not unions of the bands {list(fine)}." + ) + return assign_bands(fine, coarse) diff --git a/tools/uk_scotland_property_income/measures.py b/tools/uk_scotland_property_income/measures.py new file mode 100644 index 000000000..a6d6ed6c1 --- /dev/null +++ b/tools/uk_scotland_property_income/measures.py @@ -0,0 +1,106 @@ +"""Person contributions to the SPI facts, on the build's own bindings (microcosm#1326 A4). + +``uk_population_targets.json`` binds each fact as a weighted sum over persons: + +- property amount: ``property_income`` of taxpayers (``income_tax > 0``); +- property count: taxpayers with ``property_income > 0``; +- regional taxpayers: persons with ``income_tax > 0``; +- regional total income: ``total_income_after_loss_relief`` of everyone in the + band (the binding has no taxpayer filter); +- regional income tax: ``income_tax`` of everyone in the band. + +Every band reads ``total_income_after_loss_relief``: the engine's +``total_income`` with self-employment profit after losses brought forward +(``uk_runtime/measure_simulation.py``, microcosm#1095). +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np +import pandas as pd +from scipy import sparse + +from tools.uk_scotland_property_income.registry import SpiFact + +#: Person columns the contributions read. ``household_row`` is the position of +#: the person's household in the household table. +PERSON_COLUMNS = ( + "household_row", + "region", + "income_tax", + "property_income", + "total_income_after_loss_relief", +) + + +def total_income_after_loss_relief( + total_income: object, self_employment_income: object, loss_relief: object +) -> np.ndarray: + """Total income with the trading profit after losses brought forward.""" + + profit = np.asarray(self_employment_income, dtype=float) + after = np.maximum(profit - np.asarray(loss_relief, dtype=float), 0.0) + return np.asarray(total_income, dtype=float) - profit + after + + +def contribution(persons: pd.DataFrame, fact: SpiFact) -> np.ndarray: + """Each person's unweighted contribution to the fact.""" + + income = persons["total_income_after_loss_relief"].to_numpy(dtype=float) + selected = (income >= fact.lower) & (income < fact.upper) + if fact.region is not None: + selected &= persons["region"].to_numpy() == fact.region + tax = persons["income_tax"].to_numpy(dtype=float) + profit = persons["property_income"].to_numpy(dtype=float) + if fact.measure == "property_amount": + value = np.where(tax > 0, profit, 0.0) + elif fact.measure == "property_count": + value = ((tax > 0) & (profit > 0)).astype(float) + elif fact.measure == "taxpayers": + value = (tax > 0).astype(float) + elif fact.measure == "total_income": + value = income + elif fact.measure == "income_tax": + value = tax + else: + raise ValueError(f"Unknown measure {fact.measure!r}.") + return np.where(selected, value, 0.0) + + +def household_matrix( + persons: pd.DataFrame, facts: Sequence[SpiFact], n_households: int +) -> sparse.csr_matrix: + """Households by facts: each household's summed contribution. + + A household's weight multiplies every member, so the fact's estimate is + ``matrix.T @ household_weights``. + """ + + households = persons["household_row"].to_numpy(dtype=np.int64) + rows, cols, data = [], [], [] + for column, fact in enumerate(facts): + values = contribution(persons, fact) + nonzero = values != 0 + rows.append(households[nonzero]) + cols.append(np.full(int(nonzero.sum()), column, dtype=np.int64)) + data.append(values[nonzero]) + matrix = sparse.coo_matrix( + (np.concatenate(data), (np.concatenate(rows), np.concatenate(cols))), + shape=(n_households, len(facts)), + ) + return matrix.tocsr() + + +def region_matrix( + regions: Sequence[str], household_regions: object +) -> sparse.csr_matrix: + """Households by regions: one where the household sits in the region.""" + + household_regions = np.asarray(household_regions) + columns = [ + sparse.csr_matrix((household_regions == region).astype(float)[:, None]) + for region in regions + ] + return sparse.hstack(columns, format="csr") diff --git a/tools/uk_scotland_property_income/rake.py b/tools/uk_scotland_property_income/rake.py new file mode 100644 index 000000000..2ce904aff --- /dev/null +++ b/tools/uk_scotland_property_income/rake.py @@ -0,0 +1,110 @@ +"""Bounded reweighting of household weights to a handful of facts (microcosm#1326 A4). + +Each household's weight is scaled by a ratio ``g`` in ``[lower, upper]``. The +ratios minimise the squared relative misses on the facts plus a small +divergence from the starting weights, +``sum(w0 * (g log g - g + 1)) / sum(w0)``. When the facts can be met, the +divergence picks the solution closest to the starting weights. When they +can't be met within the bounds, it gives the least distorting compromise +instead of failing. L-BFGS-B searches over ``log g``. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from scipy import sparse +from scipy.optimize import Bounds, minimize + + +@dataclass(frozen=True) +class RakeResult: + weights: np.ndarray + ratio: np.ndarray + targets: np.ndarray + before: np.ndarray + after: np.ndarray + converged: bool + message: str + + @property + def relative_error_before(self) -> np.ndarray: + return self.before / self.targets - 1.0 + + @property + def relative_error_after(self) -> np.ndarray: + return self.after / self.targets - 1.0 + + @property + def effective_sample_size_ratio(self) -> float: + """Kish effective sample size after, over before.""" + + def kish(weights: np.ndarray) -> float: + return float(weights.sum() ** 2 / np.square(weights).sum()) + + return kish(self.weights) / kish(self.weights / self.ratio) + + +def rake( + matrix: sparse.spmatrix, + initial_weights: object, + targets: object, + *, + lower: float = 0.25, + upper: float = 4.0, + divergence_weight: float = 1e-3, + max_iterations: int = 5_000, +) -> RakeResult: + """Ratios on ``initial_weights`` so that ``matrix.T @ weights`` meets ``targets``. + + ``matrix`` is households by facts. + """ + + a = sparse.csr_matrix(matrix, dtype=float) + w0 = np.asarray(initial_weights, dtype=float) + t = np.asarray(targets, dtype=float) + if a.shape != (w0.size, t.size): + raise ValueError( + f"matrix {a.shape} does not match {w0.size} weights and {t.size} facts." + ) + if not 0 < lower <= 1 <= upper: + raise ValueError(f"Ratio bounds must contain 1: [{lower}, {upper}].") + if np.any(w0 <= 0): + raise ValueError("Starting weights must be positive.") + if np.any(t <= 0): + raise ValueError("Every fact must be positive.") + + scale = 1.0 / t + total = w0.sum() + a_t = a.T.tocsr() + + def objective(log_ratio: np.ndarray) -> tuple[float, np.ndarray]: + ratio = np.exp(log_ratio) + weights = w0 * ratio + miss = (a_t @ weights - t) * scale + divergence = np.sum(w0 * (ratio * log_ratio - ratio + 1.0)) / total + value = float(miss @ miss + divergence_weight * divergence) + gradient = weights * (a @ (2.0 * miss * scale)) + gradient += divergence_weight * w0 * ratio * log_ratio / total + return value, gradient + + solution = minimize( + objective, + np.zeros(w0.size), + jac=True, + method="L-BFGS-B", + bounds=Bounds(np.full(w0.size, np.log(lower)), np.full(w0.size, np.log(upper))), + options={"maxiter": max_iterations, "maxcor": 30, "ftol": 1e-15, "gtol": 1e-12}, + ) + ratio = np.clip(np.exp(solution.x), lower, upper) + weights = w0 * ratio + return RakeResult( + weights=weights, + ratio=ratio, + targets=t, + before=a_t @ w0, + after=a_t @ weights, + converged=bool(solution.success), + message=str(solution.message), + ) diff --git a/tools/uk_scotland_property_income/reform.py b/tools/uk_scotland_property_income/reform.py new file mode 100644 index 000000000..281a22963 --- /dev/null +++ b/tools/uk_scotland_property_income/reform.py @@ -0,0 +1,82 @@ +"""An illustrative change to the Scottish rates on property income (microcosm#1326 B). + +policyengine-uk taxes a Scottish taxpayer's property income at the Scottish +rates, stacked on top of their other non-savings income +(``property_income_tax``). It has no separate Scottish property rate. Adding +``delta`` to every Scottish rate on that slice adds exactly ``delta`` times the +slice, ``taxable_property_income_after_allowances``, whichever bands it +crosses. That is what this reform charges. + +The finance-cost reduction is given at the UK property basic rate +(``property_finance_cost_relief``). By default it stays there. With +``relief_follows``, Scottish taxpayers' reduction rises by ``delta`` too, +within the same Step 6 limit the engine applies. +""" + +from __future__ import annotations + +from policyengine_core.reforms import Reform +from policyengine_uk.model_api import ( + GBP, + YEAR, + Person, + Variable, + add, + max_, + min_, + where, +) +from policyengine_uk.variables.gov.hmrc.income_tax.liability.property_income_tax import ( + property_income_tax as _engine_property_income_tax, +) + + +def scottish_property_rate_reform( + delta: float = 0.01, *, relief_follows: bool = False +) -> type[Reform]: + """A reform adding ``delta`` to the Scottish rates charged on property income.""" + + class property_income_tax(Variable): # noqa: N801 (engine variable names) + value_type = float + entity = Person + label = "Income tax on property income" + definition_period = YEAR + unit = GBP + + def formula(person, period, parameters): # noqa: N805 + charged = _engine_property_income_tax.formula(person, period, parameters) + scottish = person("pays_scottish_income_tax", period) + property_slice = person("taxable_property_income_after_allowances", period) + return charged + where(scottish, delta * property_slice, 0) + + class property_finance_cost_relief(Variable): # noqa: N801 + value_type = float + entity = Person + label = "residential property finance costs tax reduction" + definition_period = YEAR + unit = GBP + + def formula(person, period, parameters): # noqa: N805 + income_tax = parameters(period).gov.hmrc.income_tax + scottish = person("pays_scottish_income_tax", period) + rate = income_tax.rates.property.basic + where(scottish, delta, 0) + reduction = rate * person("property_finance_costs_relieved", period) + other_reductions = [ + variable + for variable in income_tax.income_tax_subtractions + if variable != "property_finance_cost_relief" + ] + tax_left = max_( + 0, + person("income_tax_pre_charges", period) + - add(person, period, other_reductions), + ) + return min_(reduction, tax_left) + + class ScottishPropertyRateReform(Reform): + def apply(self): + self.update_variable(property_income_tax) + if relief_follows: + self.update_variable(property_finance_cost_relief) + + return ScottishPropertyRateReform diff --git a/tools/uk_scotland_property_income/registry.py b/tools/uk_scotland_property_income/registry.py new file mode 100644 index 000000000..e6c5a2560 --- /dev/null +++ b/tools/uk_scotland_property_income/registry.py @@ -0,0 +1,101 @@ +"""The SPI 2023-24 facts the national build binds, read back to 2023-24 (microcosm#1326 A4). + +The national target registry holds each SPI fact uprated to the calibration +period. Its metadata records the factor (``uprating_factor``) and, for the +regional rows, the value before the cross-grain control +(``cross_grain_value_before``), so the published 2023-24 value is the value +before the control divided by the factor. +""" + +from __future__ import annotations + +import json +import re +from collections.abc import Mapping +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +#: The measures this tool evaluates, as ``measures.contribution`` defines them. +MEASURES = ( + "property_amount", + "property_count", + "taxpayers", + "total_income", + "income_tax", +) + +_PROPERTY = re.compile( + r"^hmrc/property_income_(?Pcount_)?income_band_" + r"(?P[0-9_]+)_to_(?P[0-9_]+|inf)$" +) +_REGION = re.compile( + r"^hmrc\.spi_region\.(?Ptaxpayers|total_income|income_tax)_by_region_" + r"(?P\d+)_(?P\d+|plus)@" +) + + +@dataclass(frozen=True) +class SpiFact: + """One SPI 2023-24 fact at its published value.""" + + name: str + measure: str + lower: float + upper: float + region: str | None + value: float + + +def load_registry(path: str | Path) -> dict[str, Any]: + return json.loads(Path(path).read_text()) + + +def published_value(spec: Mapping[str, Any]) -> float: + """The fact's 2023-24 value, before uprating and any cross-grain control.""" + + metadata = spec.get("metadata") or {} + if str(metadata.get("uprating_from_period")) != "2023": + raise ValueError(f"{spec['name']} is not a 2023-24 fact.") + if "uprating_factor" not in metadata: + raise ValueError(f"{spec['name']} records no uprating factor.") + base = float(metadata.get("cross_grain_value_before", spec["value"])) + return base / float(metadata["uprating_factor"]) + + +def _amount(text: str) -> float: + return float("inf") if text in {"inf", "plus"} else float(text.replace("_", "")) + + +def _region(spec: Mapping[str, Any]) -> str: + predicate = json.loads((spec.get("metadata") or {})["geography_predicate"]) + if predicate.get("variable") != "region" or predicate.get("operator") != "==": + raise ValueError(f"{spec['name']} has an unexpected geography predicate.") + return str(predicate["value"]) + + +def spi_2023_facts(registry: Mapping[str, Any]) -> list[SpiFact]: + """Every SPI property and regional fact in the registry, at 2023-24 values.""" + + facts = [] + for spec in registry["specs"]: + name = str(spec["name"]) + if match := _PROPERTY.match(name): + measure = "property_count" if match["count"] else "property_amount" + region = None + elif match := _REGION.match(name): + measure = match["measure"] + region = _region(spec) + else: + continue + facts.append( + SpiFact( + name=name, + measure=measure, + lower=_amount(match["lower"]), + upper=_amount(match["upper"]), + region=region, + value=published_value(spec), + ) + ) + return facts diff --git a/tools/uk_scotland_property_income/release.py b/tools/uk_scotland_property_income/release.py new file mode 100644 index 000000000..efef22461 --- /dev/null +++ b/tools/uk_scotland_property_income/release.py @@ -0,0 +1,363 @@ +"""The release tables, the records behind them and their disclosure statuses (microcosm#1326 A6). + +Rows repeat records: CGT clones, support copies and, in a dense file, area +clones all carry one source's values. So every count here is of distinct +sources, every interval resamples sources with their clones together, and the +disclosure rules read sources, never rows. +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np +import pandas as pd + +from tools.uk_scotland_property_income.bands import ( + assign_bands, + band_labels, + coarse_index_of_fine, +) + +#: The stored incomes that identify an SPI-channel record. Its tape record isn't +#: carried (microcosm#1287), so records with the same income vector count once. +SOURCE_INCOME_COLUMNS = ( + "property_income", + "employment_income", + "self_employment_income", + "private_pension_income", + "savings_interest_income", + "dividend_income", +) + +MIN_SOURCES = 10 +DOMINANCE_SHARE = 0.5 + +PUBLISHED = "published" +FEWER_THAN_FLOOR = "withheld: fewer than 10 records" +DOMINATED = "withheld: one record dominates" +SECONDARY = "withheld: protects another band" + + +def _text(values: pd.Series) -> pd.Series: + return values.map( + lambda value: value.decode() if isinstance(value, bytes) else str(value) + ) + + +def source_keys(person_table: pd.DataFrame) -> np.ndarray: + """One key per independent record. + + FRS persons are keyed by ``person_source_id``. SPI-channel persons carry + their FRS host's id, not their tape record's, so they are keyed by their + stored income vector. Read it from the released file, before any + back-cast. + """ + + channel = _text(person_table["person_support_channel"]).to_numpy() + frs = "frs:" + person_table["person_source_id"].astype(np.int64).astype(str) + incomes = ( + person_table[list(SOURCE_INCOME_COLUMNS)].astype(float).round(2).astype(str) + ) + spi = "spi:" + incomes.agg("|".join, axis=1) + return np.where(channel == "frs", frs.to_numpy(), spi.to_numpy()) + + +def audit( + cells: pd.DataFrame, edges: Sequence[int], amount_columns: Sequence[str] +) -> pd.DataFrame: + """Per band: rows, sources, effective sample size by source, weighted totals + and the largest source's share of the weight and of each amount. + + ``cells`` has one row per selected person, with ``w``, ``nsnd``, ``source`` + and the amount columns. + """ + + band = assign_bands(cells["nsnd"], edges) + rows = [] + for index, label in enumerate(band_labels(edges)): + members = cells[band == index] + weighted = members[["source"]].assign( + w=members["w"], + **{name: members["w"] * members[name] for name in amount_columns}, + ) + by_source = weighted.groupby("source").sum() + weight = by_source["w"] + row = { + "band": label, + "rows": len(members), + "sources": len(by_source), + "source_ess": float(weight.sum() ** 2 / np.square(weight).sum()) + if len(weight) + else 0.0, + "count": float(weight.sum()), + "top_source_weight_share": float(weight.max() / weight.sum()) + if len(weight) + else 0.0, + } + for name in amount_columns: + amount = by_source[name] + row[name] = float(amount.sum()) + row[f"top_source_{name}_share"] = ( + float(amount.max() / amount.sum()) + if len(amount) and amount.sum() > 0 + else 0.0 + ) + rows.append(row) + return pd.DataFrame(rows) + + +def _multinomial_draws(n_sources: int, replicates: int, seed: int) -> np.ndarray: + rng = np.random.default_rng(seed) + return rng.multinomial( + n_sources, np.full(n_sources, 1.0 / n_sources), size=replicates + ) + + +def bootstrap( + cells: pd.DataFrame, + edges: Sequence[int], + amount_columns: Sequence[str], + *, + replicates: int = 1_000, + seed: int = 0, + level: float = 0.9, +) -> pd.DataFrame: + """Percentile intervals and coefficients of variation per band. + + Sources are resampled with replacement across all selected persons, so + band membership varies between replicates, and a source's clones move + together. Weights stay fixed, so calibration error is left out and the + intervals are too narrow. + """ + + codes, sources = pd.factorize(cells["source"]) + draws = _multinomial_draws(len(sources), replicates, seed) + band = assign_bands(cells["nsnd"], edges) + tail = (1.0 - level) / 2.0 + weights = cells["w"].to_numpy(dtype=float) + rows = [] + for index, label in enumerate(band_labels(edges)): + members = band == index + multiplier = draws[:, codes[members]] + row = {"band": label} + series = {"count": weights[members]} + for name in amount_columns: + series[name] = weights[members] * cells.loc[members, name].to_numpy( + dtype=float + ) + for name, values in series.items(): + replicate_totals = multiplier @ values + estimate = values.sum() + row[f"{name}_low"], row[f"{name}_high"] = np.quantile( + replicate_totals, [tail, 1 - tail] + ) + row[f"{name}_cv"] = ( + float(replicate_totals.std() / estimate) if estimate else float("nan") + ) + rows.append(row) + return pd.DataFrame(rows) + + +def bootstrap_total( + values: object, + sources: object, + *, + replicates: int = 1_000, + seed: int = 0, + level: float = 0.9, +) -> dict[str, float]: + """The same source bootstrap for one weighted total (``values`` already weighted).""" + + values = np.asarray(values, dtype=float) + codes, unique = pd.factorize(pd.Series(sources)) + draws = _multinomial_draws(len(unique), replicates, seed) + by_source = np.bincount(codes, weights=values, minlength=len(unique)) + totals = draws @ by_source + tail = (1.0 - level) / 2.0 + low, high = np.quantile(totals, [tail, 1 - tail]) + estimate = float(values.sum()) + return { + "estimate": estimate, + "low": float(low), + "high": float(high), + "cv": float(totals.std() / abs(estimate)) if estimate else float("nan"), + "sources": int(len(unique)), + "top_source_share": float(np.abs(by_source).max() / np.abs(by_source).sum()) + if len(unique) + else 0.0, + } + + +def _primary( + table: pd.DataFrame, + amount_columns: Sequence[str], + min_sources: int, + dominance: float, +) -> list[str]: + statuses = [] + for _, row in table.iterrows(): + shares = [row["top_source_weight_share"]] + [ + row[f"top_source_{name}_share"] for name in amount_columns + ] + if row["sources"] < min_sources: + statuses.append(FEWER_THAN_FLOOR) + elif max(shares) > dominance: + statuses.append(DOMINATED) + else: + statuses.append(PUBLISHED) + return statuses + + +def _passes( + cells: pd.DataFrame, + mask: np.ndarray, + amount_columns: Sequence[str], + min_sources: int, + dominance: float, +) -> bool: + """Whether the persons under ``mask``, taken as one cell, pass the primary rules.""" + + members = cells[mask] + weighted = members[["source"]].assign( + w=members["w"], + **{name: members["w"] * members[name] for name in amount_columns}, + ) + by_source = weighted.groupby("source").sum() + if len(by_source) < min_sources: + return False + for column in ("w", *amount_columns): + total = by_source[column].sum() + if total > 0 and by_source[column].max() / total > dominance: + return False + return True + + +def disclosure_statuses( + cells: pd.DataFrame, + fine_edges: Sequence[int], + coarse_edges: Sequence[int], + amount_columns: Sequence[str], + *, + total_published: bool = True, + min_sources: int = MIN_SOURCES, + dominance: float = DOMINANCE_SHARE, +) -> tuple[list[str], list[str]]: + """Statuses for Release 1 (the ``fine`` bands) and Release 2 (the ``coarse`` + bins) together, read as one hierarchy: total, bins, bands. + + Primary: fewer than ``min_sources`` sources, or one source above + ``dominance`` of the cell's weight or of any amount. + + Secondary: whatever a reader can rebuild by subtraction must itself pass + the primary rules, so cells are withheld, fewest sources first, until: + + - the withheld bands inside each published bin pass, taken together; + - no withheld bin can be rebuilt from its bands; + - with ``total_published``, the withheld bins pass together, and so do the + withheld bands inside them. The 1p costing is about one per cent of the + total property slice, so it gives the total away, and it is published. + """ + + fine = audit(cells, fine_edges, amount_columns) + coarse = audit(cells, coarse_edges, amount_columns) + fine_status = _primary(fine, amount_columns, min_sources, dominance) + coarse_status = _primary(coarse, amount_columns, min_sources, dominance) + band = assign_bands(cells["nsnd"], fine_edges) + holder = coarse_index_of_fine(fine_edges, coarse_edges) + + def passes(bands: Sequence[int]) -> bool: + mask = np.isin(band, list(bands)) + return _passes(cells, mask, amount_columns, min_sources, dominance) + + def fewest(candidates: Sequence[int], table: pd.DataFrame) -> int: + return min( + candidates, key=lambda i: (table.loc[i, "sources"], table.loc[i, "count"]) + ) + + def bands_in(bins: Sequence[int]) -> list[int]: + return [i for i in range(len(fine_edges)) if holder[i] in bins] + + def withhold_more() -> bool: + withheld_bins = [j for j, s in enumerate(coarse_status) if s != PUBLISHED] + published_bins = [j for j, s in enumerate(coarse_status) if s == PUBLISHED] + if ( + total_published + and withheld_bins + and published_bins + and not passes(bands_in(withheld_bins)) + ): + coarse_status[fewest(published_bins, coarse)] = SECONDARY + return True + for j in withheld_bins: + members = bands_in([j]) + if all(fine_status[i] == PUBLISHED for i in members): + fine_status[fewest(members, fine)] = SECONDARY + return True + groups = [bands_in([j]) for j in published_bins] + if total_published and withheld_bins: + groups.append(bands_in(withheld_bins)) + for members in groups: + withheld = [i for i in members if fine_status[i] != PUBLISHED] + published = [i for i in members if fine_status[i] == PUBLISHED] + if withheld and published and not passes(withheld): + fine_status[fewest(published, fine)] = SECONDARY + return True + return False + + while withhold_more(): + pass + return fine_status, coarse_status + + +def _round(value: float, step: float) -> float: + return float(np.round(value / step) * step) + + +def render( + table: pd.DataFrame, + intervals: pd.DataFrame, + statuses: Sequence[str], + amount_columns: Sequence[str], + *, + mode: str, + count_step: float = 100.0, + amount_step: float = 100_000.0, + min_sources: int = MIN_SOURCES, +) -> pd.DataFrame: + """A release table. + + ``restricted`` shows every band in full, with exact record counts and the + status as a note; it is for recipients covered by the licence. ``public`` + blanks withheld bands and shows record counts under the floor as ``<10``. + """ + + if mode not in {"restricted", "public"}: + raise ValueError(f"Unknown mode {mode!r}.") + rows = [] + for index, row in table.iterrows(): + interval = intervals.loc[index] + status = statuses[index] + shown = mode == "restricted" or status == PUBLISHED + out = {"band": row["band"]} + if mode == "public" and row["sources"] < min_sources: + out["distinct_records"] = f"<{min_sources}" + else: + out["distinct_records"] = str(int(row["sources"])) + out["effective_records"] = round(float(row["source_ess"]), 1) if shown else None + for name, step in [ + ("count", count_step), + *[(n, amount_step) for n in amount_columns], + ]: + for suffix, value in ( + ("", row[name]), + ("_low", interval[f"{name}_low"]), + ("_high", interval[f"{name}_high"]), + ): + out[f"{name}{suffix}"] = _round(value, step) if shown else None + out[f"{name}_cv"] = ( + round(float(interval[f"{name}_cv"]), 2) if shown else None + ) + out["status"] = status + rows.append(out) + return pd.DataFrame(rows) From 142acc507461cdcd54a02b8cb6df183262fb0dd3 Mon Sep 17 00:00:00 2001 From: Maria Juaristi Date: Wed, 21 Oct 2026 12:22:48 +0100 Subject: [PATCH 2/4] Address the #1327 review: hidden record counts, headline count, checks MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit From Vahid's review at c4cddcdc9: 1. Public tables show a band's record count only when the band is published, so a thin band's count can't be rebuilt from its bin. 2. The headline counts taxpayers with taxable property income; landlords wholly under the property allowance are counted beside it in the restricted and private tables. 3. The private audit carries each band before and after the reweight, and the release note says what is calibrated and what is a test. 4. The tool refuses policyengine-uk older than 2.134.4. 5. SPI-channel rows that keep their FRS person's values are keyed by that person; the manifest counts them and any FRS person also hosting tape records. 6. A reweight that doesn't converge, or misses a fact or region total by more than 0.5%, stops the run. 7. The bootstrap resamples FRS households and tape records across all Scottish persons. 8. The back-cast lists the inputs it leaves at the file's year, and the round-trip test compares every column. 9. Read-back facts are checked against HMRC's printed SPI 2023-24 tables (Tables 3.7 and 3.11, shipped as a CSV). 10. The relief variant must reproduce the engine at zero on the file. 11. The costing marks rates-only as the headline and the relief variant as a sensitivity; the release note says static, 2023-24, no response. 12. The release note footnotes PA as £12,570. 13. Two records over 85% also withhold a cell. 14. The manifest moves to private/, a release note goes with each version, and the text helper is shared. Co-Authored-By: Claude Opus 5.5 --- ...test_uk_scotland_property_income_engine.py | 23 +- .../uk/test_uk_scotland_property_income.py | 160 ++++++- tools/uk_scotland_property_income/__main__.py | 453 ++++++++++++++---- tools/uk_scotland_property_income/backcast.py | 41 +- tools/uk_scotland_property_income/registry.py | 37 ++ tools/uk_scotland_property_income/release.py | 280 +++++++---- .../spi_2023_24_published.csv | 57 +++ 7 files changed, 830 insertions(+), 221 deletions(-) create mode 100644 tools/uk_scotland_property_income/spi_2023_24_published.csv diff --git a/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py b/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py index acfa33a1e..97776d0ce 100644 --- a/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py +++ b/packages/microcosm-build/tests/engine/uk/test_uk_scotland_property_income_engine.py @@ -12,7 +12,9 @@ from policyengine_uk.data.economic_assumptions import extend_single_year_dataset from tools.uk_scotland_property_income.backcast import ( + NOT_INVERTED, backcast_one_year, + held_at_file_year, refuse_earlier_period, with_household_weights, ) @@ -34,6 +36,8 @@ def _dataset(year: int) -> UKSingleYearDataset: "property_income": [10_000.0, 14_000.0, 10_000.0], "property_finance_costs": [0.0, 3_000.0, 0.0], "savings_interest_income": [200.0, 0.0, 200.0], + # Unmapped: the forward step leaves it flat, and so does the back-cast. + "gift_aid": [100.0, 0.0, 0.0], } ), benunit=pd.DataFrame({"benunit_id": ids}), @@ -67,18 +71,29 @@ def test_the_backcast_inverts_the_engines_forward_step(): < released.person["employment_income"].iloc[0] ) + assert held_at_file_year(released, divided)["person"] == ["gift_aid"] + + # Every column the forward step touches comes back, not only the divided ones. forward = extend_single_year_dataset(earlier, parameters, end_year=2024).datasets[ 2024 ] - for table in ("person", "household"): + for table in ("person", "benunit", "household"): for column in getattr(released, table).columns: - if column in divided: + if column in NOT_INVERTED: + continue + after, before = ( + getattr(forward, table)[column], + getattr(released, table)[column], + ) + if pd.api.types.is_numeric_dtype(before.dtype): np.testing.assert_allclose( - getattr(forward, table)[column].to_numpy(dtype=float), - getattr(released, table)[column].to_numpy(dtype=float), + after.to_numpy(dtype=float), + before.to_numpy(dtype=float), rtol=1e-12, err_msg=column, ) + else: + assert after.tolist() == before.tolist(), column def test_an_earlier_period_is_refused(): diff --git a/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py b/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py index edc3c1836..037d9565d 100644 --- a/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py +++ b/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py @@ -185,23 +185,44 @@ def test_rake_stops_at_its_bounds_when_a_fact_is_out_of_reach(): assert result.after[0] == pytest.approx(20.0) -def test_source_keys_count_clones_and_copies_once(): +def _persons_table(): table = pd.DataFrame( { - "person_support_channel": [b"frs", b"frs", b"spi", b"spi", b"spi"], - "person_source_id": [1001, 1001, 2001, 2001, 2001], - **{name: [0.0] * 5 for name in release.SOURCE_INCOME_COLUMNS}, + "person_support_channel": [b"frs", b"frs", b"spi", b"spi", b"spi", b"spi"], + "person_source_id": [1001, 1001, 2001, 2001, 2001, 1001], + "person_is_spi_income_recipient": [0, 0, 1, 1, 1, 0], + **{name: [0.0] * 6 for name in release.SOURCE_INCOME_COLUMNS}, } ) + table.loc[[0, 1, 5], "property_income"] = 8_000.0 table.loc[[2, 3], "property_income"] = 12_000.0 table.loc[4, "property_income"] = 13_000.0 - keys = release.source_keys(table) + return table + + +def test_source_keys_count_clones_and_copies_once(): + keys = release.source_keys(_persons_table()) assert keys[0] == keys[1] == "frs:1001" - # SPI rows share their FRS host's id, so their incomes tell them apart. + # Tape rows share their FRS host's id, so their incomes tell them apart. assert keys[2] == keys[3] != keys[4] + # An SPI-channel row that keeps its FRS person's values is that person. + assert keys[5] == "frs:1001" + +def test_clusters_are_frs_households_and_tape_records(): + table = _persons_table() + clusters = release.cluster_keys(table, ["h1", "h1", "h2", "h2", "h3", "h4"]) + assert clusters[0] == clusters[1] == "frs-hh:h1" + assert clusters[5] == "frs-hh:h4" + assert clusters[2] == clusters[3] == release.source_keys(table)[2] -def _cells(sources_per_band, rows_per_source=3, band_values=(13_000, 16_000, 21_000)): + +def _cells( + sources_per_band, + rows_per_source=3, + band_values=(13_000, 16_000, 21_000), + weight=10.0, +): rows = [] for band_value, n_sources in zip(band_values, sources_per_band, strict=True): for source in range(n_sources): @@ -210,7 +231,8 @@ def _cells(sources_per_band, rows_per_source=3, band_values=(13_000, 16_000, 21_ { "nsnd": band_value, "source": f"{band_value}:{source}", - "w": 10.0, + "cluster": f"{band_value}:{source}", + "w": weight, "taxable_property_income": 1_000.0 + source, } ) @@ -226,13 +248,22 @@ def test_audit_counts_sources_not_rows(): assert table["source_ess"].tolist() == pytest.approx([4, 12, 20]) -def test_bootstrap_moves_clones_together(): +def test_bootstrap_moves_clones_together_and_lets_the_domain_vary(): cells = _cells((1, 0, 0), rows_per_source=5) - intervals = release.bootstrap( - cells, (12_570, 15_000, 20_000), AMOUNTS, replicates=50 + edges = (12_570, 15_000, 20_000) + alone = release.bootstrap( + cells, edges, AMOUNTS, universe=cells["cluster"], replicates=50 + ) + # A single cluster is drawn once in every replicate, so its band never varies. + assert alone.loc[0, "count_low"] == alone.loc[0, "count_high"] == 50.0 + # Among many non-landlord clusters, it is drawn 0, 1, 2… times. + universe = [*cells["cluster"].unique(), *[f"other:{i}" for i in range(200)]] + varying = release.bootstrap( + cells, edges, AMOUNTS, universe=universe, replicates=200 ) - # A single source is drawn once in every replicate, so its band never varies. - assert intervals.loc[0, "count_low"] == intervals.loc[0, "count_high"] == 50.0 + assert varying.loc[0, "count_low"] == 0.0 and varying.loc[0, "count_high"] >= 100.0 + with pytest.raises(ValueError, match="outside"): + release.bootstrap(cells, edges, AMOUNTS, universe=["other"], replicates=5) FEW, SECOND, PUBLISHED = release.FEWER_THAN_FLOOR, release.SECONDARY, release.PUBLISHED @@ -266,14 +297,51 @@ def test_one_dominant_record_withholds_its_band(): assert fine_status[1] == release.DOMINATED -def test_public_tables_blank_withheld_bands_and_floor_record_counts(): - edges = (12_570, 15_000, 20_000) +def test_two_dominant_records_withhold_their_band(): + cells = _cells((0, 12, 0)) + # Each of the two is under half, but together they are over 85%. + cells.loc[cells["source"].isin(["16000:0", "16000:1"]), "w"] = 300.0 + table = release.audit(cells, FINE_EDGES, AMOUNTS) + assert table.loc[1, "top_source_weight_share"] < release.DOMINANCE_SHARE + assert table.loc[1, "top_two_source_weight_share"] > release.TOP_TWO_SHARE + fine_status, _ = release.disclosure_statuses(cells, FINE_EDGES, (12_570,), AMOUNTS) + assert fine_status[1] == release.DOMINATED + + +def _public_tables(cells, fine_edges, coarse_edges): + statuses = release.disclosure_statuses(cells, fine_edges, coarse_edges, AMOUNTS) + tables = [] + for edges, status in zip((fine_edges, coarse_edges), statuses, strict=True): + table = release.audit(cells, edges, AMOUNTS) + intervals = release.bootstrap( + cells, edges, AMOUNTS, universe=cells["cluster"], replicates=20 + ) + tables.append(release.render(table, intervals, status, AMOUNTS, mode="public")) + return tables + + +def test_public_record_counts_do_not_rebuild_a_withheld_band(): + # Vahid's case on #1327: a thin band beside two bands of 30 in one bin. + cells = _cells((6, 30, 30), band_values=(13_000, 16_000, 21_000)) + fine, coarse = _public_tables(cells, (12_570, 15_000, 20_000), (12_570,)) + shown = pd.to_numeric(fine["distinct_records"], errors="coerce") + bin_records = int(coarse.loc[0, "distinct_records"]) + hidden = fine.loc[shown.isna(), "band"].tolist() + # Every band whose count is hidden shares the remainder with another one. + assert len(hidden) >= 2 + assert bin_records - shown.sum() == 36 + assert fine["distinct_records"].tolist()[:2] == [release.WITHHELD_RECORDS] * 2 + + +def test_public_tables_blank_withheld_bands_and_their_record_counts(): cells = _cells((4, 12, 20)) - table = release.audit(cells, edges, AMOUNTS) - intervals = release.bootstrap(cells, edges, AMOUNTS, replicates=20) - statuses = [release.FEWER_THAN_FLOOR, release.SECONDARY, release.PUBLISHED] + table = release.audit(cells, FINE_EDGES, AMOUNTS) + intervals = release.bootstrap( + cells, FINE_EDGES, AMOUNTS, universe=cells["cluster"], replicates=20 + ) + statuses = [FEW, SECOND, PUBLISHED] public = release.render(table, intervals, statuses, AMOUNTS, mode="public") - assert public["distinct_records"].tolist() == ["<10", "12", "20"] + assert public["distinct_records"].tolist() == ["withheld", "withheld", "20"] assert pd.isna(public.loc[0, "count"]) and pd.isna( public.loc[1, "taxable_property_income"] ) @@ -281,3 +349,57 @@ def test_public_tables_blank_withheld_bands_and_floor_record_counts(): restricted = release.render(table, intervals, statuses, AMOUNTS, mode="restricted") assert restricted["distinct_records"].tolist() == ["4", "12", "20"] assert restricted.loc[0, "count"] == 100.0 + + +def test_a_single_number_needs_records_and_no_dominant_ones(): + shares = release.record_shares([5.0, 5.0, 0.0, 90.0], ["a", "a", "b", "c"]) + assert shares == { + "records": 2, + "top_record_share": 0.9, + "top_two_record_share": 1.0, + } + assert release.passes_primary(12, [0.2], [0.3]) == PUBLISHED + assert release.passes_primary(9, [0.2], [0.3]) == FEW + assert release.passes_primary(12, [0.2], [0.9]) == release.DOMINATED + + +def test_read_back_facts_are_checked_against_the_printed_tables(tmp_path): + printed = tmp_path / "printed.csv" + printed.write_text( + "fact,table,published,unit\n" + "hmrc/property_income_income_band_20_000_to_30_000,3.7,3820,1000000\n" + "hmrc.spi_region.taxpayers_by_region_20000_30000@S92000003,3.11,841,1000\n" + ) + facts = [ + registry.SpiFact( + "hmrc/property_income_income_band_20_000_to_30_000", + "property_amount", + 20_000, + 30_000, + None, + 3_820.3e6, + ), + registry.SpiFact( + "hmrc.spi_region.taxpayers_by_region_20000_30000@S92000003", + "taxpayers", + 20_000, + 30_000, + "SCOTLAND", + 843_061.0, + ), + registry.SpiFact("not_printed", "income_tax", 0, 1, "SCOTLAND", 1.0), + ] + rows = registry.check_against_published(facts, printed) + assert [(row["fact"][:8], row["agrees"]) for row in rows] == [ + ("hmrc/pro", True), + ("hmrc.spi", False), + ] + + +def test_the_shipped_printed_tables_cover_the_reweighted_facts(): + lines = registry.PUBLISHED_PATH.read_text().splitlines()[1:] + names = {line.split(",")[0] for line in lines} + assert "hmrc/property_income_income_band_1_000_000_to_inf" in names + assert "hmrc/property_income_count_income_band_300_000_to_500_000" in names + assert "hmrc.spi_region.total_income_by_region_200000_plus@S92000003" in names + assert len(names) == len(lines) == 56 diff --git a/tools/uk_scotland_property_income/__main__.py b/tools/uk_scotland_property_income/__main__.py index 2b11a622b..f85b7dfed 100644 --- a/tools/uk_scotland_property_income/__main__.py +++ b/tools/uk_scotland_property_income/__main__.py @@ -7,16 +7,20 @@ The output directory holds three kinds of file: -- ``private/``: the per-band audit (rows, exact sources, shares), the reweight - report and the drift on every SPI 2023-24 fact. Licensed working files; - never shared. -- ``restricted/``: both releases in full, with exact record counts. Only for - recipients whose licence covers the FRS and the SPI tape. -- ``public/``: both releases with bands under the floor withheld and flagged, - and secondary suppression across the two. +- ``private/``: the per-band audit (rows, exact records, households, shares, + before and after the reweight), the reweight report, the drift on every SPI + 2023-24 fact, the check against HMRC's printed tables, the costing detail + and the run manifest. Licensed working files; never shared. +- ``restricted/``: both releases in full, with exact record counts, the + costing and a release note. Only for recipients whose licence covers the + FRS and the SPI tape. +- ``public/``: both releases with withheld cells blank, the costing and a + release note. Safe to send as a folder. Everything is derived from UK Data Service licensed microdata, so the tool -refuses an output directory inside a git working tree. +refuses an output directory inside a git working tree. That check needs git on +the PATH, and it can't see a cloud-synced folder: choose a local, unsynced +private directory. """ from __future__ import annotations @@ -31,6 +35,7 @@ import numpy as np import pandas as pd +from packaging.version import Version from scipy import sparse from tools.uk_scotland_property_income import ( @@ -48,9 +53,14 @@ ) REPO_ROOT = Path(__file__).resolve().parents[2] +#: The headline amount first; profit before the property allowance beside it. AMOUNTS = ("taxable_property_income", "property_income") -#: Pounds of float32 rounding allowed in a non-Scottish taxpayer's tax. -OUTSIDE_SCOTLAND_TOLERANCE = 0.05 +#: Pounds of float32 rounding allowed in a tax the reform shouldn't touch. +ROUNDING_TOLERANCE = 0.05 +#: The engine the tests passed on. Older ones lack the private uprating map +#: the back-cast reads, and before it ``net_income`` didn't deduct trade-loss +#: relief, which moves the NSND bands. +MINIMUM_ENGINE = Version("2.134.4") ENGINE_VARIABLES = ( "person_id", "pays_scottish_income_tax", @@ -90,6 +100,19 @@ def _refuse_working_tree(out: Path) -> None: ) +def _refuse_old_engine() -> dict[str, str]: + versions = { + name: importlib.metadata.version(name) + for name in ("policyengine-uk", "policyengine-core") + } + if Version(versions["policyengine-uk"]) < MINIMUM_ENGINE: + raise SystemExit( + f"policyengine-uk {versions['policyengine-uk']} is older than " + f"{MINIMUM_ENGINE}, the engine this tool was tested on." + ) + return versions + + def _git_state() -> dict[str, Any]: def git(*args: str) -> str: return subprocess.run( @@ -104,12 +127,6 @@ def git(*args: str) -> str: } -def _text(values: pd.Series) -> np.ndarray: - return values.map( - lambda value: value.decode() if isinstance(value, bytes) else str(value) - ).to_numpy() - - def person_frame(simulation: Any, dataset: Any, year: int) -> pd.DataFrame: """The person columns every step reads, at ``year``.""" @@ -132,7 +149,9 @@ def person_frame(simulation: Any, dataset: Any, year: int) -> pd.DataFrame: {name: values[name] for name in ENGINE_VARIABLES if name != "person_id"} ) frame.insert(0, "household_row", household_row) - frame.insert(1, "region", _text(household["region"])[household_row]) + frame.insert( + 1, "region", release.as_text(household["region"]).to_numpy()[household_row] + ) frame["scottish"] = frame.pop("pays_scottish_income_tax") > 0 frame["nsnd"] = ( frame["net_income"] @@ -150,6 +169,77 @@ def _write(frame: pd.DataFrame, path: Path) -> None: frame.to_csv(path, index=False) +def _record_keys(released: Any) -> tuple[np.ndarray, np.ndarray, pd.DataFrame]: + """Records and bootstrap clusters for every person in the released file.""" + + person = released.person + household = released.household + key_by_household = pd.Series( + release.as_text(household["source_household_key"]).to_numpy(), + index=household["household_id"].to_numpy(), + ) + household_keys = key_by_household.loc[person["person_household_id"]].to_numpy() + return ( + release.source_keys(person), + release.cluster_keys(person, household_keys), + pd.DataFrame( + { + "person_source_id": person["person_source_id"] + .astype(np.int64) + .to_numpy(), + "channel": release.as_text(person["person_support_channel"]).to_numpy(), + } + ), + ) + + +def _reweight( + persons: pd.DataFrame, + earlier: Any, + facts: list[registry.SpiFact], + args: argparse.Namespace, +) -> tuple[rake.RakeResult, list[registry.SpiFact], list[str], dict[str, float]]: + raked_facts = [fact for fact in facts if fact.region is None] + [ + fact + for fact in facts + if fact.region == "SCOTLAND" and fact.measure in {"taxpayers", "total_income"} + ] + w0 = earlier.household["household_weight"].to_numpy(dtype=float) + household_regions = release.as_text(earlier.household["region"]).to_numpy() + regions = sorted(set(household_regions)) + matrix = sparse.hstack( + [ + measures.household_matrix(persons, raked_facts, len(w0)), + measures.region_matrix(regions, household_regions), + ], + format="csr", + ) + targets = np.concatenate( + [ + [fact.value for fact in raked_facts], + [w0[household_regions == region].sum() for region in regions], + ] + ) + raked = rake.rake( + matrix, w0, targets, lower=args.ratio_bounds[0], upper=args.ratio_bounds[1] + ) + errors = np.abs(raked.relative_error_after) + misses = { + "facts": float(errors[: len(raked_facts)].max()), + "region_households": float(errors[len(raked_facts) :].max()), + } + if not raked.converged or max(misses.values()) > args.max_miss: + raise SystemExit( + f"The reweight didn't settle: converged={raked.converged}, " + f"largest misses {misses}, allowed {args.max_miss}." + ) + return raked, raked_facts, regions, misses + + +def _with_weights(cells: pd.DataFrame, weights: np.ndarray) -> pd.DataFrame: + return cells.assign(w=weights[cells["household_row"].to_numpy()]) + + def run(args: argparse.Namespace) -> dict[str, Any]: from policyengine_uk import CountryTaxBenefitSystem, Microsimulation from policyengine_uk.data import UKSingleYearDataset @@ -158,6 +248,7 @@ def run(args: argparse.Namespace) -> dict[str, Any]: out = Path(args.out) _refuse_working_tree(out) + engine = _refuse_old_engine() h5, registry_path = Path(args.h5), Path(args.registry) h5_sha = _sha256(h5) if args.expect_sha256 and not h5_sha.startswith(args.expect_sha256): @@ -168,7 +259,7 @@ def run(args: argparse.Namespace) -> dict[str, Any]: raise SystemExit( f"{h5} is for {released.time_period}; this back-casts one year to {args.year}." ) - sources = release.source_keys(released.person) + sources, clusters, identities = _record_keys(released) parameters = CountryTaxBenefitSystem().parameters earlier, divided = backcast.backcast_one_year(released, parameters) w0 = earlier.household["household_weight"].to_numpy(dtype=float) @@ -178,35 +269,22 @@ def run(args: argparse.Namespace) -> dict[str, Any]: args.year, ) persons["source"] = sources + persons["cluster"] = clusters + persons["person_source_id"] = identities["person_source_id"] + persons["channel"] = identities["channel"] - # A4: reweight to the SPI 2023-24 facts, holding each region's households. + # A4: reweight to the SPI 2023-24 facts, checked against HMRC's tables. facts = registry.spi_2023_facts(registry.load_registry(registry_path)) - raked_facts = [fact for fact in facts if fact.region is None] + [ - fact - for fact in facts - if fact.region == "SCOTLAND" and fact.measure in {"taxpayers", "total_income"} - ] - household_regions = _text(earlier.household["region"]) - regions = sorted(set(household_regions)) - n_households = len(w0) - matrix = sparse.hstack( - [ - measures.household_matrix(persons, raked_facts, n_households), - measures.region_matrix(regions, household_regions), - ], - format="csr", - ) - targets = np.concatenate( - [ - [fact.value for fact in raked_facts], - [w0[household_regions == region].sum() for region in regions], - ] - ) - raked = rake.rake( - matrix, w0, targets, lower=args.ratio_bounds[0], upper=args.ratio_bounds[1] - ) + raked, raked_facts, regions, misses = _reweight(persons, earlier, facts, args) + published_check = pd.DataFrame(registry.check_against_published(raked_facts)) + _write(published_check, out / "private" / "published_check.csv") + if len(published_check) != len(raked_facts) or not published_check["agrees"].all(): + raise SystemExit( + "Some reweighted facts don't match HMRC's printed SPI 2023-24 tables; " + "see private/published_check.csv." + ) w1 = raked.weights - all_matrix = measures.household_matrix(persons, facts, n_households) + all_matrix = measures.household_matrix(persons, facts, len(w0)) raked_names = {fact.name for fact in raked_facts} drift = pd.DataFrame( { @@ -223,26 +301,21 @@ def run(args: argparse.Namespace) -> dict[str, Any]: ) drift["backcast_error"] = drift["backcast"] / drift["published_2023_24"] - 1 drift["reweighted_error"] = drift["reweighted"] / drift["published_2023_24"] - 1 + _write(drift, out / "private" / "drift.csv") - # A6: the two releases. + # A6: the two releases. The headline counts taxpayers with taxable property + # income; landlords wholly under the property allowance are counted beside + # it, in the restricted and private tables only. persons["w"] = w1[persons["household_row"].to_numpy()] - cells = persons[ + scottish = persons["scottish"].to_numpy() + landlords = persons[ persons["scottish"] & (persons["income_tax"] > 0) & (persons["property_income"] > 0) & (persons["nsnd"] >= PERSONAL_ALLOWANCE_2023) ] - tables = {} - for key, edges in ( - ("release_1", REQUESTED_EDGES), - ("release_2", INTENTIONAL_EDGES), - ): - tables[key] = ( - release.audit(cells, edges, AMOUNTS), - release.bootstrap( - cells, edges, AMOUNTS, replicates=args.replicates, seed=args.seed - ), - ) + cells = landlords[landlords["taxable_property_income"] > 0] + universe = persons.loc[scottish, "cluster"].to_numpy() statuses = dict( zip( ("release_1", "release_2"), @@ -252,63 +325,112 @@ def run(args: argparse.Namespace) -> dict[str, Any]: strict=True, ) ) - for key, (table, intervals) in tables.items(): - audit = table.merge(intervals, on="band").assign(status=statuses[key]) + for key, edges in ( + ("release_1", REQUESTED_EDGES), + ("release_2", INTENTIONAL_EDGES), + ): + table = release.audit(cells, edges, AMOUNTS) + intervals = release.bootstrap( + cells, + edges, + AMOUNTS, + universe=universe, + replicates=args.replicates, + seed=args.seed, + ) + before = release.audit(_with_weights(cells, w0), edges, AMOUNTS) + everyone = release.audit(landlords, edges, ()) + audit = table.merge(intervals, on="band").assign( + status=statuses[key], + count_before_reweight=before["count"], + **{f"{n}_before_reweight": before[n] for n in AMOUNTS}, + all_landlords_count=everyone["count"], + all_landlords_records=everyone["sources"], + ) _write(audit, out / "private" / f"audit_{key}.csv") - for mode in ("restricted", "public"): - _write( - release.render(table, intervals, statuses[key], AMOUNTS, mode=mode), - out / mode / f"{key}.csv", - ) + restricted = release.render( + table, intervals, statuses[key], AMOUNTS, mode="restricted" + ) + restricted.insert( + restricted.columns.get_loc("status"), + "all_landlords_count", + (everyone["count"] / 100).round() * 100, + ) + _write(restricted, out / "restricted" / f"{key}.csv") + _write( + release.render(table, intervals, statuses[key], AMOUNTS, mode="public"), + out / "public" / f"{key}.csv", + ) - # B: the illustrative 1p, with and without the finance-cost reduction following. - scottish = persons["scottish"].to_numpy() - costing = {} - for variant, relief_follows in (("rates_only", False), ("relief_follows", True)): + # B: the illustrative 1p. The relief variant copies the engine's formula, + # so check first that at zero it reproduces the engine on this file. + def tax_change(delta: float, relief_follows: bool) -> np.ndarray: reformed = Microsimulation( dataset=backcast.with_household_weights(earlier, w0), - reform=scottish_property_rate_reform( - args.delta, relief_follows=relief_follows - ), + reform=scottish_property_rate_reform(delta, relief_follows=relief_follows), ) - change = ( + return ( np.asarray(reformed.calculate("income_tax", args.year), dtype=float) - persons["income_tax"].to_numpy() ) - # The engine stores float32, and the relief variant's rate is a float64 - # array, so tax outside Scotland can move by a rounding penny. - if np.abs(change[~scottish]).max(initial=0.0) > OUTSIDE_SCOTLAND_TOLERANCE: - raise ValueError("The Scottish reform changed tax outside Scotland.") + + if np.abs(tax_change(0.0, True)).max() > ROUNDING_TOLERANCE: + raise SystemExit("The relief variant at zero doesn't reproduce the engine.") + costing = {} + for variant, relief_follows, role in ( + ("rates_only", False, "headline"), + ("relief_follows", True, "sensitivity"), + ): + change = tax_change(args.delta, relief_follows) + if np.abs(change[~scottish]).max(initial=0.0) > ROUNDING_TOLERANCE: + raise SystemExit("The Scottish reform changed tax outside Scotland.") change = np.where(scottish, change, 0.0) weighted = persons["w"].to_numpy() * change summary = release.bootstrap_total( weighted[scottish], - persons.loc[scottish, "source"], + persons.loc[scottish, "cluster"], replicates=args.replicates, seed=args.seed, ) + shares = release.record_shares( + weighted[scottish], persons.loc[scottish, "source"] + ) affected = scottish & (change > 0) - summary["taxpayers_affected"] = float(persons.loc[affected, "w"].sum()) - summary["sources_affected"] = int(persons.loc[affected, "source"].nunique()) + summary.update( + shares, + role=role, + taxpayers_affected=float(persons.loc[affected, "w"].sum()), + status=release.passes_primary( + shares["records"], + [shares["top_record_share"]], + [shares["top_two_record_share"]], + ), + ) costing[variant] = summary _write_costing(costing, out) + spi_channel = (cells["channel"] == "spi") & cells["source"].str.startswith("frs:") manifest = { "issue": "PolicyEngine/microcosm#1326", "year": args.year, "inputs": { - "h5": {"path": str(h5.resolve()), "sha256": h5_sha}, + "h5": { + "path": str(h5.resolve()), + "sha256": h5_sha, + "tag": args.release_tag, + }, "registry": { "path": str(registry_path.resolve()), "sha256": _sha256(registry_path), }, }, - "engine": { - name: importlib.metadata.version(name) - for name in ("policyengine-uk", "policyengine-core") - }, + "engine": engine, "tool": _git_state(), - "backcast": {"divided": divided, "not_inverted": list(backcast.NOT_INVERTED)}, + "backcast": { + "divided": divided, + "not_inverted": list(backcast.NOT_INVERTED), + "held_at_file_year": backcast.held_at_file_year(released, divided), + }, "reweight": { "facts": [fact.name for fact in raked_facts], "region_controls": regions, @@ -316,7 +438,7 @@ def run(args: argparse.Namespace) -> dict[str, Any]: "converged": raked.converged, "message": raked.message, "max_abs_error_before": float(np.abs(raked.relative_error_before).max()), - "max_abs_error_after": float(np.abs(raked.relative_error_after).max()), + "max_abs_error_after": misses, "ratio_quantiles": dict( zip( ("min", "p01", "median", "p99", "max"), @@ -331,6 +453,30 @@ def run(args: argparse.Namespace) -> dict[str, Any]: ) ), "effective_sample_size_ratio": raked.effective_sample_size_ratio, + "published_check": { + "facts": int(len(published_check)), + "agree": int(published_check["agrees"].sum()), + }, + }, + "records": { + "headline_rows": int(len(cells)), + "headline_records": int(cells["source"].nunique()), + "headline_households_and_tape_records": int(cells["cluster"].nunique()), + "rows_keyed_by_frs_person": int(spi_channel.sum()), + "frs_persons_also_hosting_tape_records": int( + len( + set( + cells.loc[ + cells["source"].str.startswith("frs:"), "person_source_id" + ] + ) + & set( + cells.loc[ + cells["source"].str.startswith("spi:"), "person_source_id" + ] + ) + ) + ), }, "bands": { "release_1": list(REQUESTED_EDGES), @@ -340,34 +486,34 @@ def run(args: argparse.Namespace) -> dict[str, Any]: "bootstrap": { "replicates": args.replicates, "seed": args.seed, + "clusters": "FRS households and tape records, over all Scottish persons", "weights": "fixed", }, "delta": args.delta, } - (out / "private").mkdir(parents=True, exist_ok=True) - drift.to_csv(out / "private" / "drift.csv", index=False) - (out / "manifest.json").write_text(json.dumps(manifest, indent=2) + "\n") + (out / "private" / "manifest.json").write_text( + json.dumps(manifest, indent=2) + "\n" + ) + for mode in ("restricted", "public"): + (out / mode / "release_note.md").write_text(release_note(manifest, mode)) return manifest -def _write_costing(costing: dict[str, dict[str, float]], out: Path) -> None: +def _write_costing(costing: dict[str, dict[str, Any]], out: Path) -> None: (out / "private").mkdir(parents=True, exist_ok=True) (out / "private" / "costing.json").write_text(json.dumps(costing, indent=2) + "\n") rows = [] for variant, summary in costing.items(): - published = ( - summary["sources"] >= release.MIN_SOURCES - and summary["top_source_share"] <= release.DOMINANCE_SHARE - ) rows.append( { "variant": variant, + "role": summary["role"], "change_in_scottish_income_tax": round(summary["estimate"] / 1e5) * 1e5, "low": round(summary["low"] / 1e5) * 1e5, "high": round(summary["high"] / 1e5) * 1e5, "cv": round(summary["cv"], 2), "taxpayers_affected": round(summary["taxpayers_affected"] / 100) * 100, - "status": release.PUBLISHED if published else release.DOMINATED, + "status": summary["status"], } ) table = pd.DataFrame(rows) @@ -381,6 +527,110 @@ def _write_costing(costing: dict[str, dict[str, float]], out: Path) -> None: _write(public, out / "public" / "costing.csv") +def release_note(manifest: dict[str, Any], mode: str) -> str: + """The note that goes with the tables. It states definitions, method and + caveats, and holds no figure from the tables.""" + + engine = manifest["engine"]["policyengine-uk"] + source = manifest["inputs"]["h5"] + file_label = source["tag"] or f"sha256 {source['sha256'][:8]}…" + region_miss = manifest["reweight"]["max_abs_error_after"]["region_households"] + starts = manifest["bands"]["scottish_rate_starts_2023"] + version = ( + "**Restricted version.** Every band is shown, with exact record counts. " + "It is only for recipients covered by the UK Data Service licence for " + "the Family Resources Survey and the SPI public use tape. Don't publish it." + if mode == "restricted" + else "**Public version.** Cells that fail the disclosure rules below are " + "blank, with no record count." + ) + all_landlords = ( + " The restricted tables also count them, in `all_landlords_count`." + if mode == "restricted" + else "" + ) + return f"""# Scottish taxpayers with property income, by NSND income band, 2023-24 + +Prepared by PolicyEngine from the Microcosm UK national file {file_label} +(PolicyEngine/microcosm#1326), with policyengine-uk {engine}. + +{version} + +## What the tables count + +- **Taxpayers:** people who pay Scottish income tax and have a positive income + tax liability in 2023-24, after tax reductions such as the finance-cost + reduction. +- **Property income:** the profit of a property business after allowable + expenses and before finance costs, the concept the SPI uses. The headline + amount, `taxable_property_income`, is after the £1,000 property allowance. + `property_income`, before the allowance, is shown beside it. +- **Headline count:** taxpayers with taxable property income above zero. + Landlords whose whole profit is under the property allowance are left out, + since most of them don't report it to HMRC.{all_landlords} +- **Bands:** individual non-savings, non-dividend (NSND) income before + allowances, after trade-loss relief and the property allowance. "PA" means + £12,570 for everyone, so marriage allowance transferors and blind person's + allowance holders are banded on £12,570 too. Taxpayers with NSND income + below £12,570 are not in the tables. +- **Release 1** keeps the ten requested bands. **Release 2** uses four unions + of them. Each edge is the requested edge nearest a 2023-24 Scottish rate + start: basic £{starts["basic"]:,}, intermediate £{starts["intermediate"]:,} + and higher £{starts["higher"]:,}. The bins are PA–£24,999 (starter and + basic), £25,000–£39,999 (intermediate), £40,000–£89,999 (higher) and + £90,000+ (higher and top). + +## How the 2023-24 figures are made + +- **Base:** the 2024-25 Microcosm file, back-cast one year. Each input is + divided by the growth index policyengine-uk uses to project it forward: + OBR national aggregates per person, and ONS population for the weights. + Inputs the engine doesn't uprate, such as Gift Aid and losses brought + forward, stay at their 2024-25 levels. +- **Reweight:** the back-cast weights are reweighted to the SPI 2023-24 tables + HMRC published. The targets are UK property income by range of total income + (counts and amounts, Table 3.7) and Scotland's taxpayers and total income by + range (Table 3.11). Each region's household total is held to within + {region_miss:.2%}. Because those figures are calibrated, they don't test the + data. The split of Scottish property income across NSND bands is not + targeted, so that split is what the outturn tests. +- **Law:** 2023-24 income tax. + +## Uncertainty + +Each band shows: +- the number of distinct survey records behind it; +- its effective number of records; +- a 90% interval, from resampling FRS households and SPI records across all + Scottish taxpayers. + +The weights are held fixed, so the intervals leave out calibration error and +are too narrow. + +## The 1p costing + +- **What it adds:** 1p on every Scottish rate charged on property income, + meaning the slice of property income above the allowances, with the same + band limits as the other NSND rates. +- **Basis:** static, on 2023-24 data and law, with no behavioural response. +- **Headline and sensitivity:** the headline keeps the finance-cost reduction + at the basic rate. The sensitivity raises it by 1p as well. +- **Later years:** a 2027-28 figure, when the power applies, needs a separate + run. + +## Disclosure rules + +A cell is shown only if: +- at least 10 distinct records sit behind it; +- no single record is over 50% of its weight or amounts; +- no two records together are over 85% of them. + +Further cells are withheld so that none of these can be rebuilt by +subtraction. That includes subtraction from the costing, which reveals the +total. +""" + + def parse_args(argv: list[str] | None = None) -> argparse.Namespace: parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter @@ -403,7 +653,16 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace: parser.add_argument( "--expect-sha256", help="Refuse an H5 whose sha256 doesn't start with this." ) + parser.add_argument( + "--release-tag", help="The release's cut tag, for the release note." + ) parser.add_argument("--ratio-bounds", type=float, nargs=2, default=(0.25, 4.0)) + parser.add_argument( + "--max-miss", + type=float, + default=0.005, + help="Refuse a reweight that misses any fact or region total by more (default 0.5%%).", + ) parser.add_argument("--replicates", type=int, default=1_000) parser.add_argument("--seed", type=int, default=0) parser.add_argument( @@ -419,6 +678,6 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace: manifest = run(parse_args()) reweight = manifest["reweight"] print( - f"done: reweight converged={reweight['converged']}, " - f"max |error| {reweight['max_abs_error_before']:.3f} -> {reweight['max_abs_error_after']:.4f}" + f"done: reweight converged={reweight['converged']}, largest miss " + f"{reweight['max_abs_error_before']:.3f} -> {reweight['max_abs_error_after']}" ) diff --git a/tools/uk_scotland_property_income/backcast.py b/tools/uk_scotland_property_income/backcast.py index ef0491195..fd794ed95 100644 --- a/tools/uk_scotland_property_income/backcast.py +++ b/tools/uk_scotland_property_income/backcast.py @@ -7,15 +7,26 @@ ``apply_single_year_uprating`` would multiply it by. Household weights move with ONS population growth, as they do going forward. -The forward step's custom pieces are not inverted. Council tax and rent move -by region, and student-loan plans are re-labelled; none of them enters income -tax, and they stay at the file's own year. +What the back-cast leaves at the file's own year: + +- the forward step's custom pieces. Council tax and rent move by region, and + student-loan plans are re-labelled. None of them enters income tax. +- every stored input that ``uprating_indices.yaml`` doesn't map, which the + forward step leaves flat too. Some feed income tax, such as ``gift_aid``, + ``loss_relief``, ``capital_allowances``, ``employment_expenses`` and + salary-sacrificed pension contributions. They are small, but not zero. + ``held_at_file_year`` lists them for the run's manifest. + +The indices are OBR national aggregates per person (per-capita GDP for +property income, mortgage interest for finance costs, ONS population for +weights). Differences between income bands are left to the reweight. """ from __future__ import annotations from typing import Any +import pandas as pd from policyengine_uk.data import UKSingleYearDataset from policyengine_uk.data.economic_assumptions import ( LEGACY_PROJECTION_END_YEAR, @@ -84,3 +95,27 @@ def with_household_weights( household=household, fiscal_year=int(dataset.time_period), ) + + +#: Identifier and weight columns that carry no money amount. +_NOT_AMOUNTS = ("_id", "_weight", "_index", "_build_id", "_source_id", "_key") + + +def held_at_file_year( + dataset: UKSingleYearDataset, divided: dict[str, Any] +) -> dict[str, list[str]]: + """The stored numeric columns, by table, that the back-cast leaves at the + file's own year: neither divided nor among the custom pieces.""" + + held = {} + for name in ("person", "benunit", "household"): + table = getattr(dataset, name) + held[name] = sorted( + column + for column in table.columns + if pd.api.types.is_float_dtype(table[column].dtype) + and column not in divided + and column not in NOT_INVERTED + and not column.endswith(_NOT_AMOUNTS) + ) + return held diff --git a/tools/uk_scotland_property_income/registry.py b/tools/uk_scotland_property_income/registry.py index e6c5a2560..a335dce52 100644 --- a/tools/uk_scotland_property_income/registry.py +++ b/tools/uk_scotland_property_income/registry.py @@ -99,3 +99,40 @@ def spi_2023_facts(registry: Mapping[str, Any]) -> list[SpiFact]: ) ) return facts + + +#: SPI 2023-24 as HMRC published it: Table 3.7 (property income by range of +#: total income) and Table 3.11 (Scotland), from Collated_Tables_3_1_to_3_11_2324.ods +#: (https://www.gov.uk/government/statistics/personal-incomes-statistics-for-the-tax-year-2023-to-2024, +#: sha256 ad063b06…), in thousands and £ millions as printed. +PUBLISHED_PATH = Path(__file__).with_name("spi_2023_24_published.csv") + + +def check_against_published( + facts: list[SpiFact], path: str | Path = PUBLISHED_PATH +) -> list[dict[str, Any]]: + """Each fact read back from the registry beside HMRC's printed figure. + + A fact agrees when it is within half the printed unit, the rounding of + the published table. Facts with no printed figure here are skipped. + """ + + published = {} + for line in Path(path).read_text().splitlines()[1:]: + fact, table, value, unit = line.split(",") + published[fact] = (table, float(value) * float(unit), float(unit)) + rows = [] + for fact in facts: + if fact.name not in published: + continue + table, value, unit = published[fact.name] + rows.append( + { + "fact": fact.name, + "table": table, + "read_back": fact.value, + "published": value, + "agrees": abs(fact.value - value) <= unit / 2, + } + ) + return rows diff --git a/tools/uk_scotland_property_income/release.py b/tools/uk_scotland_property_income/release.py index efef22461..ed2c29004 100644 --- a/tools/uk_scotland_property_income/release.py +++ b/tools/uk_scotland_property_income/release.py @@ -2,8 +2,9 @@ Rows repeat records: CGT clones, support copies and, in a dense file, area clones all carry one source's values. So every count here is of distinct -sources, every interval resamples sources with their clones together, and the -disclosure rules read sources, never rows. +records, the disclosure rules read records, never rows, and the intervals +resample whole clusters: an FRS household, whose members often share a let +property, or a tape record. """ from __future__ import annotations @@ -19,8 +20,8 @@ coarse_index_of_fine, ) -#: The stored incomes that identify an SPI-channel record. Its tape record isn't -#: carried (microcosm#1287), so records with the same income vector count once. +#: The stored incomes that identify a tape record. Its id isn't carried +#: (microcosm#1287), so tape rows with the same income vector count once. SOURCE_INCOME_COLUMNS = ( "property_income", "employment_income", @@ -29,110 +30,171 @@ "savings_interest_income", "dividend_income", ) +#: The income stage's flag for rows whose incomes come from the tape +#: (``uk_runtime/spi_income.py``). Every other row, in either channel, keeps +#: its FRS person's values (``uk_runtime/property_components.py``). +TAPE_RECIPIENT_COLUMN = "person_is_spi_income_recipient" MIN_SOURCES = 10 DOMINANCE_SHARE = 0.5 +TOP_TWO_SHARE = 0.85 PUBLISHED = "published" FEWER_THAN_FLOOR = "withheld: fewer than 10 records" -DOMINATED = "withheld: one record dominates" +DOMINATED = "withheld: one or two records dominate" SECONDARY = "withheld: protects another band" +WITHHELD_RECORDS = "withheld" -def _text(values: pd.Series) -> pd.Series: +def as_text(values: pd.Series) -> pd.Series: + """Strings, whether the store kept them as bytes or not.""" + return values.map( lambda value: value.decode() if isinstance(value, bytes) else str(value) ) +def _flag(values: pd.Series) -> np.ndarray: + if values.dtype == object: + return as_text(values).isin({"True", "true", "1"}).to_numpy() + return values.to_numpy(dtype=bool) + + +def _keeps_frs_values(person_table: pd.DataFrame) -> np.ndarray: + if TAPE_RECIPIENT_COLUMN not in person_table: + raise ValueError(f"The person table has no {TAPE_RECIPIENT_COLUMN!r}.") + channel = as_text(person_table["person_support_channel"]).to_numpy() + return (channel == "frs") | ~_flag(person_table[TAPE_RECIPIENT_COLUMN]) + + +def _tape_keys(person_table: pd.DataFrame) -> np.ndarray: + incomes = ( + person_table[list(SOURCE_INCOME_COLUMNS)].astype(float).round(2).astype(str) + ) + return ("spi:" + incomes.agg("|".join, axis=1)).to_numpy() + + def source_keys(person_table: pd.DataFrame) -> np.ndarray: """One key per independent record. - FRS persons are keyed by ``person_source_id``. SPI-channel persons carry - their FRS host's id, not their tape record's, so they are keyed by their - stored income vector. Read it from the released file, before any + A row that keeps FRS values is keyed by ``person_source_id``. That + includes the SPI channel's dependants and twin copies with no taxable + income, which carry their FRS person's values. A row with tape incomes + carries its FRS host's id, not its tape record's, so it is keyed by its + stored income vector. Read the table from the released file, before any back-cast. """ - channel = _text(person_table["person_support_channel"]).to_numpy() frs = "frs:" + person_table["person_source_id"].astype(np.int64).astype(str) - incomes = ( - person_table[list(SOURCE_INCOME_COLUMNS)].astype(float).round(2).astype(str) + return np.where( + _keeps_frs_values(person_table), frs.to_numpy(), _tape_keys(person_table) ) - spi = "spi:" + incomes.agg("|".join, axis=1) - return np.where(channel == "frs", frs.to_numpy(), spi.to_numpy()) + + +def cluster_keys( + person_table: pd.DataFrame, source_household_keys: object +) -> np.ndarray: + """The bootstrap's resampling unit for each person: the FRS household for + rows with FRS values, so partners who share a let property move together, + and the tape record for rows with tape incomes.""" + + households = "frs-hh:" + pd.Series(source_household_keys).astype(str) + return np.where( + _keeps_frs_values(person_table), households.to_numpy(), _tape_keys(person_table) + ) + + +def _top_shares(values: pd.Series) -> tuple[float, float]: + total = values.sum() + if len(values) == 0 or total <= 0: + return 0.0, 0.0 + largest = np.sort(values.to_numpy(dtype=float))[::-1] + return float(largest[0] / total), float(largest[:2].sum() / total) + + +def _by_source(members: pd.DataFrame, amount_columns: Sequence[str]) -> pd.DataFrame: + weighted = members[["source"]].assign( + w=members["w"], + **{name: members["w"] * members[name] for name in amount_columns}, + ) + return weighted.groupby("source").sum() def audit( cells: pd.DataFrame, edges: Sequence[int], amount_columns: Sequence[str] ) -> pd.DataFrame: - """Per band: rows, sources, effective sample size by source, weighted totals - and the largest source's share of the weight and of each amount. + """Per band: rows, records, households, effective sample size by record, + weighted totals, and the largest record's and the largest two records' + shares of the weight and of each amount. - ``cells`` has one row per selected person, with ``w``, ``nsnd``, ``source`` - and the amount columns. + ``cells`` has one row per selected person, with ``w``, ``nsnd``, ``source``, + ``cluster`` and the amount columns. """ band = assign_bands(cells["nsnd"], edges) rows = [] for index, label in enumerate(band_labels(edges)): members = cells[band == index] - weighted = members[["source"]].assign( - w=members["w"], - **{name: members["w"] * members[name] for name in amount_columns}, - ) - by_source = weighted.groupby("source").sum() + by_source = _by_source(members, amount_columns) weight = by_source["w"] row = { "band": label, "rows": len(members), "sources": len(by_source), + "households": members["cluster"].nunique(), "source_ess": float(weight.sum() ** 2 / np.square(weight).sum()) if len(weight) else 0.0, "count": float(weight.sum()), - "top_source_weight_share": float(weight.max() / weight.sum()) - if len(weight) - else 0.0, } + row["top_source_weight_share"], row["top_two_source_weight_share"] = ( + _top_shares(weight) + ) for name in amount_columns: amount = by_source[name] row[name] = float(amount.sum()) - row[f"top_source_{name}_share"] = ( - float(amount.max() / amount.sum()) - if len(amount) and amount.sum() > 0 - else 0.0 + row[f"top_source_{name}_share"], row[f"top_two_source_{name}_share"] = ( + _top_shares(amount) ) rows.append(row) return pd.DataFrame(rows) -def _multinomial_draws(n_sources: int, replicates: int, seed: int) -> np.ndarray: +def _cluster_draws(n_clusters: int, replicates: int, seed: int) -> np.ndarray: rng = np.random.default_rng(seed) return rng.multinomial( - n_sources, np.full(n_sources, 1.0 / n_sources), size=replicates + n_clusters, np.full(n_clusters, 1.0 / n_clusters), size=replicates ) +def _codes(clusters: object, universe: pd.Index) -> np.ndarray: + codes = universe.get_indexer(pd.Index(clusters)) + if np.any(codes < 0): + raise ValueError("Some clusters are outside the bootstrap's universe.") + return codes + + def bootstrap( cells: pd.DataFrame, edges: Sequence[int], amount_columns: Sequence[str], *, + universe: object, replicates: int = 1_000, seed: int = 0, level: float = 0.9, ) -> pd.DataFrame: """Percentile intervals and coefficients of variation per band. - Sources are resampled with replacement across all selected persons, so - band membership varies between replicates, and a source's clones move - together. Weights stay fixed, so calibration error is left out and the - intervals are too narrow. + Clusters are drawn with replacement from ``universe``, every cluster in + the domain (all Scottish persons), so the number of landlords varies + between replicates as well as which band they fall in. Weights stay + fixed, so calibration error is left out and the intervals are too narrow. """ - codes, sources = pd.factorize(cells["source"]) - draws = _multinomial_draws(len(sources), replicates, seed) + universe = pd.Index(pd.unique(np.asarray(universe))) + codes = _codes(cells["cluster"], universe) + draws = _cluster_draws(len(universe), replicates, seed) band = assign_bands(cells["nsnd"], edges) tail = (1.0 - level) / 2.0 weights = cells["w"].to_numpy(dtype=float) @@ -161,19 +223,23 @@ def bootstrap( def bootstrap_total( values: object, - sources: object, + clusters: object, *, replicates: int = 1_000, seed: int = 0, level: float = 0.9, ) -> dict[str, float]: - """The same source bootstrap for one weighted total (``values`` already weighted).""" + """The same cluster bootstrap for one weighted total over the whole domain. + + ``values`` are already weighted, one per person, and ``clusters`` cover + every person in the domain. + """ values = np.asarray(values, dtype=float) - codes, unique = pd.factorize(pd.Series(sources)) - draws = _multinomial_draws(len(unique), replicates, seed) - by_source = np.bincount(codes, weights=values, minlength=len(unique)) - totals = draws @ by_source + codes, unique = pd.factorize(pd.Series(clusters)) + draws = _cluster_draws(len(unique), replicates, seed) + by_cluster = np.bincount(codes, weights=values, minlength=len(unique)) + totals = draws @ by_cluster tail = (1.0 - level) / 2.0 low, high = np.quantile(totals, [tail, 1 - tail]) estimate = float(values.sum()) @@ -182,55 +248,68 @@ def bootstrap_total( "low": float(low), "high": float(high), "cv": float(totals.std() / abs(estimate)) if estimate else float("nan"), - "sources": int(len(unique)), - "top_source_share": float(np.abs(by_source).max() / np.abs(by_source).sum()) - if len(unique) - else 0.0, } -def _primary( - table: pd.DataFrame, - amount_columns: Sequence[str], - min_sources: int, - dominance: float, -) -> list[str]: - statuses = [] - for _, row in table.iterrows(): - shares = [row["top_source_weight_share"]] + [ - row[f"top_source_{name}_share"] for name in amount_columns - ] - if row["sources"] < min_sources: - statuses.append(FEWER_THAN_FLOOR) - elif max(shares) > dominance: - statuses.append(DOMINATED) - else: - statuses.append(PUBLISHED) - return statuses +def record_shares(values: object, sources: object) -> dict[str, float]: + """Distinct contributing records and the top one and two records' shares + of a total, for the dominance check on a single published number.""" + + by_source = ( + pd.Series(np.asarray(values, dtype=float)).groupby(np.asarray(sources)).sum() + ) + by_source = by_source[by_source != 0].abs() + top, top_two = _top_shares(by_source) + return { + "records": int(len(by_source)), + "top_record_share": top, + "top_two_record_share": top_two, + } + + +def passes_primary( + records: int, + top_shares: Sequence[float], + top_two_shares: Sequence[float], + *, + min_sources: int = MIN_SOURCES, + dominance: float = DOMINANCE_SHARE, + top_two: float = TOP_TWO_SHARE, +) -> str: + if records < min_sources: + return FEWER_THAN_FLOOR + if ( + max(top_shares, default=0.0) > dominance + or max(top_two_shares, default=0.0) > top_two + ): + return DOMINATED + return PUBLISHED + + +def _primary(table: pd.DataFrame, amount_columns: Sequence[str], **rules) -> list[str]: + columns = ["weight", *amount_columns] + return [ + passes_primary( + int(row["sources"]), + [row[f"top_source_{name}_share"] for name in columns], + [row[f"top_two_source_{name}_share"] for name in columns], + **rules, + ) + for _, row in table.iterrows() + ] def _passes( - cells: pd.DataFrame, - mask: np.ndarray, - amount_columns: Sequence[str], - min_sources: int, - dominance: float, + cells: pd.DataFrame, mask: np.ndarray, amount_columns: Sequence[str], **rules ) -> bool: """Whether the persons under ``mask``, taken as one cell, pass the primary rules.""" - members = cells[mask] - weighted = members[["source"]].assign( - w=members["w"], - **{name: members["w"] * members[name] for name in amount_columns}, + by_source = _by_source(cells[mask], amount_columns) + shares = [_top_shares(by_source[column]) for column in ("w", *amount_columns)] + status = passes_primary( + len(by_source), [s[0] for s in shares], [s[1] for s in shares], **rules ) - by_source = weighted.groupby("source").sum() - if len(by_source) < min_sources: - return False - for column in ("w", *amount_columns): - total = by_source[column].sum() - if total > 0 and by_source[column].max() / total > dominance: - return False - return True + return status == PUBLISHED def disclosure_statuses( @@ -242,33 +321,38 @@ def disclosure_statuses( total_published: bool = True, min_sources: int = MIN_SOURCES, dominance: float = DOMINANCE_SHARE, + top_two: float = TOP_TWO_SHARE, ) -> tuple[list[str], list[str]]: """Statuses for Release 1 (the ``fine`` bands) and Release 2 (the ``coarse`` bins) together, read as one hierarchy: total, bins, bands. - Primary: fewer than ``min_sources`` sources, or one source above - ``dominance`` of the cell's weight or of any amount. + Primary: fewer than ``min_sources`` records, one record above + ``dominance`` of the cell's weight or of any amount, or the top two above + ``top_two``. Secondary: whatever a reader can rebuild by subtraction must itself pass - the primary rules, so cells are withheld, fewest sources first, until: + the primary rules, so cells are withheld, fewest records first, until: - the withheld bands inside each published bin pass, taken together; - no withheld bin can be rebuilt from its bands; - with ``total_published``, the withheld bins pass together, and so do the withheld bands inside them. The 1p costing is about one per cent of the total property slice, so it gives the total away, and it is published. + + The same holds for record counts, which ``render`` shows for published + cells only. """ + rules = {"min_sources": min_sources, "dominance": dominance, "top_two": top_two} fine = audit(cells, fine_edges, amount_columns) coarse = audit(cells, coarse_edges, amount_columns) - fine_status = _primary(fine, amount_columns, min_sources, dominance) - coarse_status = _primary(coarse, amount_columns, min_sources, dominance) + fine_status = _primary(fine, amount_columns, **rules) + coarse_status = _primary(coarse, amount_columns, **rules) band = assign_bands(cells["nsnd"], fine_edges) holder = coarse_index_of_fine(fine_edges, coarse_edges) def passes(bands: Sequence[int]) -> bool: - mask = np.isin(band, list(bands)) - return _passes(cells, mask, amount_columns, min_sources, dominance) + return _passes(cells, np.isin(band, list(bands)), amount_columns, **rules) def fewest(candidates: Sequence[int], table: pd.DataFrame) -> int: return min( @@ -323,13 +407,14 @@ def render( mode: str, count_step: float = 100.0, amount_step: float = 100_000.0, - min_sources: int = MIN_SOURCES, ) -> pd.DataFrame: """A release table. ``restricted`` shows every band in full, with exact record counts and the status as a note; it is for recipients covered by the licence. ``public`` - blanks withheld bands and shows record counts under the floor as ``<10``. + shows a band's values and record count only when it is published. A + withheld band's record count would otherwise give a thin band's count by + subtraction from its bin. """ if mode not in {"restricted", "public"}: @@ -339,12 +424,11 @@ def render( interval = intervals.loc[index] status = statuses[index] shown = mode == "restricted" or status == PUBLISHED - out = {"band": row["band"]} - if mode == "public" and row["sources"] < min_sources: - out["distinct_records"] = f"<{min_sources}" - else: - out["distinct_records"] = str(int(row["sources"])) - out["effective_records"] = round(float(row["source_ess"]), 1) if shown else None + out = { + "band": row["band"], + "distinct_records": str(int(row["sources"])) if shown else WITHHELD_RECORDS, + "effective_records": round(float(row["source_ess"]), 1) if shown else None, + } for name, step in [ ("count", count_step), *[(n, amount_step) for n in amount_columns], diff --git a/tools/uk_scotland_property_income/spi_2023_24_published.csv b/tools/uk_scotland_property_income/spi_2023_24_published.csv new file mode 100644 index 000000000..023d4a799 --- /dev/null +++ b/tools/uk_scotland_property_income/spi_2023_24_published.csv @@ -0,0 +1,57 @@ +fact,table,published,unit +hmrc/property_income_count_income_band_12_570_to_15_000,3.7,101,1000 +hmrc/property_income_income_band_12_570_to_15_000,3.7,670,1000000 +hmrc/property_income_count_income_band_15_000_to_20_000,3.7,232,1000 +hmrc/property_income_income_band_15_000_to_20_000,3.7,1690,1000000 +hmrc/property_income_count_income_band_20_000_to_30_000,3.7,431,1000 +hmrc/property_income_income_band_20_000_to_30_000,3.7,3820,1000000 +hmrc/property_income_count_income_band_30_000_to_40_000,3.7,364,1000 +hmrc/property_income_income_band_30_000_to_40_000,3.7,3800,1000000 +hmrc/property_income_count_income_band_40_000_to_50_000,3.7,325,1000 +hmrc/property_income_income_band_40_000_to_50_000,3.7,3780,1000000 +hmrc/property_income_count_income_band_50_000_to_70_000,3.7,397,1000 +hmrc/property_income_income_band_50_000_to_70_000,3.7,5400,1000000 +hmrc/property_income_count_income_band_70_000_to_100_000,3.7,226,1000 +hmrc/property_income_income_band_70_000_to_100_000,3.7,3770,1000000 +hmrc/property_income_count_income_band_100_000_to_150_000,3.7,137,1000 +hmrc/property_income_income_band_100_000_to_150_000,3.7,2550,1000000 +hmrc/property_income_count_income_band_150_000_to_200_000,3.7,54,1000 +hmrc/property_income_income_band_150_000_to_200_000,3.7,1150,1000000 +hmrc/property_income_count_income_band_200_000_to_300_000,3.7,41,1000 +hmrc/property_income_income_band_200_000_to_300_000,3.7,1050,1000000 +hmrc/property_income_count_income_band_300_000_to_500_000,3.7,23,1000 +hmrc/property_income_income_band_300_000_to_500_000,3.7,765,1000000 +hmrc/property_income_count_income_band_500_000_to_1_000_000,3.7,13,1000 +hmrc/property_income_income_band_500_000_to_1_000_000,3.7,474,1000000 +hmrc/property_income_count_income_band_1_000_000_to_inf,3.7,7,1000 +hmrc/property_income_income_band_1_000_000_to_inf,3.7,439,1000000 +hmrc.spi_region.taxpayers_by_region_12570_15000@S92000003,3.11,238,1000 +hmrc.spi_region.total_income_by_region_12570_15000@S92000003,3.11,3290,1000000 +hmrc.spi_region.income_tax_by_region_12570_15000@S92000003,3.11,47,1000000 +hmrc.spi_region.taxpayers_by_region_15000_20000@S92000003,3.11,450,1000 +hmrc.spi_region.total_income_by_region_15000_20000@S92000003,3.11,7860,1000000 +hmrc.spi_region.income_tax_by_region_15000_20000@S92000003,3.11,374,1000000 +hmrc.spi_region.taxpayers_by_region_20000_30000@S92000003,3.11,841,1000 +hmrc.spi_region.total_income_by_region_20000_30000@S92000003,3.11,20900,1000000 +hmrc.spi_region.income_tax_by_region_20000_30000@S92000003,3.11,1880,1000000 +hmrc.spi_region.taxpayers_by_region_30000_40000@S92000003,3.11,573,1000 +hmrc.spi_region.total_income_by_region_30000_40000@S92000003,3.11,19800,1000000 +hmrc.spi_region.income_tax_by_region_30000_40000@S92000003,3.11,2360,1000000 +hmrc.spi_region.taxpayers_by_region_40000_50000@S92000003,3.11,364,1000 +hmrc.spi_region.total_income_by_region_40000_50000@S92000003,3.11,16300,1000000 +hmrc.spi_region.income_tax_by_region_40000_50000@S92000003,3.11,2240,1000000 +hmrc.spi_region.taxpayers_by_region_50000_70000@S92000003,3.11,298,1000 +hmrc.spi_region.total_income_by_region_50000_70000@S92000003,3.11,17200,1000000 +hmrc.spi_region.income_tax_by_region_50000_70000@S92000003,3.11,3210,1000000 +hmrc.spi_region.taxpayers_by_region_70000_100000@S92000003,3.11,132,1000 +hmrc.spi_region.total_income_by_region_70000_100000@S92000003,3.11,10800,1000000 +hmrc.spi_region.income_tax_by_region_70000_100000@S92000003,3.11,2690,1000000 +hmrc.spi_region.taxpayers_by_region_100000_150000@S92000003,3.11,56,1000 +hmrc.spi_region.total_income_by_region_100000_150000@S92000003,3.11,6560,1000000 +hmrc.spi_region.income_tax_by_region_100000_150000@S92000003,3.11,2010,1000000 +hmrc.spi_region.taxpayers_by_region_150000_200000@S92000003,3.11,15,1000 +hmrc.spi_region.total_income_by_region_150000_200000@S92000003,3.11,2650,1000000 +hmrc.spi_region.income_tax_by_region_150000_200000@S92000003,3.11,916,1000000 +hmrc.spi_region.taxpayers_by_region_200000_plus@S92000003,3.11,16,1000 +hmrc.spi_region.total_income_by_region_200000_plus@S92000003,3.11,7120,1000000 +hmrc.spi_region.income_tax_by_region_200000_plus@S92000003,3.11,2800,1000000 From a2d6939e5c0e33f4f69e36bb9ca9c2db5c4c1106 Mon Sep 17 00:00:00 2001 From: Maria Juaristi Date: Wed, 21 Oct 2026 12:27:17 +0100 Subject: [PATCH 3/4] Say in the release note how far the reweight moves Scottish landlords Co-Authored-By: Claude Opus 5.5 --- tools/uk_scotland_property_income/__main__.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tools/uk_scotland_property_income/__main__.py b/tools/uk_scotland_property_income/__main__.py index f85b7dfed..f88cbee5d 100644 --- a/tools/uk_scotland_property_income/__main__.py +++ b/tools/uk_scotland_property_income/__main__.py @@ -591,9 +591,11 @@ def release_note(manifest: dict[str, Any], mode: str) -> str: HMRC published. The targets are UK property income by range of total income (counts and amounts, Table 3.7) and Scotland's taxpayers and total income by range (Table 3.11). Each region's household total is held to within - {region_miss:.2%}. Because those figures are calibrated, they don't test the - data. The split of Scottish property income across NSND bands is not - targeted, so that split is what the outturn tests. + {max(region_miss, 1e-4):.2%}. Because those figures are calibrated, they + don't test the data. The split of Scottish property income across NSND + bands is not targeted, so that split is what the outturn tests. It is not + untouched, though: the UK property targets move Scottish landlords' weights + along with everyone else's. - **Law:** 2023-24 income tax. ## Uncertainty @@ -614,7 +616,8 @@ def release_note(manifest: dict[str, Any], mode: str) -> str: band limits as the other NSND rates. - **Basis:** static, on 2023-24 data and law, with no behavioural response. - **Headline and sensitivity:** the headline keeps the finance-cost reduction - at the basic rate. The sensitivity raises it by 1p as well. + at 20%, as policyengine-uk reads current law for Scottish taxpayers. The + sensitivity raises it by 1p as well. - **Later years:** a 2027-28 figure, when the power applies, needs a separate run. From 1552cdc8b05601eb7852ee8cb53ccfa2690cc588 Mon Sep 17 00:00:00 2001 From: Maria Juaristi Date: Wed, 21 Oct 2026 15:26:45 +0100 Subject: [PATCH 4/4] Reweight to Scottish taxpayer counts, not Scotland's total income On 6cba5d71, meeting Scotland's total income by band (Table 3.11) moved weight onto a few very high-income records. Scotland's untargeted investment, employment and pension income all moved away from the published table, and with the UK property facts it took 13% off Scottish taxable property income. Reweighting to UK property income by band and Scottish taxpayer counts by band keeps every untargeted Scottish check at or better than the back-cast. With fewer facts the misses reach zero before the divergence stops shrinking, so the search can end at the iteration limit with the facts met. The limit rises to 20,000; the run accepts a stop there when every miss is within --max-miss, and refuses any other unconverged stop. Co-Authored-By: Claude Opus 5.5 --- .../uk/test_uk_scotland_property_income.py | 8 +++++ tools/uk_scotland_property_income/__main__.py | 30 +++++++++++++------ tools/uk_scotland_property_income/rake.py | 14 +++++++-- 3 files changed, 41 insertions(+), 11 deletions(-) diff --git a/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py b/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py index 037d9565d..9e836cc64 100644 --- a/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py +++ b/packages/microcosm-build/tests/engine_free/uk/test_uk_scotland_property_income.py @@ -185,6 +185,14 @@ def test_rake_stops_at_its_bounds_when_a_fact_is_out_of_reach(): assert result.after[0] == pytest.approx(20.0) +def test_rake_reports_the_iteration_limit_apart_from_failure(): + matrix = sparse.csr_matrix(np.ones((10, 1))) + result = rake.rake(matrix, np.ones(10), [12.0], max_iterations=1) + assert not result.converged + assert result.stopped_at_iteration_limit + assert rake.rake(matrix, np.ones(10), [12.0]).converged + + def _persons_table(): table = pd.DataFrame( { diff --git a/tools/uk_scotland_property_income/__main__.py b/tools/uk_scotland_property_income/__main__.py index f88cbee5d..24c2952c7 100644 --- a/tools/uk_scotland_property_income/__main__.py +++ b/tools/uk_scotland_property_income/__main__.py @@ -199,10 +199,14 @@ def _reweight( facts: list[registry.SpiFact], args: argparse.Namespace, ) -> tuple[rake.RakeResult, list[registry.SpiFact], list[str], dict[str, float]]: + # UK property income by band, and Scotland's taxpayers by band. Scotland's + # total income by band is left out: on 6cba5d71 meeting it moved weight onto + # a few very high-income records, so Scotland's untargeted investment, + # employment and pension income all moved away from Table 3.11. raked_facts = [fact for fact in facts if fact.region is None] + [ fact for fact in facts - if fact.region == "SCOTLAND" and fact.measure in {"taxpayers", "total_income"} + if fact.region == "SCOTLAND" and fact.measure == "taxpayers" ] w0 = earlier.household["household_weight"].to_numpy(dtype=float) household_regions = release.as_text(earlier.household["region"]).to_numpy() @@ -228,9 +232,10 @@ def _reweight( "facts": float(errors[: len(raked_facts)].max()), "region_households": float(errors[len(raked_facts) :].max()), } - if not raked.converged or max(misses.values()) > args.max_miss: + settled = raked.converged or raked.stopped_at_iteration_limit + if not settled or max(misses.values()) > args.max_miss: raise SystemExit( - f"The reweight didn't settle: converged={raked.converged}, " + f"The reweight didn't settle: {raked.message}, " f"largest misses {misses}, allowed {args.max_miss}." ) return raked, raked_facts, regions, misses @@ -436,6 +441,7 @@ def tax_change(delta: float, relief_follows: bool) -> np.ndarray: "region_controls": regions, "ratio_bounds": list(args.ratio_bounds), "converged": raked.converged, + "stopped_at_iteration_limit": raked.stopped_at_iteration_limit, "message": raked.message, "max_abs_error_before": float(np.abs(raked.relative_error_before).max()), "max_abs_error_after": misses, @@ -589,13 +595,19 @@ def release_note(manifest: dict[str, Any], mode: str) -> str: forward, stay at their 2024-25 levels. - **Reweight:** the back-cast weights are reweighted to the SPI 2023-24 tables HMRC published. The targets are UK property income by range of total income - (counts and amounts, Table 3.7) and Scotland's taxpayers and total income by + (counts and amounts, Table 3.7) and the number of Scottish taxpayers by range (Table 3.11). Each region's household total is held to within - {max(region_miss, 1e-4):.2%}. Because those figures are calibrated, they - don't test the data. The split of Scottish property income across NSND - bands is not targeted, so that split is what the outturn tests. It is not - untouched, though: the UK property targets move Scottish landlords' weights - along with everyone else's. + {max(region_miss, 1e-4):.2%}. + - **Not a test:** because those figures are calibrated, they don't test the + data. + - **Left out on purpose:** Scotland's total income by range. Meeting it moved + weight onto a few very high-income records and pulled Scotland's other + income away from the published table. + - **What the outturn tests:** the split of Scottish property income across + NSND bands, which is not targeted. It isn't untouched, though. The UK + property targets lower Scottish landlords' weights along with everyone + else's, while matching the taxpayer counts moves Scottish taxpayers up the + bands, where landlords are more common. - **Law:** 2023-24 income tax. ## Uncertainty diff --git a/tools/uk_scotland_property_income/rake.py b/tools/uk_scotland_property_income/rake.py index 2ce904aff..3c7c5ddd7 100644 --- a/tools/uk_scotland_property_income/rake.py +++ b/tools/uk_scotland_property_income/rake.py @@ -27,6 +27,13 @@ class RakeResult: after: np.ndarray converged: bool message: str + #: L-BFGS-B's status: 0 converged, 1 stopped at the iteration limit, 2 stopped + #: abnormally (a failed line search, for instance). + status: int + + @property + def stopped_at_iteration_limit(self) -> bool: + return self.status == 1 @property def relative_error_before(self) -> np.ndarray: @@ -54,11 +61,13 @@ def rake( lower: float = 0.25, upper: float = 4.0, divergence_weight: float = 1e-3, - max_iterations: int = 5_000, + max_iterations: int = 20_000, ) -> RakeResult: """Ratios on ``initial_weights`` so that ``matrix.T @ weights`` meets ``targets``. - ``matrix`` is households by facts. + ``matrix`` is households by facts. With few facts the misses reach zero + long before the divergence stops shrinking, so the search can end at + ``max_iterations`` with the facts met; the caller judges the misses. """ a = sparse.csr_matrix(matrix, dtype=float) @@ -107,4 +116,5 @@ def objective(log_ratio: np.ndarray) -> tuple[float, np.ndarray]: after=a_t @ weights, converged=bool(solution.success), message=str(solution.message), + status=int(solution.status), )