Fit a Pareto upper tail to the capital gains imputation

changelog.d/443.md

@@ -0,0 +1 @@
+- Fix the capital gains imputation discarding the upper tail of the gains distribution. The Advani & Summers percentile table stops at p95, so sampling uniform quantiles through the linear spline linearly extrapolated everything above the 95th percentile off the last segment, leaving no Pareto tail: only ~1.6% of imputed gains came from gains of £1m or more (HMRC: ~61%) and no individual gain exceeded ~£2m. Gains above p95 are now drawn from a Pareto tail anchored at each income band's p95, with the shape parameter fitted per band so the band reproduces the `mean_gains_given_gains` already published in the same source table (which the spline alone undershot by 35-50% in every band). Fitting per income band preserves the existing income-band conditioning, and the seeded generator and the spline body below p95 are unchanged. Adds HMRC's size-of-gain distribution (Capital Gains Tax statistics Table 2.1a) as a committed validation reference.

policyengine_uk_data/datasets/imputations/capital_gains.py

@@ -19,6 +19,128 @@
     -1
 ).fillna(np.inf)
 
+# --- Upper tail of the capital gains distribution -------------------------
+#
+# The Advani & Summers percentile table stops at p95, so fitting a linear
+# spline through p05..p95 and sampling uniform quantiles linearly extrapolates
+# everything above the 95th percentile off the last segment. That removes the
+# Pareto tail of capital gains entirely: gains end up spread thinly over many
+# people instead of concentrated on a few, which both misstates the shape of
+# the distribution and suppresses CGT revenue (small gains sit in lower rate
+# bands and near the annual exempt amount).
+#
+# We keep the spline for the body (below p95) -- that is the Advani & Summers
+# evidence and is not in question -- and replace the extrapolated region with a
+# Pareto tail anchored at p95.
+#
+# The tail index is *fitted*, not assumed. The same source table publishes
+# `mean_gains_given_gains` for each income band, which the spline alone
+# undershoots by 35-50% in every band. That gap is precisely the missing tail
+# mass, so for each band we solve for the Pareto shape parameter that makes the
+# full (body + tail) distribution reproduce the published band mean. Because
+# the fit is done per income band, the income-band conditioning of the original
+# imputation is preserved rather than a single global tail being imposed.
+#
+# The resulting shape parameters cluster around 1.45-1.5, consistent with the
+# heavy tail HMRC's own size-of-gain distribution implies (see
+# storage/capital_gains_size_distribution_hmrc.csv, which is used as the
+# external validation reference for this fit).
+
+TAIL_START_QUANTILE = 0.95
+
+# Shape parameters below ~1.15 give a tail so heavy that the sample mean is
+# dominated by a handful of draws and the build becomes unstable. A small
+# number of bands have an anomalously low p95 relative to their published mean
+# and would otherwise solve to alpha ~ 1.0; we floor them.
+MIN_TAIL_ALPHA = 1.15
+MAX_TAIL_ALPHA = 50.0
+
+# Largest single gain any individual may be assigned. HMRC's top published band
+# is "£5m and over"; this bounds the tail well above that while preventing a
+# single draw from dominating national totals.
+MAX_SINGLE_GAIN = 500_000_000.0
+
+
+def _band_spline(row) -> UnivariateSpline:
+    """Linear spline through the published percentiles for one income band."""
+    return UnivariateSpline(
+        [0.05, 0.1, 0.25, 0.5, 0.75, 0.90, 0.95],
+        [row.p05, row.p10, row.p25, row.p50, row.p75, row.p90, row.p95],
+        k=1,
+    )
+
+
+def _body_mean(spline: UnivariateSpline) -> float:
+    """Mean of the spline body over quantiles [0, TAIL_START_QUANTILE)."""
+    # Deterministic quadrature rather than Monte Carlo, so the fitted tail does
+    # not depend on any random draw.
+    grid = np.linspace(0, TAIL_START_QUANTILE, 10001)
+    return float(np.trapezoid(spline(grid), grid) / TAIL_START_QUANTILE)
+
+
+def _capped_pareto_mean(alpha: float, scale: float, cap: float) -> float:
+    """E[min(X, cap)] for X ~ Pareto(shape=alpha, scale=scale)."""
+    if scale >= cap:
+        return cap
+    return scale + scale**alpha * (cap ** (1 - alpha) - scale ** (1 - alpha)) / (
+        1 - alpha
+    )
+
+
+def _fit_tail_alpha(row) -> float:
+    """Fit the Pareto shape parameter for one income band.
+
+    Chosen so that the band's full distribution -- spline body below p95, Pareto
+    above -- reproduces the published `mean_gains_given_gains` for that band.
+    """
+    from scipy.optimize import brentq
+
+    scale = float(row.p95)
+    target_mean = float(row.mean_gains_given_gains)
+    body_mean = _body_mean(_band_spline(row))
+
+    if scale <= 0:
+        # Degenerate band: no positive anchor for a multiplicative tail.
+        return MAX_TAIL_ALPHA
+
+    tail_weight = 1 - TAIL_START_QUANTILE
+    required_tail_mean = (target_mean - TAIL_START_QUANTILE * body_mean) / tail_weight
+
+    def gap(alpha: float) -> float:
+        return _capped_pareto_mean(alpha, scale, MAX_SINGLE_GAIN) - required_tail_mean
+
+    # The capped Pareto mean decreases monotonically in alpha, so the solution
+    # is bracketed whenever the endpoints straddle zero.
+    if gap(MIN_TAIL_ALPHA) < 0:
+        # Even the heaviest tail we allow cannot reach the published mean.
+        return MIN_TAIL_ALPHA
+    if gap(MAX_TAIL_ALPHA) > 0:
+        # Band mean is already met by the body; use the thinnest tail.
+        return MAX_TAIL_ALPHA
+    return float(brentq(gap, MIN_TAIL_ALPHA, MAX_TAIL_ALPHA))
+
+
+def sample_band_gains(row, quantiles: np.ndarray) -> np.ndarray:
+    """Map uniform quantiles to capital gains for one income band.
+
+    Below TAIL_START_QUANTILE this is exactly the original Advani & Summers
+    spline. Above it, quantiles are mapped through a fitted Pareto tail
+    anchored at the band's p95.
+    """
+    spline = _band_spline(row)
+    gains = spline(quantiles)
+
+    in_tail = quantiles >= TAIL_START_QUANTILE
+    if in_tail.any() and row.p95 > 0:
+        alpha = _fit_tail_alpha(row)
+        # Rescale the tail quantiles to (0, 1] and invert the Pareto CDF.
+        residual = (1 - quantiles[in_tail]) / (1 - TAIL_START_QUANTILE)
+        residual = np.clip(residual, 1e-12, 1.0)
+        gains[in_tail] = np.minimum(row.p95 * residual ** (-1 / alpha), MAX_SINGLE_GAIN)
+
+    return gains
+
+
 # Silence verbose logging
 logging.getLogger("root").setLevel(logging.WARNING)
 
@@ -124,17 +246,14 @@ def loss(blend_factor):
 
     for i in range(len(capital_gains)):
         row = capital_gains.iloc[i]
-        spline = UnivariateSpline(
-            [0.05, 0.1, 0.25, 0.5, 0.75, 0.90, 0.95],
-            [row.p05, row.p10, row.p25, row.p50, row.p75, row.p90, row.p95],
-            k=1,
-        )
         lower = row.minimum_total_income
         upper = row.maximum_total_income
         ti_in_range = (ti >= lower) * (ti < upper)
         in_target_range = has_cg * ti_in_range > 0
+        # One uniform draw per person from the seeded generator, exactly as
+        # before, so the sequence of random numbers consumed is unchanged.
         quantiles = cg_rng.random(int(in_target_range.sum()))
-        pred_capital_gains = spline(quantiles)
+        pred_capital_gains = sample_band_gains(row, quantiles)
         new_cg[in_target_range] = pred_capital_gains
 
     return new_cg, new_household_weight

policyengine_uk_data/storage/capital_gains_size_distribution_hmrc.csv

@@ -0,0 +1,39 @@
+# HMRC Capital Gains Tax statistics, Table 2.1a: estimated number of taxpayers,
+# amounts of gains and tax liabilities by size of gain, for individuals,
+# 2023 to 2024 tax year (provisional).
+#
+# Collection: https://www.gov.uk/government/statistics/capital-gains-tax-statistics
+# File:       Table_2_2025_Size_of_gain.ods, sheet "2_1a_2023-24"
+# URL:        https://assets.publishing.service.gov.uk/media/6878ac562bad77c3dae4dcef/Table_2_2025_Size_of_gain.ods
+# Published:  24 July 2025 release.
+#
+# Committed rather than downloaded at build time because the asset URL embeds a
+# content hash that changes with every HMRC release, so it is not a stable
+# fetch target. This follows the precedent in
+# policyengine_uk_data/targets/sources/ons_demographics.py.
+#
+# Units: taxpayers in thousands, gains in £ millions, as published.
+#
+# IMPORTANT DEFINITIONAL NOTE: HMRC's published figures "only include taxpayers
+# who have a CGT liability". They are therefore not the same population as
+# "everyone with gains above the annual exempt amount", and are narrower than
+# the set of people carrying a positive imputed gain in the microdata. The
+# under-£12,300 rows are correspondingly small and should not be read as the
+# true count of people with small gains.
+#
+# Trailing "All" row is the published total and may not equal the sum of the
+# bands due to rounding.
+lower_limit,taxpayers_thousands,gains_gbp_millions,tax_gbp_millions
+0,2,1,9
+3000,0,1,0
+6000,61,461,18
+10000,23,251,20
+12300,79,1418,169
+25000,74,2645,423
+50000,53,3731,706
+100000,37,5649,1077
+250000,14,4766,826
+500000,8,5705,895
+1000000,5,6390,1067
+2000000,3,9189,1718
+5000000,2,22714,4532

policyengine_uk_data/tests/test_capital_gains_tail.py

@@ -0,0 +1,165 @@
+"""Tests for the fitted Pareto upper tail of the capital gains imputation.
+
+The Advani & Summers percentile table stops at p95. Before this tail was
+added, sampling uniform quantiles through the linear spline extrapolated off
+the last segment, which removed the heavy tail of capital gains entirely.
+These tests pin the three properties that matter:
+
+* the body below p95 is still exactly the Advani & Summers spline;
+* the tail above p95 reproduces the concentration implied by the published
+  band means (and, in aggregate, HMRC's size-of-gain distribution);
+* sampling stays deterministic under the build's seed.
+"""
+
+import importlib
+
+import numpy as np
+import pytest
+
+# The imputations package re-exports a DataFrame named ``capital_gains``, which
+# shadows the submodule of the same name, so import it explicitly.
+cg_module = importlib.import_module(
+    "policyengine_uk_data.datasets.imputations.capital_gains"
+)
+
+TAIL_START_QUANTILE = cg_module.TAIL_START_QUANTILE
+capital_gains = cg_module.capital_gains
+sample_band_gains = cg_module.sample_band_gains
+
+
+@pytest.fixture(scope="module")
+def bands():
+    return [capital_gains.iloc[i] for i in range(len(capital_gains))]
+
+
+def test_body_below_p95_is_unchanged_spline(bands):
+    """Below p95 the sampler must be bit-for-bit the original spline."""
+    quantiles = np.linspace(0, TAIL_START_QUANTILE, 500, endpoint=False)
+    for row in bands:
+        expected = cg_module._band_spline(row)(quantiles)
+        actual = sample_band_gains(row, quantiles)
+        np.testing.assert_allclose(actual, expected, rtol=0, atol=0)
+
+
+def test_body_tracks_published_percentiles(bands):
+    """The body still follows the published percentiles.
+
+    Note this is a *smoothing* spline (scipy's default ``s``), not an
+    interpolating one, so it does not pass exactly through the published
+    points. That is pre-existing behaviour and is deliberately preserved here;
+    the tolerance below just pins that the body has not drifted.
+    """
+    published = np.array([0.05, 0.1, 0.25, 0.5, 0.75, 0.90])
+    for row in bands:
+        expected = np.array([row.p05, row.p10, row.p25, row.p50, row.p75, row.p90])
+        actual = sample_band_gains(row, published)
+        scale = max(abs(row.p90), 1.0)
+        assert np.max(np.abs(actual - expected)) < 0.5 * scale
+
+
+def test_tail_is_anchored_at_p95(bands):
+    """The tail starts at p95 and is monotonically increasing in the quantile."""
+    for row in bands:
+        if row.p95 <= 0:
+            continue
+        q = np.linspace(TAIL_START_QUANTILE, 1 - 1e-9, 2000)
+        tail = sample_band_gains(row, q)
+        assert tail[0] == pytest.approx(row.p95, rel=1e-6)
+        assert np.all(np.diff(tail) >= 0)
+        assert tail.max() <= cg_module.MAX_SINGLE_GAIN
+
+
+def test_tail_reproduces_published_band_means(bands):
+    """The fitted tail is what makes each band hit its published mean.
+
+    This is the calibration condition the shape parameter is solved from, so
+    it is the property that would break first if the fit regressed.
+    """
+    rng = np.random.default_rng(0)
+    floored = 0
+    for row in bands:
+        alpha = cg_module._fit_tail_alpha(row)
+        sampled = sample_band_gains(row, rng.random(200_000))
+        if alpha <= cg_module.MIN_TAIL_ALPHA:
+            # A small number of bands have an anomalously low p95 relative to
+            # their published mean; the alpha floor deliberately stops us
+            # chasing those with an unstably heavy tail, at the cost of
+            # undershooting the band mean.
+            floored += 1
+            continue
+        assert sampled.mean() == pytest.approx(row.mean_gains_given_gains, rel=0.15)
+    # If many bands start hitting the floor, the fit has regressed.
+    assert floored <= 3
+
+
+def test_spline_alone_undershoots_published_means(bands):
+    """Guards the diagnosis: the body alone cannot reach the published mean.
+
+    If this ever stops holding, the tail correction is no longer justified on
+    these grounds and the fit should be revisited.
+    """
+    rng = np.random.default_rng(0)
+    shortfalls = []
+    for row in bands:
+        body_only = cg_module._band_spline(row)(rng.random(100_000))
+        shortfalls.append(body_only.mean() / row.mean_gains_given_gains)
+    # Every band undershoots, typically by 35-50%.
+    assert max(shortfalls) < 0.95
+
+
+def test_fitted_alphas_are_in_a_plausible_heavy_tail_range(bands):
+    alphas = np.array([cg_module._fit_tail_alpha(row) for row in bands])
+    assert np.all(alphas >= cg_module.MIN_TAIL_ALPHA)
+    # Capital gains tails are heavy; a median far outside this range would
+    # indicate the moment condition is being solved against bad inputs.
+    assert 1.2 < np.median(alphas) < 2.0
+
+
+def test_tail_delivers_concentration(bands):
+    """Sampled gains must concentrate, which the truncated spline never did."""
+    rng = np.random.default_rng(0)
+    sampled = np.concatenate(
+        [sample_band_gains(row, rng.random(50_000)) for row in bands]
+    )
+    positive = sampled[sampled > 0]
+    share_1m = positive[positive >= 1e6].sum() / positive.sum()
+    # HMRC put ~61% of gains in disposals of £1m+. The imputation is fitted to
+    # Advani & Summers band means rather than to HMRC directly, so it lands
+    # below that, but it must be far above the ~2% the old spline produced.
+    assert share_1m > 0.25
+
+    # The old spline could not produce a gain above ~£2m at all.
+    assert positive.max() > 5e6
+
+
+def test_sampling_is_deterministic_under_the_seed(bands):
+    """The build seeds its generator because unseeded draws moved CGT revenue."""
+    row = bands[-1]
+
+    def draw():
+        rng = np.random.default_rng(0)
+        return sample_band_gains(row, rng.random(10_000))
+
+    np.testing.assert_array_equal(draw(), draw())
+
+
+def test_hmrc_reference_table_loads_and_is_concentrated():
+    """The committed HMRC reference table is the external validation anchor."""
+    import pandas as pd
+
+    from policyengine_uk_data.storage import STORAGE_FOLDER
+
+    table = pd.read_csv(
+        STORAGE_FOLDER / "capital_gains_size_distribution_hmrc.csv",
+        comment="#",
+    )
+    assert {
+        "lower_limit",
+        "taxpayers_thousands",
+        "gains_gbp_millions",
+    } <= set(table.columns)
+
+    total = table.gains_gbp_millions.sum()
+    above_1m = table.loc[table.lower_limit >= 1_000_000, "gains_gbp_millions"].sum()
+    # HMRC 2023-24: ~61% of gains accrue to disposals of £1m or more.
+    assert 0.55 < above_1m / total < 0.70