BUG: restrict rate results to values above -100%

numpy_financial/_financial.py

@@ -891,6 +891,92 @@ def _g_div_gp(r, n, p, x, y, w) -> Any:
     return g / gp
 
 
+def _rate_residual(r, n, p, x, y, w):
+    """Evaluate the rate equation continuously at zero."""
+    if r == 0:
+        return y + x + p * n
+    if isinstance(r, Decimal):
+        growth = (1 + r) ** n
+        annuity = (growth - 1) / r
+    else:
+        exponent = n * np.log1p(r)
+        growth = np.exp(exponent)
+        annuity = np.expm1(exponent) / r
+    return y + x * growth + p * (1 + r * w) * annuity
+
+
+def _rate_bisect(n, p, x, y, w, guess, tol, maxiter):
+    """Find a sign-changing bracket above -1, then bisect it.
+
+    Search outwards from the guess. This is a fallback, not an exhaustive
+    search for all roots; a root without a sign change may be missed.
+    """
+    decimal = isinstance(p, Decimal)
+    number = Decimal if decimal else float
+    nan = number('NaN')
+    center: Any = number(guess)
+    if not np.isfinite(float(center)) or center <= -1:
+        return nan
+    # The rate is independent of the currency unit. Normalization avoids
+    # rejecting finite Decimal cash flows outside the float exponent range.
+    cash_scale = max(abs(p), abs(x), abs(y))
+    if cash_scale != 0:
+        p, x, y = p / cash_scale, x / cash_scale, y / cash_scale
+    left = right = center
+    fl = fr = _rate_residual(center, n, p, x, y, w)
+    if fl == 0:
+        return center
+    bracket = False
+    for _ in range(maxiter):
+        if not bracket:
+            # Geometric expansion in the positive compounding factor 1+r.
+            new_left = (left - 1) / 2
+            new_right = 2 * right + 1
+            # Exhausting one direction must not stop the other direction.
+            if new_left <= -1 or new_left == left:
+                new_left = left
+            if not np.isfinite(float(new_right)) or new_right == right:
+                new_right = right
+            if new_left == left and new_right == right:
+                return nan
+            fleft = _rate_residual(new_left, n, p, x, y, w)
+            fright = _rate_residual(new_right, n, p, x, y, w)
+            if fleft == 0:
+                return new_left
+            if fright == 0:
+                return new_right
+            if (fleft < 0 < fl) or (fl < 0 < fleft):
+                right, fr = left, fl
+                left, fl = new_left, fleft
+                bracket = True
+            elif (fr < 0 < fright) or (fright < 0 < fr):
+                left, fl = right, fr
+                right, fr = new_right, fright
+                bracket = True
+            else:
+                left, fl, right, fr = new_left, fleft, new_right, fright
+            continue
+        mid = (left + right) / 2
+        fm = _rate_residual(mid, n, p, x, y, w)
+        if not np.isfinite(float(fm)):
+            return nan
+        if fm == 0:
+            return mid
+        if right - left <= tol or mid == left or mid == right:
+            scale = abs(y) + abs(x * (1 + mid) ** n)
+            scale += abs(p * (1 + mid * w) * (
+                n if mid == 0 else ((1 + mid) ** n - 1) / mid))
+            if (np.isfinite(float(scale))
+                    and abs(fm) <= tol * max(number(1), scale)):
+                return mid
+            return nan
+        if (fl < 0 < fm) or (fm < 0 < fl):
+            right, fr = mid, fm
+        else:
+            left, fl = mid, fm
+    return nan
+
+
 # Use Newton's iteration until the change is less than 1e-6
 #  for all values or a maximum of 100 iterations is reached.
 #  Newton's rule is
@@ -925,7 +1011,8 @@ def rate(
     when : {{'begin', 1}, {'end', 0}}, {string, int}, optional
         When payments are due ('begin' (1) or 'end' (0))
     guess : Number, optional
-        Starting guess for solving the rate of interest, default 0.1
+        Starting guess for solving the rate of interest, default 0.1.
+        Must be finite and greater than -1.
     tol : Number, optional
         Required tolerance for the solution, default 1e-6
     maxiter : int, optional
@@ -943,7 +1030,10 @@ def rate(
 
      fv + pv*(1+rate)**nper + pmt*(1+rate*when)/rate * ((1+rate)**nper - 1) = 0
 
-    for ``rate``.
+    for ``rate`` in the domain ``rate > -1``. If Newton's method finds a
+    root outside this domain, a sign-changing bracket is sought around
+    ``guess`` and solved by bisection. This search does not find all roots
+    and cannot guarantee the root nearest to ``guess``.
 
     References
     ----------
@@ -969,16 +1059,40 @@ def rate(
 
     nper, pmt, pv, fv, when = map(np.asarray, [nper, pmt, pv, fv, when])
 
+    if not np.isfinite(float(guess)) or guess <= -1:
+        raise ValueError("guess must be finite and greater than -1")
+
     rn: Any = guess
     iterator = 0
     close: Any = False
+    converged_iteration: Any = 0
     while (iterator < maxiter) and not np.all(close):
         rnp1 = rn - _g_div_gp(rn, nper, pmt, pv, fv, when)
         diff = abs(rnp1 - rn)
         close = diff < tol
         iterator += 1
+        converged_iteration = np.where(
+            close & (converged_iteration == 0), iterator, converged_iteration)
         rn = rnp1
 
+    # A converged algebraic root below -1 is not a financial rate.
+    invalid = np.asarray(rn <= -1) & np.asarray(close)
+    if np.any(invalid):
+        rn = np.asarray(rn).copy()
+        close = np.asarray(close).copy()
+        inputs = np.broadcast_arrays(nper, pmt, pv, fv, when, rn)
+        for index in np.ndindex(rn.shape):
+            if invalid[index]:
+                n, p, x, y, w = (values[index] for values in inputs[:5])
+                with np.errstate(over="ignore", invalid="ignore"):
+                    result = _rate_bisect(
+                        n, p, x, y, w, guess, tol,
+                        max(0, maxiter - int(converged_iteration[index])))
+                rn[index] = result
+                close[index] = np.isfinite(float(result))
+        if rn.ndim == 0:
+            rn = rn.item()
+
     if not np.all(close):
         if np.isscalar(rn):
             if raise_exceptions:

numpy_financial/tests/test_financial.py

@@ -170,6 +170,123 @@ def test_pv_decimal(self):
 
 
 class TestRate:
+    @pytest.mark.parametrize("number_type", [float, Decimal])
+    @pytest.mark.parametrize("pmt,pv,expected", [
+        (-440000, 263175, 1.6711838275594646),
+        (263175, -440000, 0.5838779110248231),
+    ])
+    def test_gh126(self, number_type, pmt, pv, expected):
+        args = [number_type(v) for v in (8, pmt, pv, 25500)]
+        result = npf.rate(*args, tol=number_type("1e-10"))
+        assert result > -1
+        assert_allclose(float(result), expected, atol=1e-10, rtol=0)
+        if number_type is Decimal:
+            assert isinstance(result, Decimal)
+        # Independently discount the actual cash flows at the returned rate.
+        cashflows = [pv] + [pmt] * 7 + [pmt + 25500]
+        residual = sum(c / (1 + float(result)) ** t
+                       for t, c in enumerate(cashflows))
+        assert abs(residual) < 1e-3
+
+    def test_gh126_broadcast(self):
+        pmt = numpy.array([[-440000], [263175]])
+        pv = numpy.array([[263175], [-440000]])
+        result = npf.rate(8, pmt, pv, [25500, 25500], tol=1e-10)
+        assert_allclose(result, [[1.6711838275594646] * 2,
+                                 [0.5838779110248231] * 2], atol=1e-10)
+
+    def test_gh126_mixed_results(self):
+        result = npf.rate([8, 2, 2], [263175, 0, 0],
+                          [-440000, -100, 100], [25500, 81, 81])
+        assert_allclose(result, [0.583877911, -0.1, numpy.nan], atol=1e-6)
+
+    @pytest.mark.parametrize("when", [0, 1])
+    def test_rate_repayment_round_trip(self, when):
+        target = 0.6
+        pmt = npf.pmt(target, 8, -440000, 25500, when=when)
+        result = npf.rate(8, pmt, -440000, 25500, when=when, tol=1e-10)
+        assert_allclose(result, target, atol=1e-10)
+
+    @pytest.mark.parametrize("guess", [-1, -2, numpy.nan, numpy.inf])
+    def test_rate_invalid_guess(self, guess):
+        with pytest.raises(ValueError, match="guess"):
+            npf.rate(8, 263175, -440000, 25500, guess=guess)
+
+    def test_gh126_iteration_budget(self):
+        assert numpy.isnan(npf.rate(8, 263175, -440000, 25500, maxiter=10))
+        with pytest.raises(npf.IterationsExceededError):
+            npf.rate(8, 263175, -440000, 25500, maxiter=10,
+                     raise_exceptions=True)
+
+    @pytest.mark.parametrize("scale", ["1e-400", "1", "1e400"])
+    def test_gh126_decimal_currency_scale(self, scale):
+        factor = Decimal(scale)
+        result = npf.rate(Decimal(8), Decimal(263175) * factor,
+                          Decimal(-440000) * factor, Decimal(25500) * factor,
+                          tol=Decimal("1e-10"))
+        assert isinstance(result, Decimal)
+        assert_allclose(float(result), 0.5838779110248231,
+                        atol=1e-10, rtol=0)
+
+    @pytest.mark.parametrize("scale", [1e-200, 1.0, 1e200])
+    def test_gh126_float_currency_scale(self, scale):
+        result = npf.rate(8, 263175 * scale, -440000 * scale,
+                          25500 * scale, tol=1e-10)
+        assert_allclose(result, 0.5838779110248231, atol=1e-10, rtol=0)
+
+    def test_rate_fallback_large_root(self):
+        from numpy_financial._financial import _rate_bisect
+
+        # Analytic one-period growth. Reaching the lower domain boundary
+        # must not prevent discovery of a much larger positive root.
+        growth = 2.0 ** 60
+        result = _rate_bisect(1, 0, -1, growth, 0, 0.1, 1e-10, 200)
+        assert_allclose(result, growth - 1, rtol=1e-12)
+
+    def test_rate_fallback_near_lower_boundary(self):
+        from numpy_financial._financial import _rate_bisect
+
+        result = _rate_bisect(1, 0, -1, 1e-8, 0, 0.1, 1e-14, 200)
+        assert result > -1
+        assert_allclose(result, -1 + 1e-8, atol=1e-14, rtol=0)
+
+    def test_rate_fallback_without_sign_change(self):
+        from numpy_financial._financial import _rate_bisect
+
+        # Positive cash flows cannot cancel anywhere in r > -1.
+        assert numpy.isnan(_rate_bisect(1, 0, 1, 1, 0, 0.1, 1e-10, 100))
+
+    @pytest.mark.parametrize("principal,nper,annual,balloon", [
+        (250000, 300, "0.04", 0),
+        (250000, 300, "0.06", 0),
+        (15000, 48, "0.08", 0),
+        (20000, 48, "0.07", 8000),
+    ])
+    @pytest.mark.parametrize("when", [0, 1])
+    def test_rate_practitioner_loans(self, principal, nper, annual, balloon, when):
+        from decimal import localcontext
+
+        # Independent high-precision annuity calculation, without npf.pmt.
+        # Annual rates are nominal, compounded monthly; payments are unrounded.
+        with localcontext() as context:
+            context.prec = 60
+            monthly = Decimal(annual) / 12
+            growth = (1 + monthly) ** nper
+            payment = ((principal * growth - balloon) * monthly
+                       / ((growth - 1) * (1 + monthly * when)))
+        result = npf.rate(nper, -float(payment), principal, -balloon,
+                          when=when, tol=1e-10)
+        assert_allclose(result, float(monthly), atol=1e-10, rtol=0)
+        # Check the actual month-by-month outstanding balance independently.
+        balance = float(principal)
+        for _ in range(nper):
+            if when:
+                balance -= float(payment)
+            balance *= 1 + result
+            if not when:
+                balance -= float(payment)
+        assert_allclose(balance, balloon, atol=0.01, rtol=0)
+
     def test_rate(self):
         assert_allclose(npf.rate(10, 0, -3500, 10000), 0.1107, rtol=1e-4)