PC-1211, PC-1250, etc. TVM

[Albert Chan]

To make the post more clearly, I quoted the code used here

(06-05-2024 01:05 AM)Albert Chan Wrote:  10 INPUT "B,N,P,M,F? ";B,N,P,M,F @ B=B*M
20 DEF FNF(I) @ S=-EXPM1(LOGP1(I)*-N)
25 D=(P+F)/S*(1-N*(1-S)/(S+S/I))+(B-F)
30 FNF=(P+F)*I/S+(B-F)*I+M @ END DEF ! FNF=NPMT
40 X=(ABS(P)<ABS(F))*(P+F)-B @ I=M/(X-P) @ Y=FNF(I)
50 H=-Y/D @ I=I+H @ Y=FNF(I) ! Newton's method
60 DISP 100*I,Y @ IF I+H*H<>I THEN 50

Since Sharp does not have EXPM1, LOGP1, lets try without them.
For good measure, we also report the next iteration.

>20 DEF FNF(I) @ S=1-(1+I)^(-N)
>70 DISP (I-Y/D)*100
>run
B,N,P,M,F? 0,360,0,-1000,1e6
.313028391431      504.47069099
.471270710728      62.96905505
.497374290313      1.51995026
.498036167469      .00095736
.498036584886      -.0000025
.498036583796

Final I is an improvement ... but it is really luck. Garbage in, Garbage out.
More Newton iterations suggested we only have 8 to 9 digits accuracy.

.498036583796
.498036583626
.498036583308
.49803658459
.498036583243
...

---

One way is to correct for S (1+I) error.

(1+I)^(-n) = (1+I+ε)^(-n) * (1 + ε/(1+I)) ^n ≈ (1+I+ε)^(-n) * (1 + n*ε/(1+I))

>20 DEF FNF(I) @ X=1+I @ S=X^(-N) @ S=1-S-N*(X-1-I)/X*S
>run
B,N,P,M,F? 0,360,0,-1000,1e6
.313028391431      504.47068755
.471270710304      62.96905734
.497374290775      1.51994785
.498036166884      .00095668
.498036584005      .00000001
.498036584009

---

Much better, but S calculation is messy. A novel way is to set ε=0 !!!

>20 DEF FNF(I) @ S=1-(1+I)^(-N)
>50 H=-Y/D @ I=I+H+1-1 @ Y=FNF(I) ! Newton's method
>run
B,N,P,M,F? 0,360,0,-1000,1e6
.313028391           504.47068892
.47127071             62.96905807
.497374291           1.51994733
.498036167           .00095642
.498036584           .00000002
.498036584009

[robve]

(06-05-2024 12:34 PM)Albert Chan Wrote:  Since Sharp does not have EXPM1, LOGP1, lets try without them.

Thanks, this is more clear.

EXPM1 LNP1 are easy to implement to compute TVM more accurately using a few terms of their series when applicable (up to 4th term of the series with 10 digits precision since 1+0.01^4/4 is close to 1 no need to use 5th term):

' compute L=LN(1+J) and S=(1+J)^-N and R=S-1
40 J=.01*I,K=1+J*B,C=0 : IF ABS J>=.01 LET L=LN(1+J) : GOTO 42
41 L=J*J,L=J-L/2+J*L/3-L*L/4
42 S=-N*L : IF ABS S>=.01 LET S=EXP S,R=S-1 : RETURN
43 R=S*S,R=S+R/2-S*R/6+R*R/24,S=R+1 : RETURN

With this, the secant method works quite well and is quick, small and simple. We just need to replace LN(1+J) with L and use the improved S and R:

Code:

10 "B" B=B=0 : PRINT "BGN=";MID$("NY",B+1,1) : END
11 "C" C=C=0 : BEEP C : END
12 "N" AREAD N : IF C GOSUB 40 : N=LN((M*K/J-F)/(M*K/J+P))/L
13 PRINT "N=";N : END
14 "J" AREAD I : IF C GOSUB 30
15 PRINT "I%=";I : END
16 "V" AREAD P : IF C GOSUB 40 : P=K*M*R/J-F*S
17 PRINT "PV=";P : END
18 "M" AREAD M : IF C GOSUB 40 : M=(F*S+P)*J/K/R
19 PRINT "PMT=";M : END
20 "F" AREAD F : IF C GOSUB 40 : F=(K*M*R/J-P)/S
21 PRINT "FV=";F : END
' solve Y=NPMT=0 for rate I with the secant method
30 G=9,I=1 : GOSUB 38
31 V=I,I=I/2 : GOSUB 38
32 IF ABS(Y-W)<I GOTO 31
33 T=I,I=I-(I-V)*Y/(Y-W),V=T : GOSUB 38
34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33
35 IF Y<>0 IF W<>0 IF Y<>W IF G LET G=G-1 : GOTO 33
' pick weighted average (interpolation):
36 IF SGN Y<>SGN W LET I=(I*ABS W+V*ABS Y)/(ABS Y+ABS W)
37 RETURN
38 W=Y : GOSUB 40 : Y=K*M-(P+F*S)*J/R : RETURN
' compute L=LN(1+J) and S=(1+J)^-N and R=S-1
40 J=.01*I,K=1+J*B,C=0 : IF ABS J>=.01 LET L=LN(1+J) : GOTO 42
41 L=J*J,L=J-L/2+J*L/3-L*L/4
42 S=-N*L : IF ABS S>=.01 LET S=EXP S,R=S-1 : RETURN
43 R=S*S,R=S+R/2+S*R/6+R*R/24,S=R+1 : RETURN

Line 32 is added to avoid the initial secant slope to be zero in extreme cases when the rate is very small, such as in examples 3,4,5 below.

The program always terminates with a result or gives an error (not likely, but could be in extreme and unrealistic TVM cases). When G counts down to zero from 9, the result is likely not correct because a root wasn't reliably found (not likely, but can happen in extreme TVM problems).

1) B=1,N=40,PV=900,FV=-1000,PMT=-40
comp I% => 4.753367189% is exact

2) B=1,N=40,PV=900,FV=-1000,PMT=-400
comp I% => 80% is exact in 10 digits, but 80.0000000009845%
comp FV => -1462.657745 (should be -1000, but TI BA gives -900(*) and Sharp EL-735 gives -975.11 and HP prime gives -1001.6 so it's tricky to compute FV accurately)

3) from http://www.voidware.com/tvm.htm
B=0,N=365*24*60*60,PV=0,PMT=-0.01,I%=10/N,
comp FV => 331667.0067 is exact
comp I% => 3.170979194E-07 is almost exact 3.170979198E-07 except for the last digit is 4 instead of 8

4) same as 3 but take N=365*24*60 and I%=10/N:
comp FV => 5527.782898 is exact
comp I% => 1.902587519E-05 is exact

5) same as 3 but take N=365*24*60 and I%=10/60/N:
comp FV => 5260.377543 is exact
comp I% => 3.170979051E-07 is very close to exact 3.170979198E-07

6) from http://www.voidware.com/tvm.htm (the hard example 2)
B=0,N=32,PV=-999999,FV=1000000,PMT=0
comp I% => 3.125027014E-06 is very close to exact 3.125E-06

7) B=0,N=360,PV=0,PMT=-1000,FV=1000000
comp I% => 4.980365842E-01 is exact (note that this is 5.976% APR)

EDIT: as a sanity check, I verified this program with 27 TVM problems with solutions posted by a university finance department's lectures.

To use this program on a Sharp pocket computer:
- Switch begin mode on/off: DEF-B
- Enter values: DEF-N or DEF-J (I%) or DEF-V (PV) or DEF-M (PMT) or DEF-F (FV)
- Calculate: press DEF-C (beep) then DEF-N or DEF-J or DEF-V or DEF-M or DEF-F

EDIT 2: once the N, PV, PMT, FV and I% are entered or computed and BGN is on/off, you can obtain amortization by adding five lines to the program:

Code:

50 "A" AREAD A : J=.01*I,R=P-F,T=0,V=0
51 FOR K=1 TO N : S=R*J,Q=-M-S,R=R-Q
52 IF K>=A LET T=T+S,V=V+Q : PRINT STR$ K;" ";INT(100*S+.5)/100;"int ";INT(100*Q+.5)/100;"prn"
53 NEXT K
54 END

Enter a period value then press DEF-A. When the interest paid and principle paid values are shown, press ENTER again to show the next period and so on.

The values T, V and R are the accumulated interest, accumulated principal and remaining balance. It should be easy to amend the program for direct P1/P2 and ACC computations if so desired.

EDIT 3: updated line 30 with I=1 to improve interest rate computation accuracy

- Rob

[dm319]

(06-05-2024 07:16 PM)robve Wrote:  but can happen in extreme TVM problems

I'm sure you're aware and everyone's bored to death of my requests and updates on this collation of tricky TVM problems and results. If you feel that one of these TVM problems can highlight a particular deficiency I might be tempted to add it in, and of course curious how your different algorithms and approaches fare.

[Albert Chan]

Hi, robve

If you have good EXP, LN, you can build accurate EXPM1 LNP1

(02-04-2019 07:28 PM)Dieters formula Wrote:  expm1(x) = (u-1) - (ln(u) - x) * u, where u = exp(x), rounded

log1p( x ) = ln(u) - ((u-1) - x) / u, where u = 1+x, rounded

And, many thanks to your examples.
I actually discovered a Plus42 bug, for rate very close to 0%.

Code:

10 INPUT "B,N,P,M,F? ";B,N,P,M,F @ B=B*M
15 DEF FNZ(I,A,B,C) @ D=(A+B)/2 @ B=C*A @ FNZ=A+M+(D+B)*I @ D=D+2*B @ END DEF
20 DEF FNF(I) @ IF 1+I=1 THEN FNF=FNZ(I,(P+F)/N,P-F+B+B,(N*N-1)/12*I) @ END
25 S=-EXPM1(LOGP1(I)*-N) @ D=(P+F)/S*(1-N*(1-S)/(S+S/I))+(B-F)
30 FNF=(P+F)*I/S+(B-F)*I+M @ END DEF
40 X=(ABS(P)<ABS(F))*(P+F)-B @ I=M/(X-P) ! edge rate
50 Y=FNF(I) @ DISP 100*I,Y @ H=-Y/D @ I=I+H @ IF I+H*H<>I THEN 50
60 DISP 100*I

With FNF(I ≈ 0) added, updated version handled all your examples.
Some F values adjusted to true 10 digits accuracy, to test for I round-trip

Code:

>run ! 1
B,N,P,M,F? 1,40,900,-40,-1000
 4.16666666667       -5.1783387869
 4.7524512828        -.0080716937
 4.75336718645       -.0000000187
 4.75336718857

>run ! 2
B,N,P,M,F? 1,40,900,-400,-1000
 66.6666666667       -66.666666756
 79.9999999644       -.000000183
 80.000000001

>run ! 3
B,N,P,M,F? 0,365*24*60*60,0,-0.01,331667.0067
-3.01507228575E-6     6.2977646195E-3
-1.1099216591E-7      7.022277322E-4
 3.07577857109E-7     1.52692722E-5
 3.17093002072E-7     7.8839E-9
 3.17097920053E-7
>res * n
 10.0000000068
 
>run ! 4
B,N,P,M,F? 0,365*24*60,0,-0.01,5527.7829015
-1.80904354932E-4     6.2977549796E-3
-6.659604136E-6       7.022284099E-4
 1.84546671658E-5     1.52693549E-5
 1.90255801086E-5     .000000007884
 1.90258751918E-5
>res * n
 10.0000000008

>run ! 5
B,N,P,M,F? 0,365*24*60,0,-0.01,5260.382426
-1.90100247286E-4     .005827435363
-2.255093023E-5       6.131837705E-4
-1.23294301E-7        1.15811689E-5
 3.16928063661E-7     4.4651E-9
 3.17097921621E-7
>res * 60*n
 10.0000000563
 
>run ! 6
B,N,P,M,F? 0,32,-999999,0,1000000
 0                    .03125
 3.12500161133E-6     0
 3.12500161133E-6
 
>run ! 7
B,N,P,M,F? 0,360,0,-1000,1000000
-.1                   2306.3389698
 .313028391458        504.47068745
 .47127071031         62.96905733
 .497374290777        1.51994785
 .498036166886        .00095669
 .498036584011

Update: Code showed initial edge rate as well (previously hidden)

[dm319]

(06-05-2024 07:16 PM)robve Wrote:  2) B=1,N=40,PV=900,FV=-1000,PMT=-400
comp I% => 80% (exact)
comp FV => -1462.657745 (should be -1000, but TI BA gives -900(*) and HP prime gives -1001.6 so it's tricky to compute FV accurately)

This is a really curious example!

If you have FV = -100,000 what do you get for I%?

If you start receiving 900, immediately paying out 400, then you have 500. 80% of 500 is 400 interest, which means the balance is 900 at the next payment period, when you are paying out 400. Ad infinitum presumably.

You're paying off your interest, so it's odd that these calculators can return -1000 for FV.

It's also odd they can return you an interest rate no matter how large FV is, yet the FV can never get more than -900.

It's also quite fun to solve for N with PV = 900, PMT = -400, FV = 900. On my DM42 solver I get a div/0 error, on plus42 it correctly tells me data is invalid. But then put in FV = 1000. Plus42 returns and error, but my solver gives me a negative number.

Confused.

[dm319]

(06-06-2024 12:27 AM)Albert Chan Wrote:  And, many thanks to your examples.
I actually discovered a Plus42 bug, for rate very close to 0%.

Me too! (not that surprising). Is this re: example 5 above? My code gets stuck with f and f0 dancing around zero. Maybe we should start a thread about this? In R it only seems to happen when precision is turn up beyond 250ish.

[Albert Chan]

(06-06-2024 11:02 PM)dm319 Wrote:  You're paying off your interest, so it's odd that these calculators can return -1000 for FV.

That's because robve other calculators do not have rate exactly 80%

lua> tvm_begin(40, nil, 900, -400, -1000)
0.8000000000098451
lua> tvm_begin(40, _, 900, -400, nil)
-1000.0000635560003

Quote:If you have FV = -100,000 what do you get for I%?

lua> tvm_begin(40, nil, 900, -400, -1e5)
0.8000000097564842
lua> tvm_begin(40, _, 900, -400, nil)
-100000.00038198779

[Albert Chan]

(06-06-2024 11:34 PM)dm319 Wrote:  

(06-06-2024 12:27 AM)Albert Chan Wrote:  And, many thanks to your examples.
I actually discovered a Plus42 bug, for rate very close to 0%.

Me too! (not that surprising). Is this re: example 5 above? My code gets stuck with f and f0 dancing around zero.

Yes, code had issue with example 5 (close to 0% bug).
Plus42 should have this fixed soon. (issue sent)

Plus42 previous version, secant's method did hang for some cases.
Now, with true Newton's method and edge guess, this issue goes away.

Here is closer to Plus42 do_i_pct_yr(), termination criteria based from Y, not I
It does not relied on quadratic convergence, but finite ULP's: ULP(FNF) ≥ ULP(M)
Rate search will not get stuck, even if no solution exist.

10 INPUT "B,N,P,M,F? ";B,N,P,M,F
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=A+M+(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) @ D=(P+F)/S*(1-N*(S+1)/(S+S/I))+P
50 FNF=((P+F)/S+P)*I+M @ 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

Code:

>run ! 1
B,N,P,M,F? 1,40,900,-40,-1000
 4.16666666667       -5.1783387869
 4.7524512828        -.0080716937
 4.75336718645       -.0000000187
 4.75336718857        0
 4.75336718857
>run ! 2
B,N,P,M,F? 1,40,900,-400,-1000
 66.6666666667       -66.666666756
 79.9999999644       -.000000183
 80.000000001         0
 80.000000001
>run ! 3
B,N,P,M,F? 0,365*24*60*60,0,-0.01,331667.0067
-3.01507228575E-6     6.2977646195E-3
-1.1099216589E-7      7.022277322E-4
 3.07577857249E-7     1.52692719E-5
 3.17093002025E-7     7.8839E-9
 3.17097920006E-7     0
 3.17097920006E-7
>res * n
 10.0000000053 
>run ! 4
B,N,P,M,F? 0,365*24*60,0,-0.01,5527.7829015
-1.80904354932E-4     6.2977549796E-3
-6.659604136E-6       7.022284099E-4
 1.84546671641E-5     1.52693549E-5
 1.90255801068E-5     7.8841E-9
 1.90258751938E-5    -6.E-14
 1.90258751927E-5
>res * n
 10.0000000013
>run ! 5
B,N,P,M,F? 0,365*24*60,0,-0.01,5260.382426
-1.90100247286E-4     .005827435363
-2.2550930229E-5      6.131837704E-4
-1.232943057E-7       1.15811690742E-5
 3.16928059687E-7     4.4652E-9
 3.17097921452E-7    -8.E-14
 3.1709791993E-7
>res * 60*n
 10.0000000029
>run ! 6
B,N,P,M,F? 0,32,-999999,0,1000000
 0                    .03125
 3.12500161133E-6     0
 3.12500161133E-6
>run ! 7
B,N,P,M,F? 0,360,0,-1000,1000000
-.1                   2306.3389698
 .313028391463        504.47068744
 .471270710313        62.96905732
 .497374290776        1.51994784
 .49803616688         .0009567
 .498036584009        0
 .498036584009

Update:
1. Small rate branch test (1+I=1) --> (1+N*I*I=1)
    Small rate branch does quadratic fit, not suitable for huge compounding effect.
2. time-symmetry to keep size of P small, possibly less cancellation error.
    With this change, "smallest" edge rate is simply M/F
3. Final Newton correction cut in half, to reduce error radius.

[robve]

(06-06-2024 12:27 AM)Albert Chan Wrote:  If you have good EXP, LN, you can build accurate EXPM1 LNP1

(02-04-2019 07:28 PM)Dieters formula Wrote:  expm1(x) = (u-1) - (ln(u) - x) * u, where u = exp(x), rounded

log1p( x ) = ln(u) - ((u-1) - x) / u, where u = 1+x, rounded

Like this cryptic line to compute L=LN(1+J) and S=(1+J)^-N and R=S-1:

Code:

40 J=.01*I,K=1+J*B,C=0,L=1+J,L=LN L-(L-1-J)/L,S=-N*L,R=EXP S,R=R-1-(LN R-S)*R,S=R+1 : RETURN

Replace line 40 to 43 with the above line. With Dieters formulas there are a few minor changes for the more extreme examples, while the other examples and the 27 TVM problems I tested remain unaffected:

2)
comp FV => -900 (which is closer to -1000 and the same result as TI BA, but is a tricky one)

4)
comp FV=5527.782902 (tiny bit better in the very last digit +/- 3)
comp I%=1.902587522E-05 (slightly better in the last 2 digits)

5)
comp FV=5260.382453 (slightly better in the last 2 digits)
comp I%=3.170988253E-07 (somewhat worse, but in both cases only first 5 digits are precise and this result just differs by more)

6)
comp I%=3.1248997E-06 (somewhat worse, but in both cases only the first 3 digits are precise and this result just differs by more)

It is a bit mixed pro/cons perhaps, but I do like the improved FV accuracy.

EDIT: I should mention that SHARP PC round to 10 digits when storing values in variables. But internally 12 digits are used. This matters when calculating EXPM1 and LNP1.

EDIT 2: do not use this method to compute EXPM1 and LNP1, because for larger negative S we get R=0 and LN(R) stops with an error.

- Rob

[dm319]

(06-07-2024 12:27 AM)Albert Chan Wrote:  lua> tvm_begin(40, nil, 900, -400, -1000)
0.8000000000098451

Ah! That makes sense. I was thrown by the 80% exact. I guess this problem is very affected by the number of digits of precision the calculator can store on that interest rate.

I added two of the problems to my bank of TVM problems. The HP-80 and HP-70 return -10.7 and -11.0% instead of 0.1666 recurring for problem 5 (Just realised I have measured an annual rate by i x 60x24x365). The early HPs all the way up to the HP12c seem to really struggle with this one. In fact the 12c regresses somewhat, with the 37E and 38E doing a bit better, and the HP-92 doing quite badly.

Weirdly for this problem, I thought I was getting issues with the solver (which is the plus42 solver basically), but testing it today it seems to have gone away and I don't quite understand why. It's returning a reasonable result down to 27 decimal points, with plus42 at 30 DP.

Problem 2 below is also interesting. The first HP that can do it basically only gets ~900, same for HP-12c. But we get much better results with the HP-20b/HP30b and HP-12c platinum which have more digits of precision. There was a surprise fail with the TI-83+, and it may relate to me taking the answer out, storing in a variable, and putting it back in (which is a PITA).

[robve]

(06-06-2024 12:27 AM)Albert Chan Wrote:  If you have good EXP, LN, you can build accurate EXPM1 LNP1

To follow up on my reply on Albert's suggestion to use this method, I found out that this causes errors in some cases such as when the interest rate is very small and N is large. The problem is that for larger negative S=-N*L we get R=0 and LN(R) quits with an error (or NaN):

Code:

40 J=.01*I,K=1+J*B,C=0,L=1+J,L=LN L-(L-1-J)/L,S=-N*L,R=EXP S,R=R-1-(LN R-S)*R,S=R+1 : RETURN

With the EXPM1 and LNP1 series in the previous code lines 40 to 43, there is no problem.

I made one other minor change to my code. With this change I am now able to accurately compute interest rates on the SHARPs with only 10 digits BCD and limited range (1E-99 to 9.999999999E+99), even for the more challenging cases:

1) unchanged, is exact
2) unchanged, is exact
3) I%=3.170979194E-07 is almost exact 3.170979198E-07
4) I%=1.902587519E-05 is exact
5) I%=3.170979051E-07 is very close to exact 3.170979198E-07
6) I%=3.125027014E-06 is very close to exact 3.125E-06, but HP-12C also gives an inexact result 3.125004736E-06
7) I%=4.980365842E-01 is exact

The change is to replace I=.02 with I=1 on line 30:

Code:

30 G=9,I=1 : GOSUB 38

The number of secant method iterations is about the same as before. The number of iterations for the examples are: 6, 5, 4, 14, 12, 10, 3, 7. With secant, the number of function evaluations is just the number of iterations, which is nice.

Note that I use a loop on line 32 that finds a suitable starting rate and this may take a few iterations in the worst case. But this is usually not necessary and only happens for cases 3, 4 and 5 that return very low interest rates with very large values of N which is a challenging combination.

I will update my previous post to include the updated code and improved results.

- Rob

[Albert Chan]

(06-08-2024 08:38 PM)robve Wrote:  To follow up on my reply on Albert's suggestion to use this method, I found out that this causes errors in some cases such as when the interest rate is very small and N is large. The problem is that for larger negative S=-N*L we get R=0 and LN(R) quits with an error (or NaN)

Thanks!

lua> LOG1P = fn'x,y: y=1+x; log(y) - (y-1-x)/y'
lua> EXPM1 = fn'x,y: y=exp(x); y-1 - y*(log(y)-x)'
lua> LOG1P(-1), EXPM1(-1e6)
-nan      -nan

Dieter's formula require testing for y==0 edge case.

lua> LOG1P = fn'x,y: y=1+x; log(y) - (y==0 and 0 or (y-1-x)/y)'
lua> EXPM1 = fn'x,y: y=exp(x); y-1 - (y==0 and 0 or y*(log(y)-x))'
lua> LOG1P(-1), EXPM1(-1e6)
-inf      -1

[robve]

(06-08-2024 06:19 PM)dm319 Wrote:  

(06-07-2024 12:27 AM)Albert Chan Wrote:  lua> tvm_begin(40, nil, 900, -400, -1000)
0.8000000000098451

Ah! That makes sense. I was thrown by the 80% exact. I guess this problem is very affected by the number of digits of precision the calculator can store on that interest rate.

I added two of the problems to my bank of TVM problems.

Indeed, 80% is exact in 10 places. I'm working with 10 digits only on the SHARPs.
I've added an annotation to my old post to explain.

If you want another bewildering TVM problem:
B=1 (begin mode)
N=40
PV=900
PMT=1000
FV=-1000
What's the interest rate?

In the real world this makes no sense, i.e. getting rich for nothing in return, but mathematically it is a valid TVM problem to solve. Try it on a calculator. Is the rate you get a whole number?

- Rob

[Albert Chan]

(06-08-2024 09:21 PM)robve Wrote:  If you want another bewildering TVM problem:
B=1 (begin mode)
N=40
PV=900
PMT=1000
FV=-1000
What's the interest rate?

lua> n, pv, pmt, fv = 40, 900, 1000, -1000
lua> pv, fv = pv+pmt, fv-pmt -- Now, b=0
lua> pv, pmt, fv
1900      1000      -2000

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)

lua> i = pmt / fv -- "smallest" edge rate = -0.5

With negative i and positive n, ((1+i)^n-1) --> -1, may lose many significant digits
So, we flip for time-reversed setup.

lua> n, pv, fv = -n, -fv, -pv
lua> x = -(pv+fv) / expm1(log1p(i)*n)
lua> x
-9.094947017737554e-11

i = pmt / (x - pv) = (-pmt/pv) / (1 - (x/pv)) ≈ (-pmt/pv) * (1 + (x/pv))

Since -pv was fv before time-reversal, first term is just old i = -0.5
It also mean old x = 0. No more iteration necessary.

lua> i .. ' + ' .. i*(x/pv)
-0.5 + 2.2737367544343885e-14

rate (%) = -49.9999999999977262632455656115

---

Equivalent way, tvm() to get x*i, divide by pv to get i error

lua> n,i,pv,pmt,fv
-40      -0.5      2000      1000      -1900
lua> tvm(n, i, 0, nil, pv+fv) / pv
2.2737367544343885e-14

[robve]

(06-08-2024 11:47 PM)Albert Chan Wrote:  

(06-08-2024 09:21 PM)robve Wrote:  If you want another bewildering TVM problem:
B=1 (begin mode)
N=40
PV=900
PMT=1000
FV=-1000
What's the interest rate?

rate (%) = -49.9999999999977262632455656115

Yes, excellent! Indeed -50% is not the exact answer beyond 12 digits. It can't be, given the TVM parameters. So is the rate you get (on a calculator) a whole number? Probably yes, but it fools you.

EDIT: With my modified TVM program when written in C using double precision:

Code:

./tvm 40 900 1000 -1000 1

TVM n=40 pv=900 pmt=1000 fv=-1000 BGN
i%=0.5    y=1007.24
1 i%=-51.3178964    y=-26.35792804 (9)
2 i%=-49.99647502    y=0.07049954146 (8)
3 i%=-50    y=-1.705302566e-13 (7)
4 i%=-50    y=1.136868377e-13 (6)
5 i%=-50    y=1.136868377e-13 (5)
i%=-49.9999999999977192%

The C program:

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double n, pv, pmt, fv, i, b = 0;

double fsign(double x) { return copysign(x != 0, x); }

double fn()
{
  double j = 0.01*i;
  double k = 1.0 + j*b;
  double l = log1p(j);
  double r = expm1(-n*l);
  double s = r + 1.0;
  return k*pmt - (pv + fv*s)*j/r;
}

void solve()
{
  int g = 9;
  int k = 0;
  double v, w, y;
  i = 1;
  y = fn();
  do {
    v = i;
    i /= 2;
    w = y;
    y = fn();
  } while (fabs(y - w) < i);
  printf("i%%=%g\ty=%g\n", i, y);
  do {
    do {
      double t = i;
      i -= (i - v)*y/(y - w);
      v = t;
      w = y;
      y = fn();
      printf("%d i%%=%.10g\ty=%.10g (%d)\n", ++k, i, y, g);
    } while (fsign(y) == fsign(w) && fabs(y) < fabs(w));
  } while (y != 0.0 && w != 0.0 && y != w && g-- > 0);
  if (fsign(y) != fsign(w))
    i = (i*fabs(w) + v*fabs(y))/(fabs(y) + fabs(w));
}

int main(int argc, char **argv)
{
  n = strtod(argv[1], NULL);
  pv = strtod(argv[2], NULL);
  pmt = strtod(argv[3], NULL);
  fv = strtod(argv[4], NULL);
  if (argc > 5)
    b = strtod(argv[5], NULL) != 0.0;
  printf("\nTVM n=%g pv=%g pmt=%g fv=%g", n, pv, pmt, fv);
  if (b)
    printf(" BGN");
  printf("\n");
  solve();
  printf("i%%=%.18g%%\n", i);
}

This uses a bunch of global variables for TVM parameters to make it easier to adapt to registers in a calculator.

- Rob

[dm319]

What a delicious problem! (can a TVM problem be delicious?) It's like the runner that goes half the distance to the finish line each second, but never quite gets there...

The earliest calculator I can try it out on is the HP-92, but basically it seems all HPs give the same answer, up to the HP-30b. RPN83P gives a different answer to the native TI-83+, but doesn't quite get another digit of accuracy, though it hints that the answer may not be exactly -50. I ran the problem in R using the rmpfr library to compare with plus42 and the like!

Code:

| calculator    | answer                                | accuracy |
|---------------|---------------------------------------|----------|
| HP-92         | -50                                   | -11      |
| HP-38E        | -50                                   | -11      |
| HP-12cp       | -50                                   | -11      |
| HP-30b        | -50                                   | -11      |
| TI-83 Plus    | -50                                   | -11      |
| RPN83P v0.9.0 | -49.999999999996                      | -11      |
| NSTK TVM V8   | -49.99999999999772626324556157894821  | -31      |
| Plus42 1.1.9  | -49.99999999999772626324556157894821  | -31      |
| C47 eadc4507f | -49.99999999999772626324556157894819  | -32      |
| R Rmpfr       | -49.999999999997726263245561578948194 |          |

[/code]

[rprosperi]

(06-08-2024 09:21 PM)robve Wrote:  [snip]

If you want another bewildering TVM problem:
B=1 (begin mode)
N=40
PV=900
PMT=1000
FV=-1000
What's the interest rate?

In the real world this makes no sense, i.e. getting rich for nothing in return, but mathematically it is a valid TVM problem to solve. Try it on a calculator. Is the rate you get a whole number?

- Rob

Using an HP-27S, I get -50 exactly.

[Albert Chan]

(06-02-2024 11:18 PM)Albert Chan Wrote:  

(06-02-2024 06:37 PM)Albert Chan Wrote:  rate = pmt/(x-pv)
f = x + (PV+FV) / ((1+rate)^N-1) = 0      // flatten curve to make Newton's method efficient
...
Interestingly, converge rate seems to match NPMT formula with Newton's method.

I just checked. Barring rounding errors, both ways are equivalent!

I found a more elegant proof, using reciprocal rate r = 1/i = (x-pv)/pmt

x2 = x - f / (df/dx)
(pmt*r2 + pv) = (pmt*r + pv) - f / (df/dr / pmt)
r2 = r - f / (df/dr)

Both are equivalent. We just need to show one has same convergence rate as npmt.

f(r) = (pmt*r + pv) + (pv+fv) / ((1+1/r)^n-1)

Compare this to g = npmt function (terms rearranged to make pattern clearer)

g(i) = (pmt + pv*i) + (pv+fv) / ((1+i)^n-1) * i

We have f = r*g. Let fp = df/dr, gp = dg/di

fp = d(r*g)/dr = (g + r * (dg/di) / (dr/di)) = (g + (1/i) * gp * (-i^2)) = (g - i*gp)

Newton iteration for next estimate:

r2 = r - f/fp = 1/i - (g/i) / (g - i*gp) = 1 / (i - g/gp) = 1 / i2      QED

f(r) = (pmt*r + pv) + (pv+fv) / ((1+1/r)^n-1)

Bonus trivia, from recipocal rate formula, f(1/i) = npmt(i)/i
n assumed positive. If not, we can use time-symmetry: {n,pv,fv} --> {-n,-fv,-pv}

f(r = -1) = limit(npmt(i)/i, i = -1) = (pmt*-1+pv) + (pv+fv)/(∞^-n-1) = (-fv-pmt)
f(r = 0) = limit(npmt(i)/i, i = +∞) = (pmt*0+pv) + (pv+fv)/(∞^n-1) = pv

For extreme rate, npmt(i) is basically linear, slope = (pv if n>0 else -fv)

[robve]

An improved version of the TVM program for SHARPs.

I kept it compact (22 lines of BASIC) and made it faster than the previous version. It takes only a few seconds to compute interest rates (faster than the SHARP EL-735). This version is also fully tested (on a SHARP PC-1403H).

I wanted to keep it small while returning accurate results, at least as accurate as possible with a 10 digit BASIC calculator that has limitations. This update uses LNP1 and EXPM1 formulas by Dieter and the initial interest guess by Albert. Also I%=0 is now permitted as input and is returned when the computed interest is zero or very small. This required a couple of tricks to avoid computation errors. There is no error trapping. When a computation fails then the result is ERR. Some explanations are included in the program.

How to use:
- Switch begin mode on/off: DEF-B
- Enter values: DEF-N or DEF-J (I%) or DEF-V (PV) or DEF-M (PMT) or DEF-F (FV)
- Calculate: press DEF-C (beep) then DEF-N or DEF-J or DEF-V or DEF-M or DEF-F

Code:

10 "B" B=B=0 : PRINT "BGN=";MID$("NY",B+1,1) : END
11 "C" C=C=0 : BEEP C : END
12 "N" AREAD N : IF C GOSUB 39 : IF N=0 LET N=LN((K*M-F*J)/(K*M+P*J))/L
13 PRINT "N=";N : END
14 "J" AREAD I : IF C GOSUB 29
15 PRINT "I%=";I : END
16 "V" AREAD P : IF C GOSUB 39 : P=K*M*R/J-F*S
17 PRINT "PV=";P : END
18 "M" AREAD M : IF C GOSUB 39 : M=(P+F*S)*J/K/R
19 PRINT "PMT=";M : END
20 "F" AREAD F : IF C GOSUB 39 : F=(K*M*R/J-P)/S
21 PRINT "FV=";F : END
' check if i%=0 (eval order of P+F+N*M is important)
29 C=0 : IF P+F+N*M=0 LET I=0 : RETURN
' solve Y=NPMT=0 for rate I% with the secant method, pick initial nonzero guess
30 G=3,Y=(P+F)/N,W=F-P-Y,V=(M+Y)/W,W=(N*N-1)*Y*V/W,I=100*V*(W-3)/(W-1.5)
31 I=I+(I=0)*1E-5 : GOSUB 37 : V=I,I=I/2 : GOSUB 37
' secant loop with count down to ensure termination
33 T=I,I=I-(I-V)*Y/(Y-W),V=T : GOSUB 37
34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33
35 IF Y<>0 IF Y<>W IF G LET G=G-1 : GOTO 33
' round very small I% to zero
36 I=(ABS I>1E-9)*I : RETURN
' compute Y=NPMT and save previous Y in W
37 W=Y : GOSUB 40 : Y=K*M-(P+F*S)*J/R : RETURN
' check edge case I%=0 to set J,K,S,R and compute N when requested
39 C=0 : IF I=0 LET J=1,K=1,S=1,R=-N,N=N-(N=0)*(P+F)/(M+(M=0)) : RETURN
' compute L=LN(1+J) and S=(1+J)^-N and R=S-1 using LNP1 and EXPM1 formulas
40 J=.01*I,K=1+J*B,L=1+J,L=LN L-(L-1-J)/L,S=-N*L,R=EXP S,R=R-1-(LN R-S)*R,S=R+1 : RETURN

Notes:
- line 30: a retry counter G is still necessary on this calculator, but G=3 suffices instead of G=9 with the improved initial guess
- line 31: we must avoid initializing I%=0 even when pv+fv+n*pmt is nonzero, so we assign I%=1E-5
- line 36: small |I%|<1E-9 are zero, based on experiments, this threshold may vary on other machines' MachEps
- line 39: if I%=0 then assign variables so that P, M, F and N are computed correctly with J/R=-1/N and K=R=S=1 and direct formula for N while avoiding div by zero when M=0

Add these four lines for amortization with period DEF-A after entering a TVM problem:

Code:

' amortization for given period and next, accumulated int and prn are T and V, remaining balance R
50 "A" AREAD A : J=.01*I,R=P-F,T=0,V=0
51 FOR K=1 TO N : S=R*J,Q=-M-S,R=R-Q
52 IF K>=A LET T=T+S,V=V+Q : PRINT STR$ K;" ";INT(100*S+.5)/100;"int ";INT(100*Q+.5)/100;"prn"
53 NEXT K : END

The results for the 12 TVM test examples look alright to me:

Code:

| #  | Ref        | N            | I%YR      | PV       | PMT          | FV        | P/YR | End |
|----|------------|--------------|-----------|----------|--------------|-----------|------|-----|
| 1  | DM         | 38 x 12      | 5.25%     | 270'000  | ?            | 0         | 12   | yes |
PMT=-1368.148355

| 1b | DM         | 38 x 12      | ?         | 270'000  | -14'584/12   | 0         | 12   | yes |
I%YR=4.373218369

| 2  | SlideRule  | 360          | 15% → 12% | 100'000  | ?-?          | 0         | 12   | yes |
-1264.444022 - -1028.612597 = -235.831425 => PV=6803.092173

| 3  | Kahan 1983 | 60x60x24x365 | 10%       | 0        | -0.01        | ?         | =N   | yes |
FV=331667.0067

| 4  | DM         | 480          | 0 → ?     | 100'000  | ?→ PMT       | 0         | 12   | yes |
I%=0

| 5  | Dieter     | 10           | ?         | 50       | -30          | 400       | 1    | yes |
I%=14.43587133

| 6  | Dieter     | 10           | ?         | 50       | -30          | 80        | 1    | yes |
I%=ERR

| 7  | A Chan     | 10           | ?         | -100     | 10           | 1e-10     | 12   | yes |
I%=0

| 8  | Miguel     | 32           | ?         | -999'999 | 0            | 1e6       | 1    | yes |
I%=3.12500161E-06

| 9  | DM         | ?            | 25        | 100000   | -2083.333334 | 0         | 12   | yes |
N=1040.639029 this one is not that good due to limited precision of this machine

| 10 | DM         | ?            | 25        | 100000   | -2040.816327 | 0         | 12   | no  |
N=1139.586234 also this is is not that good due to limited precision of this machine

| 11 | robve      | 60x24x365    | 1/6% → ?  | 0        | -0.01        | ?→ FV     | =N   | yes |
FV=5260.382453 => I%N=1.666673857E-01

| 12 | robve      | 40           | ? → I%YR  | 900      | -400         | -1000 → ? | 1    | no  |
I%=80 => FV=-900

Results for the C implementation of the same program with IEEE754 double precision, where I% is returned (to get I%YR multiply I% by 12 when P/YR=12 or multiply by N for tests 3 and 11). To compute the rate, the number of npmt function evaluations for root finding and retry counter G initialized to 2 are shown as evals/G. The retry counter G matters for TVM problems that dance around the root, such as n=40 pv=900 pmt=1000 fv=-1000 in BGN mode and n=360 pv=0 pmt=-1000 fv=1e+06 (not shown):

Code:

======== 1

TVM n=456 i%=0.4375% pv=270000 pmt=0 fv=0
pmt=-1368.1483553505775

======== 1b

TVM n=456 i%=0% pv=270000 pmt=-1215.33 fv=0
i%=0.182956    y=-341.776
 3/2  i%=    0.364314888635  y=   -0.244321551002
 4/2  i%=     0.36444462776  y=   0.0198834207638
 5/1  i%=     0.36443486391  y=-9.14124257179e-07
 6/0  i%=    0.364434864359  y=-3.41060513165e-12
 7/0  i%=    0.364434864359  y=                 0
i%=0.364434864359157562%

======== 2

TVM n=360 i%=1.25% pv=100000 pmt=0 fv=0
pmt=-1264.44402156504339

TVM n=360 i%=1% pv=100000 pmt=0 fv=0
pmt=-1028.61259692550448

TVM n=36 i%=1.25% pv=0 pmt=-235.831 fv=0
pv=6803.09216198562535

======== 3

TVM n=3.1536e+07 i%=3.17098e-07% pv=0 pmt=-0.01 fv=0
fv=331667.006690777082

======== 4

TVM n=480 i%=0% pv=100000 pmt=0 fv=0
pmt=-208.333333333333343

TVM n=480 i%=0% pv=100000 pmt=-208.333 fv=0
i%=0%

======== 5

TVM n=10 i%=0% pv=50 pmt=-30 fv=400
i%=7.17794    y=5.88612
 3/2  i%=     14.4161676678  y=   0.0120456394605
 4/2  i%=     14.4310106992  y=  0.00296966661476
 5/2  i%=     14.4358673527  y= 2.42831893971e-06
 6/2  i%=     14.4358713273  y= 4.90221196969e-10
 7/2  i%=     14.4358713281  y=                 0
i%=14.4358713280799602%

======== 6

TVM n=10 i%=0% pv=50 pmt=-30 fv=80
i%=-50.9714    y=10.8304
 3/2  i%=               nan  y=               nan
i%=nan%

======== 7

TVM n=10 i%=0% pv=-100 pmt=10 fv=1e-10
i%=9.09172e-12    y=5.00044e-12
 3/2  i%= 1.81802059632e-11  y=                 0
i%=1.81802059632393454e-11%

======== 8

TVM n=32 i%=0% pv=-999999 pmt=0 fv=1e+06
i%=1.5625e-06    y=0.015625
 3/2  i%= 3.12500161133e-06  y=                 0
i%=3.12500161132921604e-06%

======== 9

TVM n=0 i%=2.08333% pv=100000 pmt=-2083.33 fv=0
n=1060.30340666380812

======== 10

TVM n=0 i%=2.08333% pv=100000 pmt=-2040.82 fv=0 BGN
n=1076.31957397602309

======== 11

TVM n=525600 i%=3.17098e-07% pv=0 pmt=-0.01 fv=0
fv=5260.38242600032663

TVM n=525600 i%=0% pv=0 pmt=-0.01 fv=5260.38
i%=1.58549e-07    y=4.1684e-06
 3/2  i%= 3.17097919838e-07  y= 1.73472347598e-18
 4/2  i%= 3.17097919838e-07  y=                 0
i%=3.17097919837667697e-07%

======== 12

TVM n=40 i%=0% pv=900 pmt=-400 fv=-1000 BGN
i%=25.7127    y=-271.439
 3/2  i%=     79.9993961671  y= -0.00301916932335
 4/2  i%=     79.9999999949  y=-3.02858325085e-08
 5/2  i%=      80.000000001  y=-1.13686837722e-13
 6/2  i%=      80.000000001  y=                 0
i%=80.0000000009845138%

TVM n=40 i%=80% pv=900 pmt=-400 fv=0 BGN
fv=-999.999523662601518

The C program and a script to run the tests:

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double fsign(double x) { return copysign(x != 0, x); }

double n, pv, pmt, fv, i, b = 0; // TVM parameters
double j, k, l, r, s;            // dependent variables

double npmt()
{
  j = 0.01*i;
  k = 1.0 + j*b;
  l = log1p(j);
  r = expm1(-n*l);
  s = r + 1.0;
  return k*pmt - (pv + fv*s)*j/r;
}

void comp()
{
  if (i == 0.0) {
    j = k = s = 1.0;
    r = -n;
    if (n == 0.0)
      n = -(pv + fv)/pmt;
    return;
  }
  npmt();
}

double iguess()
{
  // https://www.hpmuseum.org/forum/thread-12941-post-160018.html#pid139867
  double a = (fv + pv)/n;
  double b = (fv - pv) - a;
  double c = (pmt + a)/b;
  b = (n*n-1)*a*c/b;
  return 100*c*(b - 3)/(b - 1.5);
}

void solve()
{
  int m = 1; // count number of evals
  int g = 2; // retries after overshooting
  double v, w, y;
  if (pv + fv + n*pmt == 0) {
    i = 0;
    return;
  }
  i = iguess();
  if (i == 0)
    i = 1e-5;
  y = npmt();
  v = i;
  i /= 2.0;
  w = y;
  y = npmt();
  ++m;
  printf("i%%=%g\ty=%g\n", i, y);
  do {
    do {
      double t = i;
      i -= (i - v)*y/(y - w);
      v = t;
      if (i != 0) { // to avoid I%=0 (not included in the BASIC version)
        w = y;
        y = npmt();
      }
      printf("%2d/%d  i%%=%18.12g  y=%18.12g\n", ++m, g, i, y);
      if (isnan(i))
        return;
    } while (fsign(y) == fsign(w) && fabs(y) < fabs(w));
  } while (y != 0.0 && y != w && g-- > 0);
  if (fabs(i) < 1e-16)
    i = 0;
}

int main(int argc, char **argv)
{
  n = strtod(argv[1], NULL);
  i = strtod(argv[2], NULL);
  pv = strtod(argv[3], NULL);
  pmt = strtod(argv[4], NULL);
  fv = strtod(argv[5], NULL);
  if (argc > 6)
    b = strtod(argv[6], NULL) != 0.0;
  printf("\nTVM n=%g i%%=%g%% pv=%g pmt=%g fv=%g", n, i, pv, pmt, fv);
  if (b)
    printf(" BGN");
  printf("\n");
  if (*argv[2] == '?') {
    solve();
    printf("i%%=%.18g%%\n", i);
  } else {
    comp();
    if (*argv[1] == '?')
      printf("n=%.18g\n", (n != 0.0 ? n : log((k*pmt-fv*j)/(k*pmt+pv*j))/l));
    else if (*argv[3] == '?')
      printf("pv=%.18g\n", k*pmt*r/j-fv*s);
    else if (*argv[4] == '?')
      printf("pmt=%.18g\n", (fv*s+pv)*j/k/r);
    else if (*argv[5] == '?')
      printf("fv=%.18g\n", (k*pmt*r/j-pv)/s);
  }
}

Code:

./tvm 40 ? 900 -40 -1000 1
./tvm 40 ? 900 1000 -1000 1
./tvm 40 ? 900 -400 -1000 1
./tvm 31536000 ? 0 -0.01 331667.0067
./tvm 525600 ? 0 -0.01 5282.36776893622391
./tvm 525600 ? 0 -0.01 5260.38242600032663
./tvm 32 ? -999999 0 1000000
./tvm 360 ? 0 -1000 1000000
echo
echo ======== 1
./tvm 456 0.4375 270000 ? 0
echo
echo ======== 1b
./tvm 456 ? 270000 -1215.3333333333 0
echo
echo ======== 2
./tvm 360 1.25 100000 ? 0
./tvm 360 1 100000 ? 0
./tvm 36 1.25 ? -235.83142463954 0
echo
echo ======== 3
./tvm 31536000 3.1709791983765e-07 0 -0.01 ?
echo
echo ======== 4
./tvm 480 0 100000 ? 0
./tvm 480 ? 100000 -208.333333333333343 0
echo
echo ======== 5
./tvm 10 ? 50 -30 400
echo
echo ======== 6
./tvm 10 ? 50 -30 80
echo
echo ======== 7
./tvm 10 ? -100 10 1e-10
echo
echo ======== 8
./tvm 32 ? -999999 0 1e6
echo
echo ======== 9
./tvm ? 2.08333333333333333 100000 -2083.333334 0
echo
echo ======== 10
./tvm ? 2.08333333333333333 100000 -2040.816327 0 1
echo
echo ======== 11
./tvm 525600 3.1709791983765e-07 0 -0.01 ?
./tvm 525600 ? 0 -0.01 5260.38242600032663
echo
echo ======== 12
./tvm 40 ? 900 -400 -1000 1
./tvm 40 80.0000000009844996 900 -400 ? 1

- Rob

[dm319]

(06-11-2024 07:13 PM)robve Wrote:  An improved version of the TVM program for SHARPs.

Great work!! Very impressive to fit all this in these lines of code. Also a very interesting performance profile. Eyeballing it, very close to the HP-37, sometimes identical. I will put it into the results table and I'm curious to run the dimension reduction on it again and see where it pops up. I'd be curious of a video of it in action solving tvm problems if you ever find the time.