TVM solve for interest rate, revisited

[Albert Chan]

I made some fixes, and changes to the interface, here is the update
1. _find_rate() does not do time-symmetry anymore. tvm() already does.
2. _iter_i() updated with cleaner code
3. _iter_i2() removed. Halley's method may have trouble with tiny rate
4. edge_rate(n,i,pv,pmt,fv) can now return big edge, if i = false

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)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n,i,pv,pmt,fv)
    return function()
        local i0, f, eps = i
        if 1+n*i*i==1 then
            local a = (pv+fv)/n         -- = f(0) - pmt
            local b = (n*n-1)*a/12*i    -- = f''(0)/2*i
            local c = (pv-fv+a)*0.5 + b -- = f'(0) + f''(0)/2*i
            f = a + pmt + c*i
            eps = -f / (b+c)
        else
            local s = EFF(i,n)
            local k = (pv+fv)/s
            local d = k*(1-n*(s+1)/(s+s/i)) + pv
            f = (k + pv)*i + pmt
            eps = -f / d
        end
        i = i + eps
        return i0, f, eps
    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,eps in iter_i(n,i0,pv,pmt,fv) do
        if verbose then print(i,f) end
        if c>0 then c=c-1; f0=f; continue end
        if signbit(f*f0) or f==0 then return (i+eps/2) end
        if abs(f) < abs(f0) then f0=f; continue end
        return (i == i+eps*0.001) and (i+eps/2)
    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

(06-23-2022 07:04 PM)Albert Chan Wrote:  f value should read as if whole column get shifted up.

Now it is easier to read, verbose show (i, f(i)), without shifting column
It also show initial guess rate (previously hidden)

lua> tvm_begin(40, nil, 900, 1000, -1000, true)

Code:

-0.5                    -4.547473508864641e-11
-0.49999999999997724    0
-0.49999999999997724

Now, we can also supply rate guess (not recommeded for general use)
Code assumed user guess is bad, might overshoot, and use next one as initial guess.

lua> tvm_begin(40, 0, 900, 1000, -1000, true)

Code:

0                       997.5
-0.5118665811417575     -23.73316228353292
-0.4999999999999826     -1.0800249583553523e-11
-0.49999999999997724    0
-0.49999999999997724

Update 6/16/2024
Now, it can use big edge if i=false, small edge if i = nil

lua> tvm_begin(40, false, 900, 1000, -1000, true)

Code:

-0.5263157894736842     -52.63157894737378
-0.49999999999999134    -2.8194335754960775e-11
-0.49999999999997724    0
-0.49999999999997724

---

I am unsured when to use special branch (from f(0),f'(0),f''(0)) for iter_i().
Because of quadratic fit, it does not work well with huge compounding effect.
When rate is tiny, general branch (using f(i) formula) does not work.

For now, I think a compromise may be best: 1+n*i*i == 1

Update 6/25/24 patch is making more sense now, because below 2 are equivalent:

NPMT(n, i, pv, pmt, fv) == NPMT(n/2, i*(i+2), pv, pmt*(i+2), fv)

RHS rate = (1+i)^2-1 = i*(i+2)
--> number of payments cut in half. n → n/2
--> payments about doubled, pmt → pmt*(1+i) + pmt = pmt*(2+i)

LHS==RHS, we cannot base only on rate to decide when to switch.
We need to get rough invariant. If i is small, (n/2) * (i*(i+2)) ≈ (n*i)