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