diff --git a/changelog.d/fix-linear-average-rate-top-bracket.fixed.md b/changelog.d/fix-linear-average-rate-top-bracket.fixed.md new file mode 100644 index 000000000..f24a1440d --- /dev/null +++ b/changelog.d/fix-linear-average-rate-top-bracket.fixed.md @@ -0,0 +1 @@ +`LinearAverageRateTaxScale.calc` now taxes a base at or above the last threshold at the last average rate. It used to return 0 for such bases unless the last threshold was infinity, as `MarginalRateTaxScale.to_average` sets it. diff --git a/policyengine_core/taxscales/linear_average_rate_tax_scale.py b/policyengine_core/taxscales/linear_average_rate_tax_scale.py index cfe66e055..5180582d8 100644 --- a/policyengine_core/taxscales/linear_average_rate_tax_scale.py +++ b/policyengine_core/taxscales/linear_average_rate_tax_scale.py @@ -22,6 +22,28 @@ def calc( tax_base: NumericalArray, right: bool = False, ) -> numpy.float_: + """ + Compute the tax amount for the given tax bases by applying a taxscale. + + Each bracket's rate is the average rate at its threshold. Between two + thresholds the average rate moves linearly from one rate to the next. + From the last threshold up it stays at the last rate. Below the first + threshold the tax is 0, except that a scale with a single bracket + applies its rate to every tax base. + + :param ndarray tax_base: Array of the tax bases. + + :returns: Float array with tax amount for the given tax bases. + + For instance: + + >>> tax_scale = LinearAverageRateTaxScale() + >>> tax_scale.add_bracket(0, 0) + >>> tax_scale.add_bracket(100, 0.1) + >>> tax_base = array([50, 100, 150]) + >>> tax_scale.calc(tax_base) + [2.5, 10.0, 15.0] + """ if len(self.rates) == 1: return tax_base * self.rates[0] @@ -47,9 +69,15 @@ def calc( log.info(f"bracket_average_start_rate : {bracket_average_start_rate}") log.info(f"average_rate_slope: {average_rate_slope}") + # The brackets above end at the last threshold, so bases from there up + # fall in none of them. With no next rate to move toward, their average + # rate stays at the last rate. + above_last_threshold = numpy.asarray(tax_base) >= thresholds_array[-1] + return tax_base * ( +bracket_average_start_rate + (tax_base - bracket_threshold) * average_rate_slope + + above_last_threshold * rates_array[-1] ) def to_marginal(self) -> MarginalRateTaxScale: diff --git a/tests/core/tax_scales/test_linear_average_rate_tax_scale.py b/tests/core/tax_scales/test_linear_average_rate_tax_scale.py index 3e62b8987..dc711e4dd 100644 --- a/tests/core/tax_scales/test_linear_average_rate_tax_scale.py +++ b/tests/core/tax_scales/test_linear_average_rate_tax_scale.py @@ -1,7 +1,7 @@ import numpy import pytest -from policyengine_core import taxscales, tools +from policyengine_core import parameters, periods, taxscales, tools def test_bracket_indices(): @@ -69,6 +69,93 @@ def test_to_dict(): assert result == {"0": 0.0, "100": 0.1} +def test_calc(): + tax_base = numpy.array([0.5, 1, 1.5, 2, 2.5, 100]) + tax_scale = taxscales.LinearAverageRateTaxScale() + tax_scale.add_bracket(0, 0) + tax_scale.add_bracket(1, 0.1) + tax_scale.add_bracket(2, 0.2) + + result = tax_scale.calc(tax_base) + + # From the last threshold up the average rate stays at 0.2. These bases + # used to get 0 because they fall in no bracket. + tools.assert_near( + result, [0.025, 0.1, 0.225, 0.4, 0.5, 20], absolute_error_margin=1e-10 + ) + + +def test_calc_below_first_threshold(): + tax_base = numpy.array([-5, 0, 5, 10, 15, 20, 30]) + tax_scale = taxscales.LinearAverageRateTaxScale() + tax_scale.add_bracket(10, 0.1) + tax_scale.add_bracket(20, 0.2) + + result = tax_scale.calc(tax_base) + + tools.assert_near(result, [0, 0, 0, 1, 2.25, 4, 6], absolute_error_margin=1e-10) + + +@pytest.mark.parametrize("top_rate", [0, 0.2, 0.9]) +def test_calc_ignores_an_infinite_last_threshold(top_rate): + # ``MarginalRateTaxScale.to_average`` ends scales with an infinite + # threshold. The bracket it closes then has a slope of 0, which is what a + # scale without it does from its last finite threshold up. + tax_base = numpy.array([0.5, 1, 1.5, 2, 2.5, 100, 1e12]) + tax_scale = taxscales.LinearAverageRateTaxScale() + tax_scale.add_bracket(0, 0) + tax_scale.add_bracket(1, 0.1) + tax_scale.add_bracket(2, 0.2) + closed_tax_scale = tax_scale.copy() + closed_tax_scale.add_bracket(numpy.inf, top_rate) + + tools.assert_near( + closed_tax_scale.calc(tax_base), + tax_scale.calc(tax_base), + absolute_error_margin=1e-10, + ) + + +def test_calc_agrees_with_to_marginal_from_last_threshold_up(): + tax_base = numpy.array([2, 2.5, 100, 1e9]) + tax_scale = taxscales.LinearAverageRateTaxScale() + tax_scale.add_bracket(0, 0) + tax_scale.add_bracket(1, 0.1) + tax_scale.add_bracket(2, 0.2) + + result = tax_scale.calc(tax_base) + + tools.assert_near( + result, tax_scale.to_marginal().calc(tax_base), absolute_error_margin=1e-6 + ) + + +def test_calc_of_average_rate_parameter_scale(): + # A parameter scale with ``average_rate`` brackets loads as a + # LinearAverageRateTaxScale, with no infinite last threshold. + data = { + "description": "Average rate income tax", + "brackets": [ + { + "threshold": {"2020-01-01": {"value": threshold}}, + "average_rate": {"2020-01-01": {"value": rate}}, + } + for threshold, rate in [ + (1_100, 0.02), + (10_000, 0.25), + (100_000, 0.6), + ] + ], + } + scale = parameters.ParameterScale("average_rate_scale", data, "") + tax_scale = scale.get_at_instant(periods.Instant((2020, 6, 1))) + + result = tax_scale.calc(numpy.array([500, 1_100, 100_000, 250_000])) + + assert isinstance(tax_scale, taxscales.LinearAverageRateTaxScale) + tools.assert_near(result, [0, 22, 60_000, 150_000], absolute_error_margin=1e-6) + + def test_to_marginal(): tax_base = numpy.array([1, 1.5, 2, 2.5]) tax_scale = taxscales.LinearAverageRateTaxScale() diff --git a/tests/core/tax_scales/test_linear_average_rate_tax_scale_property.py b/tests/core/tax_scales/test_linear_average_rate_tax_scale_property.py new file mode 100644 index 000000000..ee2342ad1 --- /dev/null +++ b/tests/core/tax_scales/test_linear_average_rate_tax_scale_property.py @@ -0,0 +1,176 @@ +"""Properties of ``LinearAverageRateTaxScale.calc``. + +For any scale with 1 to 6 brackets and finite thresholds, and any tax base: + +* from the first threshold up, the tax is the base times the rate that + ``numpy.interp`` gives: linear between thresholds, and the last rate from the + last threshold up. Below the first threshold the tax is 0, except that a + scale with one bracket applies its rate to every base; +* at each threshold the tax is the threshold times its rate, and above the + first threshold the tax has no jump there; +* adding a last threshold of infinity, at any rate, does not change the tax; +* when the first threshold is 0, ``to_marginal`` levies the same tax at each + threshold and above the last one, within accumulated floating-point error. + +Examples are in ``test_linear_average_rate_tax_scale.py``. +""" + +from __future__ import annotations + +import numpy +import pytest + +# The smoke job installs Core without the dev extra but collects every module. +pytest.importorskip("hypothesis") + +from hypothesis import example, given, settings # noqa: E402 +from hypothesis import strategies as st # noqa: E402 + +from policyengine_core import taxscales # noqa: E402 + +RATES = st.integers(min_value=-1_000, max_value=1_000).map(lambda r: r / 1_000) +BASES = st.lists( + st.floats(min_value=-1e7, max_value=1e7, allow_nan=False), + max_size=20, +) + + +@st.composite +def _brackets( + draw, first_threshold=st.integers(min_value=0, max_value=1_000), min_size=1 +): + first = draw(first_threshold) + upper_thresholds = draw( + st.lists( + st.integers(min_value=first + 1, max_value=1_000_000), + min_size=min_size - 1, + max_size=5, + unique=True, + ) + ) + thresholds = [first] + sorted(upper_thresholds) + rates = draw(st.lists(RATES, min_size=len(thresholds), max_size=len(thresholds))) + return thresholds, rates + + +def _scale(thresholds, rates): + tax_scale = taxscales.LinearAverageRateTaxScale() + for threshold, rate in zip(thresholds, rates): + tax_scale.add_bracket(threshold, rate) + return tax_scale + + +def _bases(thresholds, extra_bases): + # Every threshold, its neighbours on both sides, and points well above + # the last threshold, where bases used to get 0. + near = [ + numpy.nextafter(threshold, direction) + for threshold in thresholds + for direction in (-numpy.inf, numpy.inf) + ] + far = [thresholds[-1] + 1, 2 * thresholds[-1] + 1, 1e9] + return numpy.array(thresholds + near + far + extra_bases, dtype=float) + + +@settings(max_examples=300, deadline=None) +@given(_brackets(), BASES) +def test_calc_interpolates_average_rate(brackets, extra_bases): + thresholds, rates = brackets + tax_base = _bases(thresholds, extra_bases) + + result = _scale(thresholds, rates).calc(tax_base) + + expected = tax_base * numpy.interp(tax_base, thresholds, rates) + if len(thresholds) > 1: + expected = numpy.where(tax_base >= thresholds[0], expected, 0) + numpy.testing.assert_allclose(result, expected, rtol=1e-9, atol=1e-6) + + +@settings(max_examples=300, deadline=None) +@given(_brackets()) +@example(brackets=([0, 999_999, 1_000_000], [0, 1, 0])) +def test_calc_is_continuous_at_thresholds(brackets): + thresholds, rates = brackets + tax_scale = _scale(thresholds, rates) + at_threshold = numpy.array(thresholds, dtype=float) + below_threshold = numpy.nextafter(at_threshold, -numpy.inf) + + result = tax_scale.calc(at_threshold) + + numpy.testing.assert_allclose( + result, at_threshold * numpy.array(rates), rtol=1e-9, atol=1e-6 + ) + # Below the first threshold is outside the scale, so start at the second. + rates_array = numpy.array(rates) + slope = numpy.diff(rates_array) / numpy.diff(at_threshold) + below = below_threshold[1:] + distance = at_threshold[1:] - below + # For T(x) = x * r(x), T(t) - T(b) = (t-b) * (r(t) + b*slope). + # The explicit example has distance 2**-33 and slope -1, so its tax + # movement is (1_000_000 - 2**-33) * 2**-33 = 0.0001164153218269. + movement = distance * (numpy.abs(rates_array[1:]) + numpy.abs(below * slope)) + # Budget interpolation arithmetic using its terms, including cancellation + # near a zero rate, rather than only the resulting tax's magnitude. + rounding = ( + 8 + * numpy.finfo(numpy.float64).eps + * at_threshold[1:] + * (numpy.abs(rates_array[:-1]) + numpy.abs((below - at_threshold[:-1]) * slope)) + ) + numpy.testing.assert_array_less( + numpy.abs(tax_scale.calc(below_threshold)[1:] - result[1:]), + 1e-6 + movement + rounding, + ) + + +@settings(max_examples=300, deadline=None) +@given(_brackets(), RATES, BASES) +def test_infinite_last_threshold_changes_nothing(brackets, top_rate, extra_bases): + thresholds, rates = brackets + tax_base = _bases(thresholds, extra_bases) + + if len(thresholds) == 1: + # One bracket applies its rate below its threshold too; two do not. + tax_base = tax_base[tax_base >= thresholds[0]] + + open_scale = _scale(thresholds, rates) + closed_scale = _scale(thresholds + [numpy.inf], rates + [top_rate]) + + numpy.testing.assert_array_equal( + closed_scale.calc(tax_base), open_scale.calc(tax_base) + ) + + +@settings(max_examples=300, deadline=None) +@given(_brackets(first_threshold=st.just(0), min_size=2)) +@example(brackets=([0, 999_999, 1_000_000], [0, 1, 0])) +def test_calc_agrees_with_to_marginal_at_thresholds_and_above(brackets): + # One bracket is left out: ``to_marginal`` cannot convert it yet. + thresholds, rates = brackets + tax_scale = _scale(thresholds, rates) + tax_base = numpy.array( + thresholds + [thresholds[-1] + 1, 2 * thresholds[-1] + 1, 1e9], dtype=float + ) + + marginal_scale = tax_scale.to_marginal() + lower = numpy.array(marginal_scale.thresholds, dtype=float) + upper = numpy.append(lower[1:], numpy.inf) + # Marginal calc shifts thresholds by 1 + eps. Each active bracket's + # endpoint magnitudes bound this shift and width subtraction; weighting + # by |marginal rate| also covers cancellation between large bracket taxes. + endpoints = lower + numpy.minimum(tax_base[:, None], upper) + endpoints = numpy.where(tax_base[:, None] >= lower, endpoints, 0) + # The stored threshold*rate taxes telescope in exact arithmetic. Allow + # rounding in secant conversion, widths, products and the bracket sum. + rounding = ( + (len(lower) + 8) + * numpy.finfo(numpy.float64).eps + * numpy.dot(endpoints, numpy.abs(marginal_scale.rates)) + ) + # In the explicit example marginal rates are [1, -999_999, 0]. Shifting + # 999_999 by 2**-32 leaves residual tax 1_000_000 * 2**-32 + # = 0.0002328306436539 at the final threshold, whose average tax is 0. + numpy.testing.assert_array_less( + numpy.abs(tax_scale.calc(tax_base) - marginal_scale.calc(tax_base)), + 1e-6 + rounding, + )