diff --git a/README.md b/README.md index b7da84050..754f1c68e 100644 --- a/README.md +++ b/README.md @@ -87,6 +87,7 @@ before running UK population examples. To download the raw `.h5` file directly, - `examples/income_distribution_us.py`: Analyse benefit distribution by decile - `examples/employment_income_variation_uk.py`: Model employment income phase-outs - `examples/policy_change_uk.py`: Analyse policy reform impacts +- `examples/household_targeting_uk.py`: Replicate a published household-level targeting rule (UK) - `examples/paper_repro_uk.py`: Reproduce the UK reform analysis used in the JOSS paper draft ## Installation diff --git a/changelog.d/uk-household-targeting-example.added.md b/changelog.d/uk-household-targeting-example.added.md new file mode 100644 index 000000000..897b63be1 --- /dev/null +++ b/changelog.d/uk-household-targeting-example.added.md @@ -0,0 +1 @@ +Added `examples/household_targeting_uk.py`, a worked replication of a published household targeting rule on the UK microdata. diff --git a/docs/examples.md b/docs/examples.md index fba2cb0ca..6f63b7bcb 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -23,6 +23,13 @@ Estimating results by area from per-household arrays, checking the release, with ```{.python include="../examples/uk_geography_disclosure.py"} ``` +## UK household targeting rule + +Replicating a published analysis whose eligibility rule is not a PolicyEngine variable: the highest individual income in a household, passporting by reported or modelled benefit receipt, after-housing-costs deciles, benefit units joined to households, and the cost of a fixed budget. Needs access to the private UK data. + +```{.python include="../examples/household_targeting_uk.py"} +``` + ## UK income bands Calculating net income and tax by income decile using representative microdata and `Aggregate` with quantile filters. diff --git a/examples/household_targeting_uk.py b/examples/household_targeting_uk.py new file mode 100644 index 000000000..fa6806f0f --- /dev/null +++ b/examples/household_targeting_uk.py @@ -0,0 +1,202 @@ +"""Example: replicate a published household targeting analysis on UK microdata. + +Many published analyses define eligibility with a rule that is not a +PolicyEngine variable. This script works through one: the energy bill discount +options in the Resolution Foundation's "Billing me softly" (August 2026, +https://www.resolutionfoundation.org/publications/billing-me-softly/), which +tests the highest individual taxable income in a household against GBP 24,000, +and passports households that report a means-tested benefit. + +It shows the steps such a replication needs: + +1. Load the certified UK microdata with ``managed_microsimulation``. +2. Read person, benefit unit and household arrays with their weights. +3. Aggregate a person variable to the household (here, the maximum). +4. Join benefit units to households through the person-level id columns. +5. Choose between benefit receipt reported in the survey and modelled + entitlement (``*_reported`` versus the modelled benefit variables). +6. Rank households into income deciles on the basis the publication uses. +7. Turn eligibility into counts, shares and per-household amounts. + +The UK data is private. Set HUGGING_FACE_TOKEN to a token from a Hugging Face +account with access (see docs/microsim.md), so this script does not run in CI. + +A full implementation of these options as a policyengine-uk reform, with the +published figures beside PolicyEngine's, is in +https://github.com/PolicyEngine/uk-energy-reforms. + +Run: python examples/household_targeting_uk.py +""" + +import numpy as np +import pandas as pd + +import policyengine as pe + +YEAR = 2024 # fiscal year 2024-25, the survey year of the publication +INDIVIDUAL_THRESHOLD = 24_000 +LOWER_TIER_THRESHOLD = 18_000 +HOUSEHOLD_THRESHOLD = 30_000 +BUDGET = 2e9 + +# The Warm Home Discount passporting benefits. +PASSPORTING_BENEFITS = [ + "universal_credit", + "pension_credit", + "housing_benefit", + "esa_income", + "jsa_income", + "income_support", +] + + +def values(sim, variable: str) -> np.ndarray: + """Return a variable's values as a plain array on its own entity. + + ``sim.calculate`` returns a weighted MicroSeries. Take ``.values`` only + when you apply the weights yourself, as this script does. + """ + return np.asarray(sim.calculate(variable, YEAR).values) + + +def household_max(person_values, person_household_id, household_id): + """Maximum of a person-level value within each household.""" + by_household = pd.Series(person_values).groupby(person_household_id).max() + return by_household.reindex(household_id).to_numpy() + + +def household_any(person_values, person_household_id, household_id): + """Whether any person in each household has a positive value.""" + by_household = pd.Series(person_values > 0).groupby(person_household_id).any() + return by_household.reindex(household_id).to_numpy(dtype=bool) + + +def weighted_decile(income, weight): + """Decile (1 to 10) of each row, with boundaries set by ``weight``.""" + order = np.argsort(income, kind="stable") + cumulative = np.cumsum(weight[order]) + decile = np.empty(len(income)) + decile[order] = np.ceil(10 * cumulative / cumulative[-1]) + return np.clip(decile, 1, 10).astype(int) + + +def main(): + # The certified default dataset, pinned by the policyengine.py release. + sim = pe.uk.managed_microsimulation() + + # Households. The publication covers Great Britain. + household_id = values(sim, "household_id") + household_weight = values(sim, "household_weight") + in_gb = values(sim, "country").astype(str) != "NORTHERN_IRELAND" + weight = household_weight * in_gb + + # Persons. person_household_id links each person to a household. + person_household_id = values(sim, "person_household_id") + person_benunit_id = values(sim, "person_benunit_id") + taxable_income = values(sim, "total_income") + + # 3. Highest individual income in the household. + highest_income = household_max(taxable_income, person_household_id, household_id) + individual_test = highest_income < INDIVIDUAL_THRESHOLD + + # The household option: total taxable income, equivalised. policyengine-uk's + # household_equivalisation_bhc is the modified OECD scale rescaled so that a + # couple is 1.0, as in HBAI. Check which form a publication uses: on the + # unrescaled form (single adult 1.0) the same threshold passes far more + # households. + household_income = ( + pd.Series(taxable_income) + .groupby(person_household_id) + .sum() + .reindex(household_id) + .to_numpy() + ) + equivalised = household_income / values(sim, "household_equivalisation_bhc") + household_test = equivalised < HOUSEHOLD_THRESHOLD + + # 5. Passporting. The *_reported variables hold receipt reported in the + # survey; the modelled variables hold entitlement after take-up. Published + # analyses of survey data usually mean reported receipt. + reported = np.zeros(len(person_household_id)) + for benefit in PASSPORTING_BENEFITS: + reported += values(sim, f"{benefit}_reported") + passported_reported = household_any(reported, person_household_id, household_id) + + modelled = sum( + np.asarray(sim.calculate(benefit, YEAR, map_to="household").values) + for benefit in PASSPORTING_BENEFITS + ) + passported_modelled = modelled > 0 + + # 6. Deciles. household_income_decile ranks people by equivalised income + # before housing costs. HBAI headline distributions, and many publications, + # use income after housing costs, so build that ranking explicitly and say + # whether households or people set the boundaries. + income_ahc = values(sim, "equiv_hbai_household_net_income_ahc") + people = values(sim, "household_count_people") + decile_ahc = weighted_decile(income_ahc, weight * people) + lowest_four = decile_ahc <= 4 + + def share(mask, among=None): + among = np.ones(len(mask), dtype=bool) if among is None else among + return (weight * mask * among).sum() / (weight * among).sum() + + print(f"GB households: {weight.sum() / 1e6:.1f}m") + print(f"Highest individual income below £24,000: {share(individual_test):.0%}") + print( + " share of the four lowest AHC deciles: " + f"{share(individual_test, lowest_four):.0%}" + ) + print(f"Equivalised household income below £30,000: {share(household_test):.0%}") + print(f"Passported, reported receipt: {share(passported_reported):.0%}") + print(f"Passported, modelled entitlement: {share(passported_modelled):.0%}") + + # 7. A fixed budget shared among eligible households. + eligible = individual_test | passported_reported + recipients = (weight * eligible).sum() + print( + f"£{BUDGET / 1e9:.0f}bn among {recipients / 1e6:.1f}m eligible households: " + f"£{BUDGET / recipients:.0f} each" + ) + upper = (highest_income < LOWER_TIER_THRESHOLD) | passported_reported + lower = individual_test & ~upper + cost = (weight * (220 * upper + 85 * lower)).sum() + print(f"Two tiers of £220 and £85 cost £{cost / 1e9:.2f}bn") + + # 4. Benefit units to households. Join through the person-level ids. + # Do not use sim.calculate("household_id", map_to="benunit"): mapping from a + # larger entity to a smaller one shares the value among the smaller units, + # which is right for amounts and wrong for ids and flags. + benunit_id = values(sim, "benunit_id") + benunit_household = ( + pd.Series(person_household_id) + .groupby(person_benunit_id) + .first() + .reindex(benunit_id) + .to_numpy() + ) + + def by_benunit(household_values): + return ( + pd.Series(household_values, index=household_id) + .reindex(benunit_household) + .to_numpy() + ) + + benunit_weight = values(sim, "benunit_weight") * by_benunit(in_gb) + couple_with_children = values(sim, "is_couple").astype(bool) & ( + values(sim, "benunit_count_children") > 0 + ) + loses = ( + couple_with_children + & by_benunit(household_test).astype(bool) + & ~by_benunit(individual_test).astype(bool) + ) + print( + "Couples with children passing the household test but not the " + f"individual one: {(benunit_weight * loses).sum() / 1e3:,.0f}k" + ) + + +if __name__ == "__main__": + main()