TVM solve for interest rate, revisited

[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