Anchor bus fare distribution to DfT BUS05ai and NTS0705a

changelog.d/bus-fare-distribution-anchoring.md

@@ -0,0 +1 @@
+- Anchor the distribution of bus_fare_spending to DfT BUS05ai (London/outside-London receipts split) and NTS0705a (local-bus trip rates by income quintile), preserving the national BUS05aii total, with post-build validation tests.

policyengine_uk_data/datasets/imputations/consumption.py

@@ -797,12 +797,53 @@ def calibrate_dataset_fuel_litre_proxies_to_road_fuel(
     2025: 3.4e9 * ENGLAND_TO_UK_POPULATION_UPLIFT,
 }
 
+# DfT BUS05ai, year ending March 2025: England passenger fare receipts by area
+# — England GBP 3,417m, of which London GBP 1,347m and England outside London
+# GBP 2,070m. The LCFS imputation under-captures London fares (Oyster/
+# contactless spend is poorly separated in household diaries), and a
+# national-only calibration smears the shortfall across the other English
+# regions, overstating outside-London fare exposure by ~47%.
+# https://www.gov.uk/government/statistical-data-sets/bus-statistics-data-tables
+BUS_FARE_LONDON_SHARE_OF_ENGLAND = 1347 / 3417
+
+# NTS0705a 2024, 'Other local bus' (outside London) average trips per person
+# per year by household income quintile (England). Anchors the income
+# distribution of household fare spending, assuming a common average fare per
+# trip across quintiles: the raw LCFS spend gradient runs the wrong way
+# (richer households record more bus & coach spend, while NTS shows the lowest
+# quintile makes ~3.7x the local bus trips of the highest).
+# https://www.gov.uk/government/statistical-data-sets/nts07-car-ownership-and-access
+NTS_BUS_TRIPS_BY_INCOME_QUINTILE = {1: 47.6, 2: 37.5, 3: 24.8, 4: 18.5, 5: 13.0}
+
+ENGLAND_REGIONS_OUTSIDE_LONDON = {
+    "NORTH_EAST",
+    "NORTH_WEST",
+    "YORKSHIRE",
+    "EAST_MIDLANDS",
+    "WEST_MIDLANDS",
+    "EAST_OF_ENGLAND",
+    "SOUTH_EAST",
+    "SOUTH_WEST",
+}
+
 
 def calibrate_bus_fare_spending(
     dataset: UKSingleYearDataset,
     time_period: int,
 ) -> float | None:
-    """Scale bus_fare_spending to the DfT passenger-fare total (BUS_FARE_TARGETS)."""
+    """Calibrate bus_fare_spending to DfT and NTS anchors.
+
+    Three constraints, applied in order (the distributional anchors preserve
+    the national total):
+
+    1. National total: DfT BUS05aii passenger fare receipts (BUS_FARE_TARGETS).
+    2. Within-England regional split: DfT BUS05ai London / outside-London
+       receipts (BUS_FARE_LONDON_SHARE_OF_ENGLAND).
+    3. Income distribution: England quintile fare shares proportional to
+       population-weighted NTS0705a local-bus trip rates
+       (NTS_BUS_TRIPS_BY_INCOME_QUINTILE), alternated with (2) in a small
+       IPF loop.
+    """
     target = BUS_FARE_TARGETS.get(time_period)
     if target is None:
         return None
@@ -811,18 +852,86 @@ def calibrate_bus_fare_spending(
     dataset.time_period = str(original_time_period)
     try:
         simulation = Microsimulation(dataset=dataset)
-        actual = simulation.calculate(
+        fares = simulation.calculate(
             "bus_fare_spending",
             period=time_period,
             map_to="household",
-        ).sum()
+        )
+        actual = fares.sum()
+        region = np.array(
+            simulation.calculate(
+                "region", period=time_period, map_to="household"
+            ).values
+        ).astype(str)
+        decile = np.array(
+            simulation.calculate(
+                "household_income_decile", period=time_period, map_to="household"
+            ).values,
+            dtype=float,
+        )
+        people = np.array(
+            simulation.calculate(
+                "people", period=time_period, map_to="household"
+            ).values,
+            dtype=float,
+        )
+        weights = np.array(fares.weights, dtype=float)
     finally:
         dataset.time_period = original_time_period
     if actual <= 0:
         raise ValueError(f"Cannot calibrate bus_fare_spending: aggregate is {actual}.")
 
     scale = target / actual
-    dataset.household["bus_fare_spending"] *= scale
+    values = np.array(dataset.household["bus_fare_spending"], dtype=float) * scale
+
+    is_london = region == "LONDON"
+    outside_london = np.isin(region, list(ENGLAND_REGIONS_OUTSIDE_LONDON))
+    in_england = is_london | outside_london
+    quintile = np.where(
+        decile >= 1, np.clip(((decile - 1) // 2 + 1).astype(int), 1, 5), 0
+    )
+
+    def _wsum(mask_or_values):
+        return float((np.asarray(mask_or_values, dtype=float) * weights).sum())
+
+    england_total = _wsum(values * in_england)
+    quintile_population = {
+        q: _wsum(people * (in_england & (quintile == q)))
+        for q in NTS_BUS_TRIPS_BY_INCOME_QUINTILE
+    }
+    target_norm = sum(
+        quintile_population[q] * rate
+        for q, rate in NTS_BUS_TRIPS_BY_INCOME_QUINTILE.items()
+    )
+    for _ in range(3):
+        for q, rate in NTS_BUS_TRIPS_BY_INCOME_QUINTILE.items():
+            mask = in_england & (quintile == q)
+            current = _wsum(values * mask)
+            if current > 0:
+                quintile_target = (
+                    quintile_population[q] * rate / target_norm * england_total
+                )
+                values = np.where(mask, values * quintile_target / current, values)
+        london_actual = _wsum(values * is_london)
+        outside_actual = _wsum(values * outside_london)
+        if london_actual > 0 and outside_actual > 0:
+            values = np.where(
+                is_london,
+                values
+                * BUS_FARE_LONDON_SHARE_OF_ENGLAND
+                * england_total
+                / london_actual,
+                np.where(
+                    outside_london,
+                    values
+                    * (1 - BUS_FARE_LONDON_SHARE_OF_ENGLAND)
+                    * england_total
+                    / outside_actual,
+                    values,
+                ),
+            )
+
+    dataset.household["bus_fare_spending"] = values
     return scale
 
 

policyengine_uk_data/tests/test_bus_fare_distribution.py

@@ -0,0 +1,103 @@
+"""Post-build validation of the bus fare spending distribution.
+
+These tests run against the *built* Enhanced FRS (via the ``baseline``
+fixture, which skips when the dataset has not been rebuilt), so they
+validate the dataset after a rebuild rather than the calibration code in
+isolation. They lock in the anchors applied by
+``calibrate_bus_fare_spending``:
+
+1. National total: DfT BUS05aii passenger fare receipts, England → UK.
+2. Within-England regional split: DfT BUS05ai London / outside-London
+   receipts (London GBP 1,347m of England's GBP 3,417m). Without this,
+   the LCFS imputation under-captures London and overstates
+   outside-London fare exposure by ~47%.
+3. Income gradient: England quintile fare shares anchored to NTS0705a
+   local-bus trip rates (Q1 47.6 ... Q5 13.0 trips/person/year) — the
+   lowest quintile makes ~3.7x the trips of the highest, so fare
+   spending must decline with income, not rise with it as the raw LCFS
+   spend gradient does.
+"""
+
+from __future__ import annotations
+
+import os
+
+import numpy as np
+
+from policyengine_uk_data.datasets.frs_release import CURRENT_FRS_RELEASE
+from policyengine_uk_data.datasets.imputations.consumption import (
+    BUS_FARE_LONDON_SHARE_OF_ENGLAND,
+    BUS_FARE_TARGETS,
+    ENGLAND_REGIONS_OUTSIDE_LONDON,
+    NTS_BUS_TRIPS_BY_INCOME_QUINTILE,
+)
+
+PERIOD = CURRENT_FRS_RELEASE.calibration_year
+# The anchoring is a deterministic post-processing rescale, but PR CI builds
+# with TESTING=1 (reduced calibration epochs), which shifts weights and decile
+# boundaries slightly relative to a full build — widen tolerances there, per
+# repo convention (cf. test_cgt_band_donors, test_salary_sacrifice_headcount).
+_REDUCED_BUILD = os.environ.get("TESTING") == "1"
+TOTAL_TOLERANCE = 0.10 if _REDUCED_BUILD else 0.05
+SHARE_TOLERANCE = 0.05 if _REDUCED_BUILD else 0.03  # abs points, London/England
+
+
+def _household_arrays(baseline):
+    fares = baseline.calculate("bus_fare_spending", PERIOD, map_to="household")
+    values = np.array(fares.values) * np.array(fares.weights)
+    region = np.array(
+        baseline.calculate("region", PERIOD, map_to="household").values
+    ).astype(str)
+    decile = np.array(
+        baseline.calculate(
+            "household_income_decile", PERIOD, map_to="household"
+        ).values,
+        dtype=float,
+    )
+    return values, region, decile
+
+
+def test_bus_fare_total_matches_dft_target(baseline):
+    """Weighted bus fare spending matches the DfT BUS05aii UK target."""
+    target = BUS_FARE_TARGETS[PERIOD]
+    total = baseline.calculate("bus_fare_spending", PERIOD).sum()
+    assert abs(total / target - 1) < TOTAL_TOLERANCE, (
+        f"bus_fare_spending total {total / 1e9:.2f}bn is >{TOTAL_TOLERANCE:.0%} "
+        f"from the DfT target {target / 1e9:.2f}bn."
+    )
+
+
+def test_london_fare_share_matches_bus05ai(baseline):
+    """London's share of England fares matches DfT BUS05ai receipts."""
+    values, region, _ = _household_arrays(baseline)
+    london = values[region == "LONDON"].sum()
+    outside = values[np.isin(region, list(ENGLAND_REGIONS_OUTSIDE_LONDON))].sum()
+    share = london / (london + outside)
+    assert abs(share - BUS_FARE_LONDON_SHARE_OF_ENGLAND) < SHARE_TOLERANCE, (
+        f"London fare share {share:.3f} is >{SHARE_TOLERANCE} from the "
+        f"BUS05ai share {BUS_FARE_LONDON_SHARE_OF_ENGLAND:.3f}."
+    )
+
+
+def test_fare_spending_declines_with_income(baseline):
+    """England quintile fare totals follow the NTS trip gradient direction."""
+    values, region, decile = _household_arrays(baseline)
+    in_england = np.isin(region, list(ENGLAND_REGIONS_OUTSIDE_LONDON) + ["LONDON"])
+    quintile = np.where(
+        decile >= 1, np.clip(((decile - 1) // 2 + 1).astype(int), 1, 5), 0
+    )
+    totals = {
+        q: values[in_england & (quintile == q)].sum()
+        for q in NTS_BUS_TRIPS_BY_INCOME_QUINTILE
+    }
+    assert totals[1] > totals[5], (
+        f"Q1 fare spending ({totals[1] / 1e9:.2f}bn) should exceed Q5 "
+        f"({totals[5] / 1e9:.2f}bn) per the NTS0705a trip gradient."
+    )
+    # The full anchored gradient: per-person trip rates decline monotonically,
+    # so quintile totals must not increase from Q3 upwards.
+    assert totals[3] >= totals[4] >= totals[5], (
+        "Upper-quintile fare spending should decline: "
+        f"Q3 {totals[3] / 1e9:.2f}bn, Q4 {totals[4] / 1e9:.2f}bn, "
+        f"Q5 {totals[5] / 1e9:.2f}bn."
+    )