TVM solve for interest rate, revisited

[Albert Chan]

Instead of cheap tricks to get good rate guess, why not start from the edge, and iterate with Halley's method ?

For NPMT=0 with 2 rate solutions, this picked the "smaller" edge.
Note: NaN edge moved up front, for rate search early termination.

Code:

function edge_i(n,pv,pmt,fv)
    local a, b = pmt/fv, -pmt/pv
    if abs(a)>abs(b) and b>-1 or b~=b then a,b = b,a end
    return a, b  -- a = small edge or NaN
end

We solve NPMT=0, with PV=0, for rate x

x = loan_rate(n, pv, pmt, fv) = loan_rate(n, 0, pv*x+pmt, pv+fv)

f = (pv+fv)/n * g + pv * x + pmt = 0      // same formula used for Plus42

g = n*x / ((1+x)^n-1)
g'/g = -(g-1+(n-1)*x) / (x+x^2) = -num / den
g''/g = (g+n*x)*(2*(g-1)+(n-1)*x) / (x+x^2)^2 = (num+x+1)*(num+(g-1)) / den^2

Let k = (pv+fv)/n * g

f = k + pv*x + pmt
f' = k*(g'/g) + pv
f'' = k*(g''/g)

Using first derivative to get to second, Halley's correction is cheap.
To make code more robust, we special cased with limits when rate goes 0

Code:

function loan_rate2(n, pv, pmt, fv, i)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n, pv, pmt, fv)
    return function()
        if i==0 then
            local a, b = (pv+fv)/n, pv-fv
            local f, fp = (a+pmt), (a+b)*0.5
            i = -f*fp/(fp*fp-f*(n*n-1)/12*a)
            return i, i
        end
        local x = i/expm1(n*log1p(i))
        local k, y = (pv+fv)*x, n*x-1
        local f = k + pv*i + pmt
        local num, den = y+(n-1)*i, i+i*i
        x = k*num - pv*den      -- f' * -den
        y = k*(num+i+1)*(num+y) -- f'' * den^2
        x = f*x*den / (x*x-.5*f*y)
        i = i + x   -- Halley's method
        return i, x
    end
end

Previous post example:

lua> n, pv, pmt, fv = 10, 50, -30, 100
lua> g = loan_rate2(n,pv,pmt,fv)
lua> g()
-0.2844446981069045        0.015555301893095532
lua> g()
-0.28443599888025756      8.699226646944535e-006
lua> g()
-0.28443599888025595      1.6279605736411871e-015

Fun Math Algorithms, car lease APR example.

lua> n, pv, pmt, fv = 36, 30000, -550, -15000
lua> edge_i(n,pv,pmt,fv) -- Note: 2nd edge can be removed
0.018333333333333333       0.03666666666666667
lua>
lua> g = loan_rate2(n,pv,pmt,fv)
lua> g()
0.005817405851432355      -0.012515927481900979
lua> g()
0.005805072819430545      -1.2333032001809594e-005
lua> g()
0.0058050728194201295    -1.0415537616188019e-014

Note sign of errors.
Iteration does not over-shoot to the other side.
Next iteration guaranteed better estimate (cubic convergence)

Bonus: with one-sided convergence pattern, we can detect for no solution.
If error sign changes is not due to simple rounding errors, no solution.

lua> g = loan_rate2(10,50,-10,100)
lua> g()
0.3175682256333693      0.4175682256333693
lua> g()
-0.0740994732881643     -0.3916676989215336 --> no solution