Tax bases at or above the last LinearAverageRateTaxScale threshold

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.

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:

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()

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,
+    )