TVM solve for interest rate, revisited
|
06-23-2022, 04:36 PM
Post: #21
|
|||
|
|||
RE: TVM solve for interest rate, revisited
(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) 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 |
|||
06-23-2022, 07:04 PM
(This post was last modified: 07-18-2022 02:20 PM by Albert Chan.)
Post: #22
|
|||
|
|||
RE: TVM solve for interest rate, revisited
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) 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 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) |
|||
06-23-2022, 10:21 PM
Post: #23
|
|||
|
|||
RE: TVM solve for interest rate, revisited
(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 |
|||
06-25-2022, 01:39 AM
(This post was last modified: 06-28-2022 01:18 PM by Albert Chan.)
Post: #24
|
|||
|
|||
RE: TVM solve for interest rate, revisited
(06-10-2022 08:25 PM)Albert Chan Wrote: Proof: f is either concave-up, or concave-down (f'' have same sign) 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) = + / + = + |
|||
06-28-2022, 12:55 PM
Post: #25
|
|||
|
|||
RE: TVM solve for interest rate, revisited
(06-25-2022 01:39 AM)Albert Chan Wrote: h is a polynomial, so does all its derivatives. 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' 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'' 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] |
|||
06-29-2022, 07:26 PM
(This post was last modified: 06-30-2022 12:16 PM by Albert Chan.)
Post: #26
|
|||
|
|||
RE: TVM solve for interest rate, revisited
This solved TVM calculations, 5 variables, 1 unknown
Code: function tvm(n,i,pv,pmt,fv) 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 |
|||
07-11-2022, 09:08 PM
Post: #27
|
|||
|
|||
RE: TVM solve for interest rate, revisited
(06-23-2022 07:04 PM)Albert Chan Wrote: NPV = NFV / K = NPMT / C 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 In this situation, we can use log version, from TVM formula to solve for n (*) (06-29-2022 07:26 PM)Albert Chan Wrote: If n ≠ -(pv+fv)/pmt, n has the form logab → 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. |
|||
07-13-2022, 03:34 PM
Post: #28
|
|||
|
|||
RE: TVM solve for interest rate, revisited
Code: D.__index = D 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 logab → n*ln(a) - ln(b) = 0 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 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 |
|||
07-16-2022, 02:58 PM
Post: #29
|
|||
|
|||
RE: TVM solve for interest rate, revisited
(07-11-2022 09:08 PM)Albert Chan Wrote: (*) Another way is simply get a better 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 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 |
|||
11-28-2023, 05:43 AM
Post: #30
|
|||
|
|||
RE: TVM solve for interest rate, revisited
tvm program assumed END mode. Here is tvm wrapper for BEGIN mode.
Code: function tvm_begin(n,i,pv,pmt,fv) > tvm_begin(n,i,360, .05/12, 100e3, nil, 0) -534.5941473979808 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)