From 3c6157bf55acba4d6987310aa7e1294eeac05ff5 Mon Sep 17 00:00:00 2001 From: Vahid Ahmadi Date: Mon, 20 Jul 2026 20:13:51 +0100 Subject: [PATCH] Fit a Pareto upper tail to the capital gains imputation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Advani & Summers percentile table stops at p95, so fitting a linear spline through p05..p95 and sampling uniform quantiles extrapolated everything above the 95th percentile off the last segment. That left no Pareto tail: in the published enhanced_frs_2024_25.h5 only 1.6% of gains came from gains of £1m or more (HMRC: ~61%) and no gain exceeded ~£2m. Gains above p95 are now drawn from a Pareto tail anchored at each income band's p95. The shape parameter is fitted per band, not assumed: the same source table publishes mean_gains_given_gains, which the spline alone undershoots by 35-50% in every band, and that gap is the missing tail mass. We solve for the shape parameter that makes body + tail reproduce the published band mean. Fitted values cluster around 1.46. Fitting per income band preserves the existing income-band conditioning. The seeded generator and the spline body below p95 are unchanged. Adds HMRC Capital Gains Tax statistics Table 2.1a as a committed validation reference. Co-Authored-By: Claude Opus 4.8 (1M context) --- changelog.d/443.md | 1 + .../datasets/imputations/capital_gains.py | 131 +++++++++++++- .../capital_gains_size_distribution_hmrc.csv | 39 +++++ .../tests/test_capital_gains_tail.py | 165 ++++++++++++++++++ 4 files changed, 330 insertions(+), 6 deletions(-) create mode 100644 changelog.d/443.md create mode 100644 policyengine_uk_data/storage/capital_gains_size_distribution_hmrc.csv create mode 100644 policyengine_uk_data/tests/test_capital_gains_tail.py diff --git a/changelog.d/443.md b/changelog.d/443.md new file mode 100644 index 00000000..d0aa3a13 --- /dev/null +++ b/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. diff --git a/policyengine_uk_data/datasets/imputations/capital_gains.py b/policyengine_uk_data/datasets/imputations/capital_gains.py index 22fe38d0..10c15d5c 100644 --- a/policyengine_uk_data/datasets/imputations/capital_gains.py +++ b/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 diff --git a/policyengine_uk_data/storage/capital_gains_size_distribution_hmrc.csv b/policyengine_uk_data/storage/capital_gains_size_distribution_hmrc.csv new file mode 100644 index 00000000..804b8669 --- /dev/null +++ b/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 diff --git a/policyengine_uk_data/tests/test_capital_gains_tail.py b/policyengine_uk_data/tests/test_capital_gains_tail.py new file mode 100644 index 00000000..6c42dab3 --- /dev/null +++ b/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