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

[Albert Chan]

f = ((pv+fv)/(K-1) + pv)*x + pmt      , where K = (1+x)^n

To improve accuracy, we noticed ULP(f) ≥ ULP(pmt)
If pv=0, we have ULP(f) = ULP(pmt)

We like to make |pv| smaller, if possible, to increase f precision.
We can do this with time-symmetry (formula does not care time direction):

NPV = NFV / K  = NPMT / C

K = (1+x)^n
C = (x*n)/(1-1/K) = (x*n)/(K-1) * K = (x*n)/(K-1) + n*x

NPV(n,i,pv,pmt,fv) = NFV(-n,i,fv,-pmt,pv)
NFV(n,i,pv,pmt,fv) = NPV(-n,i,fv,-pmt,pv)
NPMT(n,i,pv,pmt,fv) = NPMT(-n,i,fv,-pmt,pv)

f = NPMT/n      → f(n,i,pv,pmt,fv) = -f(-n,i,fv,-pmt,pv) = f(-n,i,-fv,pmt,-pv)

Code:

function _find_rate(iter_i, found, n,pv,pmt,fv,i0, verbose)
    if abs(pv) > abs(fv) then n,pv,fv = -n,-fv,-pv end
    local c, f0 = 2
    for i,eps,f in iter_i(n,pv,pmt,fv,i0) do
        if verbose then print(i,f) end
        if c>0 then c=c-1; f0=f; continue end
        if found(f,f0) then return i-eps/2, true end
        if abs(f) < abs(f0) then f0=f; continue end
        return i-eps/2, (i == i-eps*0.001)
    end
end

For convenience, I made function fn, that act like Python's lambda.
Note: Halley's method for f=0, convergence may not be 1-sided.

Code:

function find_rate(...) return _find_rate(iter_i, fn'a,b: a*b<=0', ...) end
function find_rate2(...) return _find_rate(iter_i2, fn'a,b: a==0', ...) end

Update, Jul 11,2022:

Rate iterators actually returned (i+eps), eps, f
In other words, (i,f) is the point, (i+eps) extrapolated value.

Final f of only a few ULP's, (i+eps) may not be that good.
A reasonable guess is true rate between i and (i+eps).
Patched code returned (i+eps/2) = (i+eps) - eps/2
This is patched result:

lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> find_rate2(n,pv,pmt,fv,nil,true)
0.3749999875677088      999999995.5
0.7951945637769868      161944286.30682382
1.3009220024991999      22940767.56093494
1.9230635724160559      3128383.8915700433
2.694052015186147       422117.19292708224
3.6428250853054216      56675.26724831162
4.70010922758184        7473.075259553417
5.291135425226106       837.9935225675524
5.32155201095865        25.713350899574266
5.321554259156788       0.0018612216055089448
5.321554259156787       -1.2505552149377763e-012
5.321554259156789       2.3874235921539366e-012
5.3215542591567875      true

Above verbose output, f value should read as if whole column get shifted up.

(i1, f1) = (0.375, 162.e6)
(i2, f2) = (0.795, 22.9e6)
...

It showed a sign flip when i ends in 788, flipped again when i ends in 787
With f down to few ULP's, correction is worse, with final i ends in 789, outside bracketed rate.
Patched code assumed true i ends in middle of (787, 789)

Update, Jul 18,2022:

fuzz factor for rate found widened, with less false negatives.

< return i-eps/2, (i == i-eps*0.01)
> return i-eps/2, (i == i-eps*0.001)

[Albert Chan]

(06-23-2022 07:04 PM)Albert Chan Wrote:  if abs(pv) > abs(fv) then n,pv,fv = -n,-fv,-pv end

An example to show improvement by keeping |pv| small.

f = 0      ⇒ -pmt/x = (pv+fv)/expm1(log1p(x)*n) + pv

Plus42 TVM, P/YR=1, END mode. Solve for I, then use it to solve back PMT

N, PV, PMT, FV = 10, 1E10, 10, -100
I%/YR --> -84.34981140003005709231099197156087
PMT -----> 9.999999999999999999999999899268213

N, PV, PMT, FV = -10, 100, 10, -1E10
I%/YR --> -84.34981140003005709231099197367593      // 1 ULP error
PMT -----> 9.999999999999999999999999999999781

[Albert Chan]

(06-10-2022 08:25 PM)Albert Chan Wrote:  Proof: f is either concave-up, or concave-down (f'' have same sign)

g = n*x / ((1+x)^n-1) = 1 - (n-1)/2*x + (n²-1)/12*x² - (n²-1)/24*x³ + ...

Previous proof showed g'' > 0

In order to do that, it also have to show g above x=0 tangent line.
I had wrongly assumed second part true, for my first proof attempt.

This proof is more direct, by showing g'' has no roots.
As long as g'' has the same sign, f'' also have same sign throughout.

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

XCAS> g := x*n/((1+x)^n-1)
XCAS> g := g(x=x-1); // shifted. Now, domain is x > 0
XCAS> h := normal(n/g)

(x^n-1) / (x-1)      // = x^(n-1) + x^(n-2) + ... + x + 1 > 1, if n>1

h is a polynomial, so does all its derivatives.

g'' = n*(h^-1)''

(h^-1)'' = (-h^-2 * h')' = (-h^-2 * h'' + 2*h^-3 * h'^2)
(h^-1)'' * h^3 = 2*h'^2 - h*h''

RHS looks promising to be positive; it is Halley's correction denominator.
But, all we needed is to show RHS never touch 0.

XCAS> m := factor(simplify(2*h'^2 - h*h''))

\(\displaystyle \frac {n\;x^{n-2} \cdot \left[
(n\!-\!1)\; x^{n+1}
\;-\; (n\!+\!1)\; x^{n}
\;+\; (n\!+\!1)\; x
\;-\; (n\!-\!1)
\right]}{\left(x-1\right)^{3}}\)

Above expression touched up a bit, to show polynomial sign changes.
For n>1, numerator has 3 sign changes, thus 1 or 3 positive roots.
Denominator had factor (x-1)^3, cancelled all numerator roots.

For n>1, m has no roots ⇒ g'' has no roots ⇒ f'' has no roots. QED

---

Not needed for the proof, but we might as well get sign(g'')

XCAS> factor(limit(m, x=1))      // (n>1) ⇒ (m>0) ⇒ (g''>0)

n^2*(n+1)*(n-1)/6

Or, we just imagine x ≫ n > 1, sign(m) = + / + = +

[Albert Chan]

(06-25-2022 01:39 AM)Albert Chan Wrote:  h is a polynomial, so does all its derivatives.

g'' = n*(h^-1)''

(h^-1)'' = (-h^-2 * h')' = (-h^-2 * h'' + 2*h^-3 * h'^2)
(h^-1)'' * h^3 = 2*h'^2 - h*h''

RHS looks promising to be positive; it is Halley's correction denominator.
But, all we needed is to show RHS never touch 0.

Direct way is to expand RHS polynomial, and make sure no negative coefficients.
Example, say n=5

h = 1 + x + x^2 + x^3 + x^4
h' = 1 + 2 x + 3 x^2 + 4 x^3
h'' = 2 + 6 x + 12 x^2

kth column = kth term coefs. (polynomial with increasing power).

Code:

    1   2   3   4       = h'
*   1   2   3   4       = h'
-----------------   
  1*1 1*2 1*3 1*4
      2*1 2*2 2*3 2*4
          3*1 3*2 3*3 3*4
              4*1 4*2 4*3 4*4
-----------------------------
    1   4  10  20  25  24  16 = h'^2
    2   8  20  40  50  48  32 = 2*h'^2

h' has (n-1) terms ⇒ (2*h'^2) has 2*(n-1)-1 = (2n-3) terms.
For k ≤ (n-1), (2*h'^2) k-th coeff = 2*sum(j*(k+1-j), j=1..k) = 2*comb(k+2,3)

Code:

  1*2 2*3 3*4           = h''
*   1   1   1   1   1   = h
---------------------
  1*2 2*3 3*4
      1*2 2*3 3*4
          1*2 2*3 3*4
              1*2 2*3 3*4
                  1*2 2*3 3*4
-----------------------------
    2   8  20  20  20  18  12 = h*h''

For k ≤ (n-2), (h*h'') k-th coeff = sum(j*(j+1), j=1..k) = 2*comb(k+2,3)

2 sides coefs matched, (n-2) terms, from the front.
Lets flip the order, and check differences from the back, (2n-3) - (n-2) = (n-1) terms.

(2*h'^2 - h*h''), k-th coeff, from the back
= sum(2*(n-j)*(n-(k+1-j)) - (n-j)*(n-(j+1)), j=1..k)
= sum(n^2 - n*(1+2*(k-j)) + (2*j*(k+1-j) - j*(j+1)), j=1..k)
= n^2*k - n*(k^2) + 0
= n*k*(n-k)

With k = 1 .. n-1, RHS coeffs all positive; other coeffs all zero.
With no sign changes, RHS have no positive roots.      QED

Just to confirm, with bigger n

XCAS> h := (x^n - 1)/(x - 1)
XCAS> n := 10
XCAS> e2r(2*h'^2 - h*h'')

[90, 160, 210, 240, 250, 240, 210, 160, 90, 0, 0, 0, 0, 0, 0, 0, 0]

XCAS> makelist(k -> n*k*(n-k), 1, n-1)

[90, 160, 210, 240, 250, 240, 210, 160, 90]

[Albert Chan]

This solved TVM calculations, 5 variables, 1 unknown

Code:

function tvm(n,i,pv,pmt,fv)
    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 = expm1(log1p(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)/expm1(log1p(i)*n)-pv)*i
    elseif not i then
        i, n = find_rate(n,pv,pmt,fv)
        return n and i  -- false = no solution
    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
    end
end

lua> n,pv,pmt,fv = 36, 30000, -550, -15000
lua> tvm(n,nil,pv,pmt,fv) -- i
0.005805072819420131
lua> i = _
lua> tvm(nil,i,pv,pmt,fv) -- n
36
lua> tvm(n,i,nil,pmt,fv) -- pv
30000
lua> tvm(n,i,pv,nil,fv) -- pmt
-550
lua> tvm(n,i,pv,pmt,nil) -- fv
-15000.000000000002

[Albert Chan]

(06-23-2022 07:04 PM)Albert Chan Wrote:  NPV = NFV / K  = NPMT / C

K = (1+x)^n
C = (x*n)/(1-1/K) = (x*n)/(K-1) * K = (x*n)/(K-1) + n*x

NPV(n,i,pv,pmt,fv) = NFV(-n,i,fv,-pmt,pv)
NFV(n,i,pv,pmt,fv) = NPV(-n,i,fv,-pmt,pv)
NPMT(n,i,pv,pmt,fv) = NPMT(-n,i,fv,-pmt,pv)

I just realized there is no financial function called NPMT (we do have NPV, NFV)
We can define NPMT as a TVM function that work the same way, if time direction reversed.

NPMT = C*pv + C(n=-n)*fv + n*pmt

NPMT tends to be more "straight" than (NPV, NFV), with values inside (NPV, NFV)
Thus, f = NPMT/n = 0 is great for solving rate x, using Newton's method.

But, with extreme compounding effect, f might still be too curvy.
Below is using Halley's method. Newton's method will be twice as long!

(06-23-2022 07:04 PM)Albert Chan Wrote:  lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> find_rate2(n,pv,pmt,fv,nil,true)
0.3749999875677088      999999995.5
0.7951945637769868      161944286.30682382
1.3009220024991999      22940767.56093494
1.9230635724160559      3128383.8915700433
2.694052015186147       422117.19292708224
3.6428250853054216      56675.26724831162
4.70010922758184        7473.075259553417
5.291135425226106       837.9935225675524
5.32155201095865        25.713350899574266
5.321554259156788       0.0018612216055089448
5.321554259156787       -1.2505552149377763e-012
5.321554259156789       2.3874235921539366e-012
5.3215542591567875      true

In this situation, we can use log version, from TVM formula to solve for n (*)

(06-29-2022 07:26 PM)Albert Chan Wrote:  

Code:

    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
    end

If n ≠ -(pv+fv)/pmt, n has the form log_ab      → n*ln(a) - ln(b) = 0

We have to use a guess better than edge rate, otherwise ln(b) → NaN
Let's use first iteration of Halley's method (from small edge) as rate guess

lua> D = require'dual'.D
lua> function L(x) return n*D.log1p(x) - D.log1p(-(pv+fv)/(pmt/x+pv)) end

lua> n,pv,pmt,fv = D.new(10), D.new(-100), D.new(10), D.new(1e10)
lua> x = D.new(0.375,1)
lua> for i=1,7 do print(D.newton(x,L(x))) end
2.2611317877571713     -15.546298358404194
4.413843278673193      -6.645171069338584
5.25277183822195       -1.5540048931122605
5.321177493828256      -0.10965326382959262
5.321554247893293      -0.0005973748729495298
5.3215542591567875     -1.7858138079418495e-008
5.3215542591567875     0

Note: outputs are (extrapolated x, L(x)); Shift up 2nd column to get points (x, L(x))

(*) Another way is simply get a better guess.
For above example, we may ignore pmt, and directly solve for rate.

pv*(1+x)^n + fv = 0      → x = expm1(log(-fv/pv)/n)

We have x = 5.310, which is a great starting guess.

[Albert Chan]

Code:

D.__index = D
function D.add(f,a) return D.new(f[1]+a,f[2]) end
function D.sub(f,a) return D.new(f[1]-a,f[2]) end
function D.mul(f,a) return D.new(f[1]*a,f[2]*a) end
function D.div(f,a) return D.new(f[1]/a,f[2]/a) end
function D.rcp(f,a) a = (a or 1)/f[1]; return D.new(a,f[2]*a/-f[1]) end
function D.newton(f,x) f=f(x); x=x:add(-f[1]/f[2]); return x,x[1],f[1] end

I updated above code to dual number, for object oriented, functional style.
Newton's method returned extrapolated x, without updated it back to x.

"D.__index = D" allowed we write D.log1p(x) as x:log1p()
With formula not tied to D, it can be used for other types.
add/sub/mul/div/rcp code are to allow constants also not tied to type.

---

(07-11-2022 09:08 PM)Albert Chan Wrote:  If n ≠ -(pv+fv)/pmt, n has the form log_ab      → n*ln(a) - ln(b) = 0

We have to use a guess better than edge rate, otherwise ln(b) → NaN
Let's use first iteration of Halley's method (from small edge) as rate guess

lua> D = require'dual'.D
lua> function L(x) return n*D.log1p(x) - D.log1p(-(pv+fv)/(pmt/x+pv)) end

L(x) = n*log1p(x) - ln(b)

b = 1 - (pv+fv)/(pmt/x+pv) = (pmt - fv*x) / (pmt + pv*x)

L(0) = n*log1p(0) - ln(pmt/pmt) = 0 - 0 = 0
L(x) introduced extra root at x=0, thus required pv+fv+n*pmt ≠ 0 check.
Bad guess might leads to this, with next iteration turns NaN, wasting time and effort.

Also, if rate guess is close to edge, ln(b) blows up, turning L(x) very curvy.
For Newton's method, better formula is to place ln(b) as denominator:

M(x) := log1p(x) / log1p(-(pv+fv)/(pmt/x+pv)) - 1/n

This also removed L "false" root of 0; M roots now matched f roots.

lua> D = require'dual'.D -- with above updated code
lua> function L(x) return x:log1p():mul(n) - x:rcp(pmt):add(pv):rcp(-(pv+fv)):log1p() end
lua> function M(x) return (x:log1p() / x:rcp(pmt):add(pv):rcp(-(pv+fv)):log1p()):sub(1/n) end
lua> n,pv,pmt,fv = 10, -100, 10, 1e10 -- constants are just numbers!
lua> a, b = pmt/fv, pmt/-pv -- rate edges
lua> a, b
1e-009      0.1

lua> require'fun'()
lua> eps = 1e-6
lua> x = D.new(b+eps, 1)
lua> zip(iter(D.newton,L,x), iter(D.newton,M,x)) :take(8) :each(print)

Code:

0.1000299805316535      0.10091091761367948
0.10079690287260726     0.5091108558877786
0.1185800651649927      2.6664564468120817
0.4711351213239042      4.663249410730011
2.4792994528821524      5.28599641783414
4.552404426854714       5.321453996394143
5.272567034486379       5.3215542583611155
5.321363352132326       5.3215542591567875

Except for the hole [a,b], M shape seems continuous.
Newton's method for M=0, we can even use the wrong edge!
Actually, this is even better guess than from the correct edge!

lua> x = D.new(a-eps, 1)
lua> zip(iter(D.newton,L,x), iter(D.newton,M,x)) :take(8) :each(print)

Code:

5.90875527827794e-006   0.80756269309747
-NaN                    3.0682528156090147
NaN                     4.8610143159046855
NaN                     5.3043434970095396
NaN                     5.321530792686004
NaN                     5.321554259113202
NaN                     5.321554259156787
NaN                     5.321554259156788

[Albert Chan]

(07-11-2022 09:08 PM)Albert Chan Wrote:  (*) Another way is simply get a better guess.
For above example, we may ignore pmt, and directly solve for rate.

pv*(1+x)^n + fv = 0      → x = expm1(log(-fv/pv)/n)

We have x = 5.310, which is a great starting guess.

Had we put back pmt term, we get:

newx = function(x) return expm1(log1p(-(pv+fv)/(pmt/x+pv))/n) end

newx(x = ±Inf) reduced to quoted pmt=0 approximation.
x ← newx(x) work well for huge x; it is almost as good as Newton's method.

x - (newx - x) / (newx - x)' ≈ x - (newx - x) / (0 - 1) = newx

lua> require'fun'()
lua> genx = function(f,x) f=f(x); return f, f, f-x end
lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> iter(genx,newx,Inf) :take(8) :each(print)

Code:

5.309573444801933       -Inf
5.321581577924096       0.012008133122162867
5.321554197006317       -2.738091777931828e-005
5.321554259298183       6.229186588768698e-008
5.321554259156466       -1.4171686046893228e-010
5.321554259156789       3.232969447708456e-013
5.3215542591567875      -1.7763568394002505e-015
5.3215542591567875      0

Some references, for how Lua iterator work

https://luafun.github.io/under_the_hood.html#iterators
http://lua-users.org/wiki/IteratorsTutorial
https://www.lua.org/pil/7.1.html

[Albert Chan]

tvm program assumed END mode. Here is tvm wrapper for BEGIN mode.

Code:

function tvm_begin(n,i,pv,pmt,fv)
    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)
end

> tvm_begin(n,i,360, .05/12, 100e3, nil, 0)
-534.5941473979808

[Albert Chan]

I made some fixes, and changes to the interface, here is the update
Note1: _find_rate() does not do time-symmetry anymore. tvm() already does.
Note2: _iter_i() updated with cleaner code

Code:

local expm1, log1p = mathx.expm1, mathx.log1p

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,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

function iter_i2(n, i, pv, pmt, fv)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n, pv, pmt, fv)
    return function()
        local i0, f, eps = i
        if 1+n*i*i==1 then
            local a, b = (pv+fv)/n, pv-fv
            local fp, fpp = (a+b)*0.5, (n*n-1)*a/6
            f = a+pmt + i*(fp + i/2*fpp)
            fp, fpp = fp+fpp*i, fpp*(1-1.5*i)
            eps = -f*fp / (fp*fp - .5*f*fpp)
        else
            local x = i/EFF(i,n)
            local k, y = (pv+fv)*x, n*x-1
            local num, den = y+(n-1)*i, i+i*i
            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
        return i0, f, eps
    end
end

function iter_i(n, i, pv, pmt, fv)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n, 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(iter_i, found, 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 found(f,f0) 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 find_rate(...) return _find_rate(iter_i, fn'a,b: a*b<=0', ...) end
function find_rate2(...) return _find_rate(iter_i2, fn'a,b: a==0', ...) 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 n*i < -0.5 then n = pmt/i; n=log((n-fv)/(n+pv))/log1p(i)
        elseif 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

We can also supply rate guess.
Code assumed supplied 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

---

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

[dm319]

I did* translate your code into R for use with the (R)mpfr library (in a very clunky manual way - it was only to serve the basic purpose of getting more precise results than the DM42/Plus42 can output).

But just for fun, here is it with the problem above, initial guess of 10%, solving with 1000 bits of precision:

Code:

  iter         i              f             f0            eps
1    1 -0.524179   1.189774e+03   0.000000e+00  -6.241790e-01
2    2 -0.500000  -4.835808e+01   1.189774e+03   2.417904e-02
3    3 -0.500000  -2.586899e-11  -4.835808e+01   1.293450e-14
4    4 -0.500000  -2.251961e-35  -2.586899e-11   1.125981e-38
5    5 -0.500000  -1.706569e-83  -2.251961e-35   8.532843e-87
6    6 -0.500000 -9.800526e-180  -1.706569e-83  4.900263e-183
7    7 -0.500000   0.000000e+00 -9.800526e-180   0.000000e+00
8    8 -0.500000   0.000000e+00   0.000000e+00   0.000000e+00

With i as:

Code:

1 'mpfr' number of precision  1000   bits
[1] -0.4999999999999772626324556157894819480015601928133859453506649561424082849
50469649920734546388757663388950218950670857339905568445420132396647739217005786​
23316123318913635334886714726022977608726558724596098528753080203217196000004117​
8985411578382337606374842897420492915931795226261455970009197588717915

* probably before your changes above, though I don't think you've changed the algorithm of the newton method?

[Albert Chan]

(06-09-2024 04:46 PM)dm319 Wrote:  But just for fun, here is it with the problem above, initial guess of 10%, solving with 1000 bits of precision ...
* probably before your changes above, though I don't think you've changed the algorithm of the newton method?

It is the same Newton's method, f(i) = npmt(n,i,pv,pmt,fv)
But now {i, f(i)} is on the same row, for easy reading.

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

Code:

0.1                     1189.774058558563
-0.5241790424643634     -48.35808492873343
-0.4999999999999902     -2.580691216280684e-11
-0.4999999999999773     -1.1368683772161603e-13
-0.49999999999997724    0
-0.49999999999997724

Last rate is really i - (f/f')/2, to reduce error radius.
Since f=0, there is no correction, last = previous

[dm319]

Yes that code does look cleaner. Looks like you're getting similar results, the differences in the later f is probably due to the differences in precision.

Not sure if this will work, but hopefully a pretty graph:

   

It was quite hard to get the scales working on this graph...

[dm319]

Another bit of fun, this is for problem 4, which switches to the secant method as it hits zero.

Code:

   iter             i             f            f0            eps   type
1     1  1.532628e-03  6.408126e+02  0.000000e+00  -6.800705e-03 newton
2     2  1.469909e-04  8.609612e+01  6.408126e+02  -1.385637e-03 newton
3     3  1.684695e-06  7.451270e+00  8.609612e+01  -1.453062e-04 newton
4     4  2.265215e-10  8.442160e-02  7.451270e+00  -1.684469e-06 newton
5     5  4.096406e-18  1.134967e-05  8.442160e-02  -2.265215e-10 newton
6     6  1.339647e-33  2.052470e-13  1.134967e-05  -4.096406e-18 newton
7     7  1.432731e-64  6.712188e-29  2.052470e-13  -1.339647e-33 newton
8     8 1.638754e-126  7.178582e-60  6.712188e-29  -1.432731e-64 newton
9     9 2.143937e-250 8.210842e-122  7.178582e-60 -1.638754e-126 newton
10   10  0.000000e+00 1.074202e-245 8.210842e-122 -2.143937e-250 secant
11   11  0.000000e+00  0.000000e+00 1.074202e-245   0.000000e+00 secant
12   12  0.000000e+00  0.000000e+00  0.000000e+00   0.000000e+00 secant
13   13  0.000000e+00  0.000000e+00  0.000000e+00   0.000000e+00 secant
14   14  0.000000e+00  0.000000e+00  0.000000e+00   0.000000e+00 secant

I'm not quite sure why it seems to do an extra few rounds in the secant loop. It might be that the high precision can't be represented using the doubles here. I converted to double so that dealing with dataframes and graphs was familiar.

[Albert Chan]

(06-09-2024 06:42 PM)dm319 Wrote:  Another bit of fun, this is for problem 4, which switches to the secant method as it hits zero.

When rate is really tiny, there is almost no compounding effect.
We can treat curve as linear, perhaps quadratic.
iter_i() top branch does exactly that, f(ε) ≈ f(0) + f'(0)*ε + f''(0)/2!*ε²

Problem 4:

lua> n, pv, fv = 480, 1e5, 0
lua> pmt = -pv/n -- assumed 0% interest
lua> pmt
-208.33333333333334

C*pv + n*pmt = 0
C = -n*pmt / pv
C-1 = (n*pmt+pv) / -pv = (n+1)/2*i + (n^2-1)/12*i² + ...

Even if we use full precision, in binary or decimal, pmt has some error. (C-1 ≠ 0)
Assume C-1 as linear is good enough. Solving Quadratic give same rate.

lua> C1 = fma(n,pmt,pv) / -pv
lua> C1
4.5474735088646414e-17
lua> i = C1 / ((n+1)/2)
lua> i
1.8908413758272938e-19

[dm319]

Thanks Albert. Yes, this was the reason I switched over to your solver for the TVM - it was failing at i = 0 to converge on 0 using the built-in solver. (PS do you know how the built-in solver works?).

I thought it was neat to see it switch over to the top branch halfway through solving this.

I don't quite understand why if there is little/no compounding effect the Newton solver doesn't quite home in properly.

What would be the equation to view the curve? I'm not quite understanding how if the curve is straight that would be a problem for Newton?

PS - also you mentioned in another thread that one of Rob's examples had pointed out a bug in the Plus42 code - I'm now realising you must have come across the example and bug a while ago, and this isn't a recent thing? I was saying that I experienced the bug in my code, but it turns out I had set it to 'end' rather than 'begin', and the code was running fine.

[Albert Chan]

(06-10-2024 08:57 AM)dm319 Wrote:  I don't quite understand why if there is little/no compounding effect the Newton solver doesn't quite home in properly.

Because f(ε) is hard to compute accurately. Even worse for f'(ε) ... Suggestions welcome!
This is why top branch is needed, for tiny rate.

Let's try Problem 4 with tvm, but tiny rate top branch disabled.

Code:

< if 1+n*i*i==1 then
> if false then

lua> tvm(n,nil,pv,pmt,fv, true)
0.0020833333333333333   121.44489395195453
0.0002491081681869074   12.7294855996644
4.760073393189775e-06    0.2385901432158164
1.8075046870790729e-09   9.056352914171839e-05
2.6088558403073705e-16   1.3073986337985843e-11
-4.488825652889797e-19   0
-4.488825652889797e-19

f(ε) is bad ... f'(ε) not shown, but likely worse.

Lets make bottom branch more accurate, with log1p_sub / expm1_sub
More specifically, we wanted accurate y = D-1, where D = n / USFV(i,n)
Let x = log1p_sub(i) = log1p(i) - i

y ← (1+i)^n - 1 - n*i
= expm1(n*log1p(i)) - n*i
= expm1(n*(x+i)) - n*i
= expm1_sub(n*(x+i)*n) + n*(x+i) - n*i
= expm1_sub((x+i)*n) + n*x

y ← -y / (y+n*i)      -- = D-1

old Wrote:local x = i/expm1(n*log1p(i))
local k, y = (pv+fv)*x, n*x-1
local num, den = y+(n-1)*i, i+i*i
f = k + pv*i + pmt
eps = f / (k*num/den - pv)

new Wrote:local x = log1p_sub(i)
local y = expm1_sub(n*(x+i)) + n*x
local k = y + n*i
y, k = -y/k, (pv+fv)/k
f = fma(k+pv,i,pmt)
eps = (f+f*i) / (fma(n-1,i,y)*k-pv*i-pv)

lua> tvm(n,nil,pv,pmt,fv, true)
0.0020833333333333333   121.44489395195444
0.0002491081681869104   12.72948559966454
4.7600733931911846e-06  0.23859014321589925
1.8075046869901344e-09  9.056352914562415e-05
2.6102777224213185e-16  1.3096744612389123e-11
-3.625568200933119e-19   -4.837592392112496e-14
1.201966826538812e-19

(log1p_sub / expm1_sub / fma) still does not help f with tiny rate!

I am not too concern with actual rate. Getting to 0% really depends on luck.
The issue is with bad Newton correction, it might mess up convergence.

Code don't require quadratic convergence, but must be one-sided!
Unanticipated overshoot might cause code to quit too early.

Quote:What would be the equation to view the curve?

{f(0), f'(0), f''(0)} evaluated symbolically using limit.

lua> f0 = fma(n,pmt,pv+fv) / n
lua> f1 = (pv-fv + (pv+fv)/n)/2
lua> f2 = (pv+fv)/n * (n*n-1)/6
lua> f0, f1, f2
-9.473903143468002e-15   50104.166666666664   7999965.277777779

f(i) ≈ p(i) = f0 + f1*i + f2/2*i²

lua> i = -f0/f1
lua> i
1.8908413758272935e-19
lua> i = -(f0 + f2/2*i*i)/f1
lua> i
1.8908413758272935e-19

Note: top branch does not actually use f(i). It only use quadratic p(i)

Quote:PS - also you mentioned in another thread that one of Rob's examples had pointed out a bug in the Plus42 code - I'm now realising you must have come across the example and bug a while ago, and this isn't a recent thing?

It was a recent thing.

Example 5 solved i ≈ 3.17e-9, 1+i*i==1 goes to top-branch. (Plus42 binary version)
But it also have huge compounding effect, quadratic fit simply does not cut it.

I then changed it to 1+i==1, but now tiny i numerical issue switched to general branch.
As a compromise, for now, I use 1+n*i*i==1 to switch to tiny rate branch

(1+i)^n = (1+APR/PYR) ^ (N*PYR)

n*i = N*APR, likely not crazy sized number

[Albert Chan]

I recently discovered split loan method. Is this new?

If loan split correctly, we can get a good guess rate, and formula to iterate for true rate.

(06-08-2024 11:47 PM)Albert Chan Wrote:  Solve with split loan method
f(x) = x + (pv+fv) / ((1+i)^n-1), where i = pmt / (x-pv)

Code:

Loan:  n  pv    pmt  fv
#1:    n  pv-x  pmt  x-pv      --> i = pmt / (x-pv)
#2:    n  x     0    (pv+fv)-x --> x = -(pv+fv) / ((1+i)^n-1)

Here is one with huge true rate (see post #22, #27, #28, #29)

lua> n, pv, pmt, fv = 10, -100, 10, 1e10
lua> tvm(n,nil,pv,pmt,fv, true)

Code:

1e-09                   999999995.5
0.22222222032098768     345130893.0191475
0.4241393063954123      127299482.5393054
0.6292750775369145      48106454.844047576
0.8468793905100045      18381292.76136444
1.0824224688331676      7062941.954833955
1.3400583139839632      2722233.2107647057
1.6234867740344179      1051054.2519718618
1.936322276082247       406212.4885288175
2.282247095736833       157054.32380092412
2.6649766778452406      60694.70062511578
3.0878290650500033      23400.014453658747
3.5521226148603975      8951.90999749287
4.052051396405174       3345.878170726684
4.560402776655654       1168.0978535777294
5.001365887839277       335.11729441390037
5.257160227141495       55.734261383772946
5.3187739972480506      2.306061198884663
5.321549001716665       0.004352499460765102
5.321554259137976       1.557395989948418e-08
5.3215542591567875      5.684341886080801e-13
5.321554259156788       -1.2505552149377763e-12
5.3215542591567875

tvm() picked time-reversed setup (x=0), edge = pmt/fv = 1e-9
The other edge (also x=0) is just as bad, edge = pmt/-pv = 0.1
Convergence is slow, because huge rate caused f(i) to be curvy.

(07-16-2022 02:58 PM)Albert Chan Wrote:  newx = function(x) return expm1(log1p(-(pv+fv)/(pmt/x+pv))/n) end

Above is from post#29, I just realize this can be easily explained as a split loan!
If we replace log1p(z) = ln(1+z), ln argument is:

-(pv+fv) / (pmt/x+pv) + 1 = -(fv-pmt/x) / (pv+pmt/x)

We might as well simplify this, with y = pmt/x

Code:

Loan:  n  pv        pmt     fv
#1     n  pv+y      0       fv-y        --> (1+i)^n = -(fv-y)/(pv+y)
#2     n  -y        pmt     y           --> i = pmt/y

We don't have to use the setup to solve rate. Guess i when y=0 is already good.

lua> y = 0
lua> i = EFF(-(pv+fv)/(y+pv), 1/n) -- guess
lua> tvm(n,i,pv,pmt,fv, true)
5.30957344513139        9.999999722758162
5.321456835153608      0.08065988603186725
5.321554252697288      5.347635692487529e-06
5.321554259156789      -1.3642420526593924e-12
5.321554259156788

Here is another way to split loan. (any more ways?)

Code:

Loan:  n  pv        pmt     fv
#1     n  pv+z      0       fv        --> (1+i)^n = -fv/(pv+z)
#2     n  -z        pmt     0         --> z = pmt * (1-(1+i)^-n)/i

Sign checks:
#1: sgn(pv+z) = sgn(-fv)
#2: sgn(z) = sgn(n*pmt), because NPMT is time-symmetrical, f=NPMT/n is not

lua> for loop=1,7 do i=APR(-fv/(pv+z)-1,n); z=pmt*-EFF(i,-n)/i; print(i,z) end
5.309573444801933      1.8833904463248339
5.321581577921809      1.8791405816968803
5.321554197006327      1.879150250410811
5.321554259298183      1.8791502284142771
5.321554259156466      1.8791502284643202
5.321554259156789      1.879150228464206
5.3215542591567875     1.8791502284642068

---

I don't know how to name this ... any idea?

EFF(i,n) = (1+i)^n - 1
APR(i,n) = (1+i)^(1/n) - 1

But google search effective rate formula, it was normally defined with mul/div.
To me this is less useful, so I added a third argument, in case true EFF, APR needed.

EFF(i,n,true) = EFF(i/n,n)
APR(i,n,true) = APR(i,n)*n

lua> EFF(0.12, 12, true)
0.12682503013196972
lua> APR(_, 12, true)
0.12

[Albert Chan]

Code:

Loan:  n  pv    pmt  fv
#1:    n  pv-x  pmt  x-pv      --> i = pmt / (x-pv)
#2:    n  x     0    (pv+fv)-x --> x = -(pv+fv) / ((1+i)^n-1)

We are not done with split loan yet ... Lets scale PV column to -1

Code:

#1:    n -1     i    1   , where i = pmt/(x-pv)
#2:    n -1     0    1+i2, where i2 = -(pv+fv)/x

(1+i)^n = 1+i2      → sgn(i2) = sgn(n*i)      (*)

i2 / i = -(pv+fv) / (pmt/i+pv) / i = -(pv+fv) / (pmt+pv*i)

sgn(-(pv+fv)/n) = sgn(pmt+pv*i) = sgn(pmt-fv*i)      // RHS from time-symmetry

Or, rewrite to show edge rate (note: ∞ is not an edge ... we might not have 2 edges)

sgn(-(pv+fv)/n) = sgn(pv*(i - pmt/-pv)) = sgn(-fv*(i - pmt/fv))

Based on signs, there is no reason to start from wrong side of edge rate.
It still work, but would converge slower than starting guess from edge.

OTTH, rate better than edge may overshoot ... and we don't know where!

---

(*) Let K_n = (1+x)^n > 0      // integer n, x = (-1, ∞)
Let S_n = K_n - 1, we want to show sgn(S_n) = sgn(n*x)

S_a+b = (S_a+1)*(S_b+1) - 1 = S_a*(1+S_b) + S_b = S_a*K_b + S_b

S1 = (1+x) - 1 = x
S2/S1 = (x*(1+x) + x)/x = (x+2) > 1
S3/S2 = (x*K2 + S2)/S2 = K2/(S2/S1) + 1 > 1
S4/S3 = (x*K3 + S3)/S3 = K3/(S3/S1) + 1 > 1
...

0 = S_n-n = S_-n K_n + S_n

(S_n≥1 / x ≥ 1) and (S_n / S_-n = -K_n < 0)      ⇒ sgn(S_n) = sgn(n*x)      QED
In other words, APY and APR have same sign, even if n is negative.

---

This perhaps is more elegant proof. Bonus: it covered real n too!

log1p(x)' = 1/(x+1) > 0
expm1(x)' = exp(x) > exp(-1) > 0

Both are increasing function. (we don't care the shape, only signs)

expm1(0) = log1p(0) = 0      --> sgn(expm1(x)) = sgn(log1p(x)) = sgn(x)

S_n = expm1(n*log1p(x))
log1p(S_n) = n*log1p(x)
sgn(S_n) = sgn(n*x)              QED