TVM solve for interest rate, revisited

[Albert Chan]

(05-13-2022 10:34 PM)Albert Chan Wrote:  lua> n,pv,pmt,fv = 10,50,-30,100
lua> pmt/fv, -pmt/pv
-0.3      0.6

Lets try sign test, to see which way possible roots are located.

(06-14-2024 03:29 PM)Albert Chan Wrote:  sgn(-(pv+fv)/n) = sgn(pv*(i - pmt/-pv)) = sgn(-fv*(i - pmt/fv))

-(pv+fv)/n = -(50+100)/10 = -15

sgn(i - pmt/-pv) = sgn(-15/50) = -1      → i < pmt/-pv = 0.6
sgn(i - pmt/fv) = sgn(-15/-100) = 1      → i > pmt/fv = -0.3

True rates are 0.582038 and -0.284436 ... all checked out.

Lets try split loan "y" method

(06-13-2024 01:05 PM)Albert Chan Wrote:  

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

(pv+y) * (fv-y) < 0
(y+pv) * (y-fv) > 0
(y+50) * (y-100) > 0

pmt = -30, i = pmt/y:

(y<-50)      → i < -30/-50 = 0.6
(y>100)      → i > -30/100 = -0.3      -- same edges and direction, without memorized formula!

[Albert Chan]

Assume n>1, (x → ∞), we wanted to show (x - f/f') → -pmt/pv

Let K = (1+x)^n, g = x / (K-1)

(x → ∞) ⇒ (K → ∞) ⇒ (g → 0)

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

(x → ∞) ⇒ (f → pv*x + pmt) ⇒ (f' → pv)

x - f/f' → x - (pv*x + pmt)/pv = -pmt/pv

---

Assume n>1, (x → -1), we wanted to show (x - f/f') → pmt/fv

(x → -1) ⇒ (K → 0) ⇒ (g → 1)

ln(g) = ln(x) - ln(K-1)
g'/g = 1/x - n*K/(1+x)/(K-1) = 1/x - n/((1+x)*(1-1/K))

(x → -1) ⇒ (g'/g → 1/-1 - n/∞ = -1) ⇒ (g' → -1)

f = (pv+fv)*(1) + pv*(-1) + pmt = fv + pmt
f' = (pv+fv)*(-1) + pv = -fv

(x → -1) ⇒ (f → -fv*x + pmt) ⇒ (f' → -fv)

x - f/f' → x - (-fv*x + pmt)/(-fv) = pmt/fv

---

This allow for getting not just edge rate, but direction where it is going!

i > pmt/fv, because previous guess was -1
i < pmt/-pv, because previous guess was ∞

If we get smaller edge by picking bigger sized denominator, pmt=0 also work.

lua> tvm(10, nil, -10, 0, 20, true)
0                                      1
0.06896551724137931      0.03766397216483768
0.07176928097229544      5.600549027886197e-05
0.07177346252703834      1.2395358963512803e-10
0.07177346253629316      0
0.07177346253629316

Above example, we pick fv, because its size is bigger --> i > 0/20 = 0
Here is a nice mnemonic: < pv pmt fv >, same order as most financial calculator.

Of course, this is assuming root in that direction exist.
We might not have 2 solutions, even with 2 good edges.
Direction info can only show where is possible, or not.

Update: 6/18/24 I made HP71B Secants Method Rate Solver using asymptotic slopes.

[Albert Chan]

(06-10-2024 03:27 PM)Albert Chan Wrote:  Lets make bottom branch more accurate, with log1p_sub / expm1_sub
...
(log1p_sub / expm1_sub / fma) still does not help f with tiny rate!

I figured out why Newton on f does not work well with tiny rate!
Say, with non-zero a, we want to solve f(x) = g(x) - a

Newton's method: x ← x - f(x)/f'(x) = x - (g(x)-a) / g'(x)

With iteration zeroed-in to true root, g(x) ≈ a
f(x) will get less and less accurate, due to cancellation errors.
It does not matter if g(x) is super-accurate, even if correctly rounded!
Cancellation errors may happen for slope too, if a is also a function of x

This is normally not an issue. If f(x) only get half-precision, iterated x holds the other half.
With final f turns into garbage, more Newton's iteration may make it worse.

lua> x = 1
lua> for i=1,6 do f=x*x-2; print(x,f); x=x-f/(2*x) end
1                       -1
1.5                     0.25
1.4166666666666667      0.006944444444444642
1.4142156862745099      6.007304882871267e-06
1.4142135623746899      4.510614104447086e-12
1.4142135623730951      4.440892098500626e-16

But, that does not hold if X is very close to zero. Iterated x unable to keep its share.
This is the reason tiny rate top branch is needed, to keep cancellation errors at bay.
The other way is to throw in more-precision to the problem, but that is expensive.

[Albert Chan]

I recently translated find_rate() for HP71B (with test samples)
But, instead of copy/paste the code, I redo the math.
HP71B code seems simpler, with just as good accuracy.

I am updating old with new, in iter_i(), post #31

Code:

OLD                                             NEW

local x = i/EFF(i,n)                            local s = EFF(i,n)
local k, y = (pv+fv)*x, n*x-1                   local k = (pv+fv)/s
local num, den = y+(n-1)*i, i+i*i               local d = k*(1-n*(s+1)/(s+s/i)) + pv
f = k + pv*i + pmt                              f = (k + pv)*i + pmt
eps = f / (k*num/den - pv)                      eps = -f / d

Just to check math again, for df/di

s = (1+i)^n - 1
k = (pv+fv)/s       // note: (k*s) = (pv+fv) = constant

ln(i/s)' = ln(i)' - ln(s)'
(i/s)' / (i/s) = 1/i - n*(s+1)/(1+i)/s
s * (i/s)' = 1 - n*(s+1)/(s+s/i)

f = (k + pv)*i + pmt = (k*s)*(i/s) + pv*i + pmt
f' = (k*s)*(i/s)' + pv = k*(1-n*(s+1)/(s+s/i)) + pv

[Albert Chan]

Looking up old notes, there was a reason iter_i() was written that way.
Variables (num,den) were shared to build iter_i2(), Halley's method.

Code:

-- Halley's method iter_i2()

f = (pv+fv)/n*g + pv*x + pmt = 0
g = n*x/((1+x)^n-1)
g1/g = -(g-1+(n-1)*x)/(x+x^2) = -num/den
g2/g = (g+n*x)*(2*(g-1)+(n-1)*x)/(x+x^2)^2 = (num+x+1)*(num+(g-1))/den^2

But, if Newton's method required sharp f, Halley's even more so.
I removed iter_i2(), and only solve rate with Newton's from now on.

(01-03-2024 08:10 PM)dm319 Wrote:  And in part VI Plus42 is the only one that finds the negative root first. Is that because it is only going left?

Plus42 start guess from small sized edge, not left or right.

There is a reason for this. Think of mortgage with tiny ± final balance.
Tiny balance should have little effect to rate. We don't want the other weird rate.
It is the same when using quadratic formula. We keep the one that make sense.

Update: I just think of another good reason.
Small sized guess at worst goes down to 0; Big sized guess can get *very* big.
We may hit with overflows errors, instead of inconsequential underflow.

But, it is trivial to get either edges. Perhaps Plus42 should add this option?

Code:

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

lua> tvm(10,false,50,-30,100, true) -- big sized edge = -30/-50
0.6                                 0.8260581870403279
0.5821069833315318      0.0031464275173469503
0.5820382979602389      4.932464037210593e-08
0.5820382968834661      -3.552713678800501e-15
0.5820382968834661

Updates to post #31

[Albert Chan]

Previously (post #31), iter_i() returned i, f, eps = -f/fp
With verbose option, eps is not shown, because it can be deduced from i changes.

It may give better picture if eps is replaced by fp, and shown in verbose output.
With {i, f, fp}, we can mentally picture how the curve behave. (flatline? no solution?)

Another issue is find_rate() returned false for no solution case.
Boolean cannot be used for arithmetic calculations. nan is better.

Fudge factor to test converged i doubled. (i==i+eps/1000) --> (i==i+eps/2000)

Update:
iter_i() produce correct i, f, fp even if s = EFF(i,n) not finite.
find_rate(): no Newton step if f=0. (it may not be equivalent, if f'=0 too)

Here is the updated version.

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)
    i = i or edge_i(n,i,pv,pmt,fv)
    local f, fp = 0, 1
    return function()
        i = i - f/fp
        if 1+n*i*i==1 then
            local a = (pv+fv)/n         -- = f(0) - pmt
            local b = (n*n-1)*a/6*i     -- = f''(0)*i
            fp = (pv-fv+a)*.5 + b       -- = f'(0) + f''(0)*i
            f = a + pmt + (fp-.5*b)*i
        elseif i > -1 then
            local s = EFF(i,n)
            local k = (pv+fv)/s
            fp = k==0 and pv or k*(1-n*(s+1)/(s+s/i)) + pv
            f = (k + pv)*i + pmt
        else
            fp = -fv
            f = fp*i + pmt
        end
        return i, f, fp
    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,fp in iter_i(n,i0,pv,pmt,fv) do
        if verbose then print(i,f,fp) end
        if f==0 then return i end
        if c>0 then c=c-1; f0=f; continue end
        if signbit(f*f0) then return i-.5*f/fp end
        if abs(f) < abs(f0) then f0=f; continue end
        c = -.5*f/fp
        return (i==i+c*.001) and i+c or nan
    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

Example, slight changes to pmt, we have 0 solution (returned nan), changed to 2

lua> tvm(10,nil,1000,-287,2000,true)

Code:

-0.1435                 116.13771101577868      -1087.5631610761484
-0.036712904618008635   29.240213109966078      -536.0217925151969
0.017837513887852915    7.537006560375005       -262.9585107578723
0.04649985510150145     1.9377897954382775      -129.06227911365863
0.06151423272404261     0.5045985944708491      -62.25029703788141
0.0696201955624815      0.14253250628991054     -27.205477349325747
0.07485930621875686     0.05843692701773762     -4.949693275931963
0.08666547758333192     0.29079749268407795     44.04186343829406
nan

lua> tvm(10,nil,1000,-288,2000)
0.05503638712931415
lua> tvm(10,false,1000,-288,2000)
0.09744316384349011

Above is not typical. fp normally does not change that much, especially at the end.

lua> n,i,pv,pmt,fv = 48,nil,19198,-450,0 -- test sample #12
lua> tvm(n,i,pv,pmt,fv,true)

Code:

0.023439941660589644    220.49758413166887      13196.166742730682
0.006730726899179217    19.37233286981302       10823.538869300919
0.004940893361022924    0.24155367517272452     10553.2866320046
0.004918004408642897    3.969645661072718e-05   10549.817969654638
0.004918000645880665    1.0231815394945443e-12  10549.81739940963
0.004918000645880568    0                       10549.817399409601
0.004918000645880568

[Albert Chan]

Let S=(1+I)^N-1, D = N*I/S, C = subst(D,N=-N) --> C/D=(S+1), C-D=N*I

If solved rate is close to zero, this may be a better formula:

NPMT = C*PV + D*FV + N*PMT
NPMT = (C-1)*PV + (D-1)*FV + (PV+FV+N*PMT)

if I=0 --> (C-1) = (D-1) = 0, we are left with last term, F0 = (PV+FV+N*PMT)
In other words, if F0=0, then rate is 0. If solved rate is tiny, so does F0 too.

Formula is usable as is, but let's adjust it to match our iter_i():

NPMT = (D-1)*(PV+FV) + PV*N*I + F0

(D-1) = N*I/S - 1 = (N*I-S) / S

We get this iter_i(). Note: f = NPMT, not npmt, but trivial to scale (f, fp) by n, if needed.

Code:

function iter_i(n,i,pv,pmt,fv)          -- f = NPMT
    i = i or edge_i(n,i,pv,pmt,fv)
    local f, fp, f0 = 0, 1, fma(n,pmt,pv+fv)
    return function()
        i = i - f/fp
        if 1+n*i*i==1 then
            local a = (pv+fv)
            local b = (n*n-1)*a/6*i     -- = f''(0)*i
            fp = (n*(pv-fv)+a)*.5 + b   -- = f'(0) + f''(0)*i
            f = f0 + (fp-.5*b)*i
        else
            local s = EFF(i,n)
            local k = (pv+fv)/s
            fp = (k*(1-n*(s+1)/(s+s/i)) + pv)*n
            f = k*(n*i-s) + pv*n*i + f0
        end
        return i, f, fp
    end
end

(06-10-2024 01:33 AM)Albert Chan Wrote:  Problem 4:

lua> n, pv, fv = 480, 1e5, 0
lua> pmt = -pv/n -- assumed 0% interest
lua> pmt
-208.33333333333334
...
lua> i
1.8908413758272938e-19

lua> find_rate(480,nil,1e5,-1e5/480, 0, true)

Code:

0.0020833333333333333   58293.549096938164      31781021.329671178
0.00024910816818690847  6110.15308783893        25005937.09558987
4.760073393190968e-06   114.52326874360672      24068278.534060135
1.8075046879930214e-09  0.043470494011605396    24050006.94928785
2.6176695717432634e-16  6.290947846533815e-09   24050000.000001006
1.89084137588195e-19    1.3144820818133929e-22  24050000
1.8908413758272938e-19  0                       24050000
1.8908413758272938e-19

Of course, it work with normal cases too. Say, previous post examples.

lua> tvm(10,nil,1000,-287,2000)
nan
lua> tvm(10,nil,1000,-288,2000)
0.05503638712931402
lua> tvm(10,false,1000,-288,2000)
0.0974431638434901

But, I like simple formula, f = (k + pv)*i + pmt ... I am keeping old code.

[Albert Chan]

(06-15-2024 03:29 PM)Albert Chan Wrote:  I recently translated find_rate() for HP71B (with test samples)

Replaced with NPMT = (D-1)*(PV+FV) + PV*N*I + F0 --> FNF = NPMT/n, we have:

10 INPUT "B,N,P,M,F? ";B,N,P,M,F @ F0=(P+F+N*M)/N
12 IF B THEN P=P+M @ F=F-M ! end mode
14 IF ABS(P)>ABS(F) THEN X=F @ F=-P @ P=-X @ N=-N ! time reversed
20 DEF FNZ(I,A,B,C) @ D=(A+B)/2 @ B=A*C @ FNZ=F0+(D+B)*I @ D=D+2*B @ END DEF
30 DEF FNF(I) @ IF 1+N*I*I=1 THEN FNF=FNZ(I,(P+F)/N,P-F,(N*N-1)/12*I) @ END
40 S=EXPM1(LOGP1(I)*N) @ K=(P+F)/S @ D=K*(1-N*(S+1)/(S+S/I))+P
50 FNF=K*(I-S/N)+P*I+F0 @ END DEF
60 I=M/F @ Y=FNF(I) ! "smallest" edge guess
70 DISP 100*I,Y @ H=-Y/D @ I=I+H @ Y0=Y @ Y=FNF(I)
80 IF SGN(Y0)=SGN(Y) AND ABS(Y0)>ABS(Y) THEN 70
90 DISP 100*I,Y @ J=I-Y/D/2 @ DISP 100*J ! half-correction

Running test samples shown similar result, except for Problem 5

Code:

B,N,P,M,F? 0,365*24*60,0,-0.01,5260.382426
-1.90100247286E-4     5.82743536295E-3
-2.255093023E-5       6.13183770435E-4
-1.232943034E-7       1.15811690244E-5
 3.16928060096E-7     4.46510114E-9
 3.17097918102E-7     3.945E-14
 3.17097919603E-7     1.586E-14
 3.17097920206E-7     3.172E-14
 3.17097920809E-7
>res * 60*n
 10.0000000306

It is interesting we have a one solution case, but FNF have no sign changes.
Final FNF value get worse, on the same side of root!

I expect this behavior from a 2 solutions to 0 solution borderline case, but not here!

The problem here is F0 is tiny (if true rate is 0, so does F0), much smaller than M
Without needed cancellation, FNF is less stable. Sometimes, cancellation is good!

If we do an extra iteration, indeed we have sign changes, i.e. guaranteed solution.

>FNF(J)
-1.546E-14
>(J-RES/D/2) * 100 ! half-correction
 3.17097920515E-7

[Albert Chan]

Code:

function npmt(n,i,pv,pmt,fv) return ((pv+fv)/EFF(i,n)+pv)*i + pmt end

Due to decay shape d = i/EFF(i,n), if we have 2 edges, I think zero slope (if exist) is on the left of edges center.
If my gut is correct, we may use center of 2 edges to reach big root. (*)

Proof welcome! Disproof with counter example also nice.

Center guess may overshoot to the right, but with non-zero slope, should not go crazy big.
If overshooted, we would at worst wasted 2 iteration. (next of next back to edge)
If not, we saved all iterations, from big edge to center of 2 edges.

If there is only 1 valid edge, it is improved with I ≈ 1/P - P/N^2 formula, where 1/P = edge.
Improved guess may not be in correct region for 1-sided convergence, but Newton would fix it.

Both cases, first Newton iteration will not be used to test termination.

10 INPUT "B,N,P,M,F? ";B,N,P,M,F
15 IF B THEN P=P+M @ F=F-M ! end mode
20 DEF FNZ(I,A,C) @ C=A*C @ D=(P-F+A)/2+C @ FNZ=A+M+D*I @ D=D+C @ END DEF
30 DEF FNF(I) @ IF 1+N*I*I=1 THEN FNF=FNZ(I,(P+F)/N,(N*N-1)/12*I) @ END
40 S=EXPM1(LOGP1(I)*N)
50 D=(P+F)/S*(1-N*(S+1)/(S+S/I))+P @ FNF=((P+F)/S+P)*I+M @ END DEF
51 I=-1 @ IF F THEN I=M/F
52 J=-1 @ IF P THEN J=M/-P
54 IF I<=-1 THEN I=J-1/(J*N*N) @ GOTO 60
55 IF J<=-1 THEN I=I-1/(I*N*N) @ GOTO 60
57 I=(I+J)/2 ! 2 good edges
60 DISP 100*I,FNF(I) @ I=I-RES/D @ Y=FNF(I)
70 DISP 100*I,Y @ H=-Y/D @ I=I+H @ Y0=Y @ Y=FNF(I)
80 IF SGN(Y0)=SGN(Y) AND ABS(Y0)>ABS(Y) THEN 70
90 DISP 100*I,Y @ IF Y THEN I=I-Y/D/2 @ DISP 100*I

Line 90 only correct i when FNF(i)≠0. This is not an optimization.
I want D=0 to give divide-by-0 warning, but not when FNF(i)=0 already.

>run
B,N,P,M,F? 0,10,50,-30,80

Code:

 11.25               -16.6938992149
 139.952615683        40.0050644805
 59.7832470725        .6147808653
 58.4651530528        .0014842067
 58.4619552859        .0000000092
 58.4619552661       -.0000000001
 58.4619552662

Center guess is on the inside, Newton pushed it out.
2 iterations less with big edge guess, but won't know until tried it.

Here is a case center guess help, by a lot. (no overshoot this time)

>run
B,N,P,M,F? 0,10,50,-12.8277219567,80

Code:

 4.81039573375        .0048896134
 4.45818094799        .0012163008
 4.28359009886        .0003008647
 4.19806387387        .0000724336
 4.15834013895        .0000156502
 4.14322447289        .0000022677
 4.14014044165        .0000000945
 4.14000024866        .0000000001
 4.14000009969        0

Same problem, but with big edge guess.

>i=m/-p @ run60

Code:

 25.6554439134        3.7843086116
 13.5327655001        .816004209
 8.60628229142        .1941800627
 6.31331024557        .0475787922
 5.20035064432        .0117855308
 4.65151470732        .0029310812
 4.37935459727        .0007284364
 4.24473730211        .000179131
 4.17956294095        .0000420926
 4.15066279831        .0000082864
 4.14127814623        .0000008742
 4.14002330217        .0000000156
 4.14000007907        0

Example with only 1 valid edge.

>run ! 5
B,N,P,M,F? 0,365*24*60,0,-0.01,5260.382426

Code:

 3.17141393E-7       -1.14278E-9
 3.1709792021E-7      0
 3.1709792021E-7      0

Same problem, but with edge guess.

>i=m/f @ run 60

Code:

-1.90100247286E-4     .005827435363
-2.2550930229E-5      6.131837704E-4
-1.232943057E-7       1.15811690742E-5
 3.16928059686E-7     4.4652E-9
 3.17097921451E-7    -5.E-14
 3.170979205E-7

(*) I was picturing f(x) = npmt(n,x,pv,pmt,fv) in my head
With 2 valid edges, we can make g(x) = f(edge1 + edge2 - x)

2 curves would not overlap, because of decay factor.
1 curve will have zero slope (if exist) on the left, the other on the right.
Center guess (where 2 curve intersect) would not hit with zero slope.

plot {(130/((1+x)^10-1)+50)*x-30, (130/((1+(.6-.375-x))^10-1)+50)*(.6-.375-x)-30}, x = -.375 .. 0.6

[Albert Chan]

Last example, if we set P=0.01, big edge = 100%, center about 50%
Both guesses will cause S=EXPM1(LOGP1(I)*N) to overflow.
This does not mean algorithm is wrong ... we just need a better FNF, for big rate.

FNF=((P+F)/S+P)*I+M --> S=inf --> FNF=P*I+M

< 40 S=EXPM1(LOGP1(I)*N)
> 40 S=EXPM1(LOGP1(I)*N) @ IF S>9E499 THEN D=P @ FNF=D*I+M @ END

>run
B,N,P,M,F? 0,365*24*60,0.01,-0.01,5260.382426

Code:

 49.9999049499       WRN L40:Overflow
-5.00000950501E-3
WRN L40:Overflow
 100                  0
WRN L40:Overflow
 100                  0

How do I remove WRN L40:Overflow?
Is there exception handling (like Python), for HP71B?

---

100% rate? Yes! Every period, 1 penny loan doubled to 2, paid back 1

With exact 100% rate, we still owe 1 penny at the end (fv = -0.01)
But we have credit balance instead --> rate i = 1 - ε, ε>0

Split loan

Code:

    PV      PMT     FV
#1  0.01+x  -0.01   -0.01-x         // i = 0.01/(0.01+x)
#2  -x      0       5260.392426+x   // x*(1+i)^n = 5260.392426+x --> x = 5260.392426/((1+i)^n-1)

#1: x=0 --> i = 1
#2: i=1 --> x = 5260.392426/(2^525600-1) ≈ 2^-525587.64

i = 0.01/(0.01+x) = 1/(1+100x) ≈ 1 - 100x ≈ 1 - 2^-525581

[Albert Chan]

(06-15-2024 03:29 PM)Albert Chan Wrote:  I am updating old with new, in iter_i(), post #31

OLD: s=i/EFF(i,n); k=(pv+fv)*s; fp = -k*((n*s-1)+(n-1)*i)/(i+i*i) + pv; f = k + pv*i + pmt
NEW: s=EFF(i,n)  ; k=(pv+fv)/s; fp =  k*(1-n*(s+1)/(s+s/i)) + pv;         f = (k+pv)*i + pmt

I tested both versions with car lease example, n,pv,pmt,fv = 36,30000,-550,-15000
Using Plus42 reference, solved i = 0.00580507281942013213837742298371094

With default setup (i=nil), old way error = 1 ULP, new = -2 ULP
Note: error = exact - approx --> positive ULP means under-estimated.

This is 1 example, 1 guess, and little differences. We can't tell which is better.

I then compare example against million random guesses = random() - 0.5 (same sets)
Also, tested whether Newton final correction make sense (default was half-correction)

Newton     none               half             full
ULP     old     new       old     new       old     new
-9      0       0         0       0         0       0
-8      0       0         0       0         0       0
-7      0       0         0       0         0       0
-6      0       0         0       0         0       0
-5      25645   31924     25645   31924     25645   31924
-4      56093   55614     56093   55614     165769  55614
-3      87030   86220     87030   86220     166564  168111
-2      119581  111874    119581  111874    164470  164696
-1      146559  142173    256235  142173    169406  169815
0       159774  155031    238969  235542    168496  166121
1       136685  132674    181565  185453    139650  136125
2       109676  107594    23186   136616    0       107594
3       79195   80511     8731    11133     0       0
4       44880   52779     2965    3451      0       0
5       22847   27642     0       0         0       0
6       8722    11090     0       0         0       0
7       2965    3451      0       0         0       0
8       339     1380      0       0         0       0
9       9       43        0       0         0       0 

Final Newton correction produce better rate, with error within ±5 ULP
Half correction produce better looking bell-shaped curve, centerd at 0 ULP
Full correction looks good too, but distribution flattened, skewed to one side.

Same problem can be solved in time-reversed way:

Newton     none               half             full
ULP     old     new       old     new       old     new
-9      0       0         0       0         0       0
-8      0       0         0       0         0       0
-7      0       0         0       0         0       0
-6      0       0         0       0         23745   27474
-5      899     0         0       0         0       0
-4      54374   60232     54374   60232     54374   60232
-3      87113   89638     87113   89638     164745  165200
-2      117888  112837    118787  112837    166925  168921
-1      144228  139168    167973  166642    144228  139168
0       157697  159256    235329  234442    166737  168343
1       164800  162206    213815  218283    167998  164772
2       110643  105890    110643  106266    111248  105890
3       77632   75186     9062    9094      0       0
4       49015   56077     2299    2566      0       0
5       23745   27474     605     0         0       0
6       9040    9087      0       0         0       0
7       2299    2566      0       0         0       0
8       605     376       0       0         0       0
9       22      7         0       0         0       0

No correction, time-reversed setup looks just as bad.
Half-correction, time-reversed or not, looks about the same.
Full correction still skewed, but now with a weird overshoot (error = -6 ULP).

Old and new way seems to be similarly accurate.
Time-reversed to limit size of pv, and final half-correction seems best == tvm() default

[Albert Chan]

(06-21-2024 03:11 PM)Albert Chan Wrote:  Last example, if we set P=0.01, big edge = 100%, center about 50%
Both guesses will cause S=EXPM1(LOGP1(I)*N) to overflow.
This does not mean algorithm is wrong ... we just need a better FNF, for big rate.

Fixed iter_i() code to handle big rate overflow issue.
If s overflowed, it use asymptote line: f = pv*i + pmt

lua> n, pv, pmt, fv = 365*24*60, 0.01, -0.01, 5260.382426
lua> tvm(n,0.5,pv,pmt,fv, true) -- guess about edges center
0.5         -0.005     0.01
1            0            0.01
1
lua> tvm(n,false,pv,pmt,fv, true) -- big edge, beyond OK too
1            0            0.01
1

Update: guess can now be -100%, or less!
Implementation is cheap, if i ≤ -1, then f = -fv*i + pmt

This made guess formulas safer to use, as long as overshoot is still finite.
Of course, edge rates are safest, with guaranteed no overshoot.

Updates in post #46

[Albert Chan]

Issues with previous update.

1. if iterated rate is huge, Newton step might not get back edge rate.
What if i = ±inf ? No amount of correction will make it back to finite!
It is more accurate and direct to set next iteration as edge instead.

2. We want code to work with any bad guess, even if rate overshooted to ±inf

3. Returning f, fp based on asymptote line might not make sense (i<-1, f and fp do not exist!)

Instead, we just use f,fp formula, whatever its results stay.
Turns out, for bad rate, fp = nan, which is a good signal to use!

Code:

local s = EFF(i,n)
local k = (pv+fv)/s
fp = k*(1-n*(s+1)/(s+s/i)) + pv

If s=inf, then (s+1)/(s+s/i) = inf/inf = nan
If i=-1, then s=-1 or inf, both cases we have fp = nan
If i<-1, then s=nan, so does fp
If i is 1 ULP above -1, (s+s/i) ≠ 0 --> i ≤ -1 ⇒ fp = nan

4. For min points to use = c+1 (c calls to get f0, and we least at least 1 f's)

a) if guess from edge, first call is good for f0 --> c=1
b) supplied guess itself caused fp = nan. Next iteration would get back edge --> c=2
c) supplied guess next iteration caused fp = nan. Next would get back edge --> c=3

This is what I come up with: (c = number of calls to get good f0)

Code:

local c, f0 = i0 and 2 or 1 -- just in case bad i0
...
if c>0 then c = (c==1 and fp~=fp and .5) or c-1; f0=f; continue end

Here is the new update, which should work with *any* rate guess!

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, sortkey)
    fv, pv = pmt/fv, pmt/-pv    -- edges from -1, inf
    sortkey = sortkey or abs    -- big size edge if i is false
    if (sortkey(fv)>sortkey(pv)) ~= (i==false) then  fv,pv = pv,fv end
    return fv>-1 and finite(fv) or pv
end

function iter_i(n,i,pv,pmt,fv)
    local i0, f, fp = i or edge_i(n,i,pv,pmt,fv)
    return function()
        i = i0 or i - f/fp
        if 1+n*i*i==1 then
            local a = (pv+fv)/n         -- = f(0) - pmt
            local b = (n*n-1)*a/6*i     -- = f''(0)*i
            fp = (pv-fv+a)*.5 + b       -- = f'(0) + f''(0)*i
            f,i0 = a + pmt + (fp-.5*b)*i
        else
            local s = EFF(i,n)
            local k = (pv+fv)/s
            f  = (k+pv)*i + pmt
            fp = k*(1-n*(s+1)/(s+s/i)) + pv + (f-f) -- edge rate if bad slope
            i0 = fp~=fp and edge_i(n,i<=-1,pv,pmt,fv, function() return 0 end)
        end
        return i, f, fp
    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,fp in iter_i(n,i0,pv,pmt,fv) do
        if verbose then print(i,f,fp) end
        if f==0 then return i end
        if c>0 then c = (c==1 and fp~=fp and .5) or c-1; f0=f; continue end
        if signbit(f*f0) then return i-.5*f/fp end
        if abs(f) < abs(f0) then f0=f; continue end
        c = -.5*f/fp
        return (i==i+c*1e-6) and i+c or nan
    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-13-2024 09:28 PM)Albert Chan Wrote:  lua> n = 365*24*60 -- = 525600
lua> i = 1/(600*n) -- = 3.1709791983764586e-09
lua> fv = tvm(n,i,0,-1,nil)
lua> fv -- unit = cents
526038.2426000326
lua> tvm(n,nil,0,-1,fv) * n -- APR
0.0016666666666665247
...
What if we supply a guess, say 0.1% ?

lua> tvm(n,1e-3,0,-1,fv, true)
0.001                                      -1
-5.137033888903107e+219     -nan
-nan                                       -nan
false

lua> tvm(n,1e-3,0,-1,fv, true)

Code:

0.001                   -1                      -1.946648633485123e-220
-5.137033888903107e+219 -nan                    -nan
-1.9010024728569764e-06 0.5827435362965103      -347804.1847830881
-2.255093022962348e-07  0.061318377038474425    -273405.44161406445
-1.2329430676915646e-09 0.0011581169113648038   -263075.4360702088
3.169280601171723e-09   4.4651455866073775e-07  -262872.5774358277
3.1709791981232053e-09  6.661338147750939e-14   -262872.4991628829
3.1709791983766108e-09  -2.220446049250313e-16  -262872.4991629354
3.1709791983761885e-09

lua> tvm(n,0,0,-1,fv, true)

Code:

0                       0.0008337949011274493   -263018.6208831188
3.1700983691872297e-09  2.3154578809858606e-07  -262872.53975251823
3.170979198308128e-09   1.7763568394002505e-14  -262872.49916287535
3.170979198375703e-09   2.220446049250313e-16   -262872.4991630108
3.1709791983765475e-09  -1.1102230246251565e-16 -262872.49916294066
3.170979198376336e-09

lua> i = guess_i(n,i,0,-1,fv)
lua> tvm(n,i,0,-1,fv, true)

Code:

3.170979198350206e-09   6.661338147750939e-15   -262872.49916295544
3.1709791983755466e-09  2.220446049250313e-16   -262872.4991630238
3.170979198376391e-09   -1.1102230246251565e-16 -262872.49916295364
3.17097919837618e-09

Update:
It may be hard to get accurate f(i), due to cancellation errors.
If above examples had no sign changes, code would consider them no solution case.
With no sign changes, we considered rate found if they have 9 digits convergence. (was 12)

[Albert Chan]

Found a nice rate formula! https://brownmath.com/bsci/loan.htm#LoanInterest
Any clue how q-to-rate method is derived? Why log2?

fv = 0
i0 = pmt/-pv
q = log2(1+1/n)                    → i ≈ ((1+i0)^(1/q)-1)^q-1

Or, let Q = 1/q ≈ n*log(2)      → (1+i)^Q ≈ (1+i0)^Q - 1

Example from previous post

lua> n, pv, pmt, fv = 525600, 0, -1, 526038.2426000326
lua> 1/(600*n)                  -- true rate
3.1709791983764586e-09

lua> n, pv, fv = -n, -fv, -pv -- time reversal to match formula
lua> i0 = pmt/-pv
lua> q = log2(1+1/n)
lua> EFF(APR(i0,q)-1, q)
3.1710141058643196e-09

lua> i0 - 1/(i0*n^2)           -- compare to I ≈ 1/P - P/N^2
3.1714139412113494e-09

[Albert Chan]

(06-26-2024 02:31 PM)Albert Chan Wrote:  Found a nice rate formula! https://brownmath.com/bsci/loan.htm#LoanInterest
Any clue how q-to-rate method is derived? Why log2?

fv = 0
i0 = pmt/-pv
q = log2(1+1/n)                    → i ≈ ((1+i0)^(1/q)-1)^q-1

I finally figured this out

C*pv + n*pmt = 0                 → C = n*(pmt/-pv) = n*i0 > 0

Let s = (1+i0)^(1/q), then i ≈ (s-1)^q-1, very tiny if n is huge --> s ≈ 2

s = ((1+i0)^(1/i0))^(i0/q) < e^(i0/q) = e^(C/(n*q))

e^ln(2) = 2, we have C/(n*q) > ln(2)      → (n*q) < C / ln(2)

s assumed i0>0. For negative i0, it is (n*q) > C / ln(2)

The rest is likely by experimentation, with typical user input.
https://brownmath.com/bsci/loan.htm#Sample2

lua> n, i, pv, pmt, fv = 360, 0.0065, 225e3, -1619.71, 0
lua> i0 = pmt/-pv
lua> i0
0.007198711111111112
lua> q = log1p(1/n) / log(2)     -- slightly better than q = 1/n / log(2)
lua> EFF(APR(i0,q)-1,q)
0.006464576864221504
lua> q = 1/(n+0.5) / log(2)      -- log1p(1/n) = 2*atanh(1/(2n+1)) ≈ 1/(n+0.5)
lua> EFF(APR(i0,q)-1,q)
0.006464578259952884

lua> C = n*i0
lua> i0 - 1/C / n                       -- compare to I ≈ 1/P - P/N^2
0.006126845708593925

Too much correction, lets try the real thing. -1/C factor was from W function.
We pick W branch so W(y) ≈ -1/C, not -C. If C<1, we want more negative W_-1
Since C<1 only occurs with negative n, unless we applied time-reversal, it is W₀

lua> y = -C*exp(-C)
lua> i0 + W(y, C<1 and -1 or 0) / n
0.0065070286651448175

[Albert Chan]

(06-26-2024 02:31 PM)Albert Chan Wrote:  let Q = 1/q ≈ n*log(2)      → (1+i)^Q ≈ (1+i0)^Q - 1

Previous post showed why ln(2) is there, because RHS is (s-1)
If RHS is (s-2), then ln(3) pops out ... prove is not complete.

Let's try again, assuming n is huge. (i is interest rate, not √(-1))

C = i*n / (1-(1+i)^-n) ≈ i*n / (1-e^(-i*n)) = 1 + (n*i)/2 + (n*i)²/12 - ...

s = (1+i0)^Q = (1+C/n)^Q ≈ (1 + 1/n + i/2 + n/12*i² - ...)^Q

(1+1/n)^Q ≈ ((1+1/n)^n) ^ ln(2) ≈ e ^ ln(2) = 2

We can pull out tiny 1/n term, as correction instead (note: n ≈ Q/ln(2))

RHS = s-1 ≈ 2 * (1 + i/2 + n/12*i² - ...)^Q - 1 ≈ 1 + Q*i + 0.490*(Q*i)² - ...

LHS = (1+i)^Q = 1 + Q*i + Q*(Q-1)/2*i² + ...  ≈ 1 + Q*i + 0.500*(Q*i)² + ...

QED

With Q=n*ln(2), RHS is still smaller. This explained Q = (n+0.5)*ln(2) is better.

[Albert Chan]

With Q-method proven, I was able to improve its accuracy.

\((1+i)^Q \,≈\, (1+i_0)^Q - 1\)


\(\displaystyle
i \,=\, i_0 - \frac{i_0}{(1+i)^n}
\,≈\, i_0 - \frac{i_0}{((1+i_0)^Q-1)^\frac{n}{Q}}
\)

This is suitable for C = n*i0 > 1, because denominator has size around C^2

(06-27-2024 09:21 PM)Albert Chan Wrote:  lua> n, i, pv, pmt, fv = 360, 0.0065, 225e3, -1619.71, 0
lua> i0 = pmt/-pv
lua> i0
0.007198711111111112
lua> q = log1p(1/n) / log(2)     -- slightly better than q = 1/n / log(2)
lua> EFF(APR(i0,q)-1,q)
0.006464576864221504

lua> i0 - i0/APR(i0,q)^(n*q)
0.006491096452188186

lua> Q = (n+.5)*log(2)        -- my preference is big Q
lua> i0 - i0/EFF(i0,Q)^(n/Q)
0.006491096805454773