TVM solve for interest rate, revisited

[Albert Chan]

Previously (post #31), iter_i() returned i, f, eps = -f/fp
With verbose option, eps is not shown, because it can be deduced from i changes.

It may give better picture if eps is replaced by fp, and shown in verbose output.
With {i, f, fp}, we can mentally picture how the curve behave. (flatline? no solution?)

Another issue is find_rate() returned false for no solution case.
Boolean cannot be used for arithmetic calculations. nan is better.

Fudge factor to test converged i doubled. (i==i+eps/1000) --> (i==i+eps/2000)

Update:
iter_i() produce correct i, f, fp even if s = EFF(i,n) not finite.
find_rate(): no Newton step if f=0. (it may not be equivalent, if f'=0 too)

Here is the updated version.

Code:

function EFF(i,n,div) i=div and i/n or i; return expm1(log1p(i)*n) end
function APR(i,n,mul) i=expm1(log1p(i)/n); return mul and i*n or i end

-- NFV = pv*K + pmt*(K-1)/i + fv            , where K = (1+i)^n
function NFV(n,i,pv,pmt,fv) return (pv+pmt/i)*EFF(i,n) + (pv+fv) end
function NPV(n,i,pv,pmt,fv) return NFV(-n,i,fv,-pmt,pv) end -- time symmetry
-- NPMT = n*npmt = C*pv + (C-n*i)*fv + n*pmt, where C = i*n/(1-1/K)
function npmt(n,i,pv,pmt,fv) return ((pv+fv)/EFF(i,n)+pv)*i + pmt end
function NPMT(n,...) return n * npmt(n,...) end

function edge_i(n,i,pv,pmt,fv)  -- big edge if i is false
    pv, fv = pmt/-pv, pmt/fv    -- edges from inf, -1
    if (abs(pv)<abs(fv)) == (i==false) then pv,fv = fv,pv end
    return pv>-1 and finite(pv) or fv
end

function iter_i(n,i,pv,pmt,fv)
    i = i or edge_i(n,i,pv,pmt,fv)
    local f, fp = 0, 1
    return function()
        i = i - f/fp
        if 1+n*i*i==1 then
            local a = (pv+fv)/n         -- = f(0) - pmt
            local b = (n*n-1)*a/6*i     -- = f''(0)*i
            fp = (pv-fv+a)*.5 + b       -- = f'(0) + f''(0)*i
            f = a + pmt + (fp-.5*b)*i
        elseif i > -1 then
            local s = EFF(i,n)
            local k = (pv+fv)/s
            fp = k==0 and pv or k*(1-n*(s+1)/(s+s/i)) + pv
            f = (k + pv)*i + pmt
        else
            fp = -fv
            f = fp*i + pmt
        end
        return i, f, fp
    end
end

function find_rate(n,i0,pv,pmt,fv,verbose)
    local c, f0 = i0 and 2 or 1 -- just in case bad i0
    for i,f,fp in iter_i(n,i0,pv,pmt,fv) do
        if verbose then print(i,f,fp) end
        if f==0 then return i end
        if c>0 then c=c-1; f0=f; continue end
        if signbit(f*f0) then return i-.5*f/fp end
        if abs(f) < abs(f0) then f0=f; continue end
        c = -.5*f/fp
        return (i==i+c*.001) and i+c or nan
    end
end

function tvm(n,i,pv,pmt,fv, verbose)
    if not pv then n,pv,pmt,fv = -n,fv,-pmt end
    if not fv then
        if i==0 then return -n*pmt - pv end
        local s = signbit(i*n)  -- force z > 0
        local z = EFF(i, s and -n or n)
        return s and (pmt/i*z-pv)/(1+z) or (-pmt/i-pv)*z-pv
    end
    local flip = abs(pv) > abs(fv)
    if flip then n,pv,fv = (n and -n),-fv,-pv end
    if not pmt then
        if i==0 then return -(pv+fv)/n end
        return (-(pv+fv)/EFF(i,n)-pv)*i
    elseif not n then
        n = -(pv+fv)/(pmt+pv*i)
        if i ~= 0 then n = log1p(n*i)/log1p(i) end
        return flip and -n or n
    else    -- i used as guess, if supplied
        return find_rate(n,i,pv,pmt,fv, verbose)
    end
end

function tvm_begin(n,i,pv,pmt,fv, verbose)
    if not fv  then return tvm(n,i,pv+pmt,pmt,fv)+pmt end
    if not pv  then return tvm(n,i,pv,pmt,fv-pmt)-pmt end
    if not pmt then return tvm(n,i,pv,pmt,fv) / (1+i) end
    return tvm(n,i,pv+pmt,pmt,fv-pmt, verbose)
end

Example, slight changes to pmt, we have 0 solution (returned nan), changed to 2

lua> tvm(10,nil,1000,-287,2000,true)

Code:

-0.1435                 116.13771101577868      -1087.5631610761484
-0.036712904618008635   29.240213109966078      -536.0217925151969
0.017837513887852915    7.537006560375005       -262.9585107578723
0.04649985510150145     1.9377897954382775      -129.06227911365863
0.06151423272404261     0.5045985944708491      -62.25029703788141
0.0696201955624815      0.14253250628991054     -27.205477349325747
0.07485930621875686     0.05843692701773762     -4.949693275931963
0.08666547758333192     0.29079749268407795     44.04186343829406
nan

lua> tvm(10,nil,1000,-288,2000)
0.05503638712931415
lua> tvm(10,false,1000,-288,2000)
0.09744316384349011

Above is not typical. fp normally does not change that much, especially at the end.

lua> n,i,pv,pmt,fv = 48,nil,19198,-450,0 -- test sample #12
lua> tvm(n,i,pv,pmt,fv,true)

Code:

0.023439941660589644    220.49758413166887      13196.166742730682
0.006730726899179217    19.37233286981302       10823.538869300919
0.004940893361022924    0.24155367517272452     10553.2866320046
0.004918004408642897    3.969645661072718e-05   10549.817969654638
0.004918000645880665    1.0231815394945443e-12  10549.81739940963
0.004918000645880568    0                       10549.817399409601
0.004918000645880568