More robustness to rate

src/functions_rate.jl

@@ -1,35 +1,45 @@
 function rate(
-    nper, pmt, pv, fv=0.0, type=0, guess=nothing; tol=1e-11, max_iter=200
+    nper, pmt, pv, fv=0.0, type=0, guess=nothing; tol=1e-11, max_iter=150
 )::Float64
     if type != 0 && type != 1
         error("type must be 0 (end of period) or 1 (beginning of period)")
     end
 
-    rate = guess === nothing ? _heuristic_rate_guess(nper, pmt, pv, fv) : guess
+    fv == pv && return NaN
 
-    for i in 1:max_iter
-        new_rate = _f_rate_newton_step(rate, nper, pmt, pv, fv, type)
+    current_rate = if guess !== nothing
+        guess
+    else
+        _rate_guess(nper, pmt, pv, fv)
+    end
+
+    for _ in 1:max_iter
+        new_rate = _f_rate_newton_step(current_rate, nper, pmt, pv, fv, type)
 
-        if abs(new_rate - rate) < tol
+        if abs(new_rate - current_rate) < tol
             return new_rate
         end
 
-        rate = new_rate
+        current_rate = new_rate
     end
 
-    println(max_iter)
-    return rate
-end
-
-function _heuristic_rate_guess(nper, pmt, pv, fv)
-    if pv - fv == 0
-        return 0.0
+    return if abs(_pv(current_rate, nper, pmt, fv, type) - pv) / max(1, abs(pv)) > 1e-6
+        NaN
     else
-        return (pmt * nper + (pv + fv)) / (-(pv + fv) * nper)
+        # pv error less than 6 significant digits
+        current_rate
     end
 end
 
-function _f_rate(rate, nper, pmt, pv, fv=0.0, type=0)
+function _rate_guess(nper, pmt, pv, fv)::Float64
+    return max(
+        -Inf,
+        (pmt * nper + (pv + fv)) / (-(pv + fv) * nper) *
+        (1 + (pv - fv) / (pv + fv) * 1 / (2 * nper)),
+    )
+end
+
+function _f_rate(rate, nper, pmt, pv, fv, type)::Float64
     if rate == 0.0
         return pv + pmt * nper + fv
     else
@@ -39,7 +49,17 @@ function _f_rate(rate, nper, pmt, pv, fv=0.0, type=0)
     end
 end
 
-function _f_rate_derivative(rate, nper, pmt, pv, fv=0.0, type=0)
+function _pv(rate, nper, pmt, fv, type)::Float64
+    if rate == 0.0
+        return -(pmt * nper + fv)
+    else
+        return -(
+            pmt * (1 + rate * type) * (1 - (1 + rate)^(-nper)) / rate + fv / (1 + rate)^nper
+        )
+    end
+end
+
+function _f_rate_derivative(rate, nper, pmt, pv, fv, type)::Float64
     if rate == 0.0
         return pmt * nper
     else
@@ -49,8 +69,8 @@ function _f_rate_derivative(rate, nper, pmt, pv, fv=0.0, type=0)
     end
 end
 
-function _f_rate_newton_step(rate_guess, nper, pmt, pv, fv, type)
+function _f_rate_newton_step(rate_guess, nper, pmt, pv, fv, type)::Float64
     f_value = _f_rate(rate_guess, nper, pmt, pv, fv, type)
     f_prime = _f_rate_derivative(rate_guess, nper, pmt, pv, fv, type)
-    return rate_guess - f_value / f_prime
+    return max(-Inf, rate_guess - f_value / f_prime)
 end

test/runtests.jl

@@ -547,7 +547,6 @@ end
 end
 
 @testitem "rate" begin
-    # Comparison to Excel
     @test rate(60, -200, 10000, 0) ≈ 0.00618341316125379
     @test rate(120, -300, 15000, 5000) ≈ 0.0163743926512608
     @test rate(48, -250, 8000, 1000) ≈ 0.0151344680950013
@@ -577,6 +576,7 @@ end
     @test rate(180, -2000, 50000, -5000) ≈ 0.0399689060288228
     @test rate(24, 500, -12000, 3000) ≈ 0.0156905080913063
     @test rate(12, 100, -1000, -200) ≈ 0 #-1.29803547364284E-09
+    @test rate(24, -500, 12000, 0) ≈ 0
     @test rate(72, -300, 15000, -3000, 0, 0.08) ≈ 0.0133056860905425
     @test rate(120, 400, -25000, 10000, 1, 0.03) ≈ 0.0143238026598161
     @test rate(84, -350, 12000, -2000, 0, 0.07) ≈ 0.0263888939219417
@@ -587,4 +587,5 @@ end
     @test rate(360, 1200, -250000, 50000, 0, 0.04) ≈ 0.00374547917992404
     @test rate(240, -800, 150000, -10000, 1, 0.05) ≈ 0.00250603712857146
     @test rate(12, 300, -3500, 200, 0, 0.1) ≈ 0.0123295838958103
+    @test rate(12, 300, 200, 200, 0, 0.1) === NaN
 end