TVM solve for interest rate, revisited

[Albert Chan]

(06-11-2022 05:34 PM)Albert Chan Wrote:  Rational Halley's modified slope = f' * (1 - (f*f'')/(f')^2/2) = f' - (f''/2)*(f/f')

With slope adjustment, (f''/2)*(f/f'), getting bigger in size relative to f',
it may turn slope to the wrong direction! Or, worse ... slope of 0!

Below example, with i0 = pmt/fv = 1e-9, 1st iteration get worse (turned negative!)
True rate = 5.32155425915678; Rate iteration should not go backward.

lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> loan_rate2(n,pv,pmt,fv)() + 0
-0.08736585216654838

Checking f and its derivative at i=i0, Halley's slope had wrong sign.

lua> f0, f1, f2 = 999999995.5, -4500000102.007542, 143515080310.56577
lua> f1 - (f2/2)*(f0/f1) -- bad slope
11446119499.270123

Technically, with guess rate from the edge, this should not happen.
The reason is a badly estimated f'', over-estimated almost 9 times!
With more precision, this is true f0, f1, f2, at i=i0

lua> f0, f1, f2 = 999999995.5, -4500000038.5, 16499999810.25
lua> f1 - (f2/2)*(f0/f1) -- true halley's slope
-2666666750.1851845
lua> 1e-9 - f0/_ -- this should be i1
0.37499998756770886

Instead of special cased at *exactly* i=0, we should widen it.
Below code treat eps=i*i, such that 1+eps==1, as special.

From above examples, slope is less affected by catastrophic cancellations.
But, I fix it anyway for Newton's method too (= Halley's method code, with f''=0)

For naming consistency, Newton's = iter_i(), Halley's = iter_i2()

Code:

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

function iter_i(n, pv, pmt, fv, i)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n, pv, pmt, fv)
    return function()
        local eps, f = i*i
        if 1+eps==1 then
            local a, b = (pv+fv)/n, pv-fv
            local fp, fpp = (a+b)*0.5, (n*n-1)*a/6
            f, fp = (a+pmt)+fp*i, fp+fpp*i
            eps = -f / fp
        else
            local x = i/expm1(n*log1p(i))
            local k, y = (pv+fv)*x, n*x-1
            local num, den = y+(n-1)*i, i+eps
            f = k + pv*i + pmt
            eps = f / (k*num/den - pv)
        end
        i = i + eps -- Newton's method
        return i, eps, f
    end
end

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

lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> iter_i(n,pv,pmt,fv)() + 0 -- Newton's method
0.22222222032098768
lua> iter_i2(n,pv,pmt,fv)() + 0 -- Halley's method converge faster
0.3749999875677088