(71B) calculate interest rate

[Albert Chan]

Calculate interest rate, given N, FV, PV, PMT
For term (1+i)ⁿ - 1, it use the more accurate expm1(n * log1p(i)) = FNE(n*FNL(i))

Code:

10 INPUT "n, fv, pv, pmt ? ";N,F,P,M
20 DEF FNE(X) @ Y=EXP(X) @ FNE=Y-1-Y*(LN(Y)-X) @ END DEF
30 DEF FNL(X) @ Y=1+X @ FNL=LN(Y)-(Y-1-X)/Y @ END DEF
40 DEF FNF(I) @ Y=F+P+(M/I+P)*FNE(N*FNL(I)) @ DISP I,Y @ FNF=Y @ END DEF
50 I1=-(M*N+P+F)/(N*((N-1)/2*M+P))
60 I2=.9*I1
70 DISP FNROOT(I1,I2,FNF(FVAR))

Example, from https://www.hpmuseum.org/forum/thread-6725.html

>RUN
n, fv, pv, pmt ? 168, 10925.76, 0, -45
4.79863130882E-3     -655.6829061
5.33181256535E-3     -1255.8749515
4.21615468581E-3     -41.4514876
4.17684617474E-3     -1.4838737
4.17538677149E-3     -.0035171
4.17538330417E-3     -.0000003
4.17538330387E-3     0
4.17538330387E-3

Note both guesses (I1,I2) over-estimated interest rate on purpose.
If they bracketed the true rate, secant's method tends to over-shoot all over.

>60 I2=.5*I1
>RUN
n, fv, pv, pmt ? 168, 10925.76, 0, -45
2.66590628268E-3     1404.76181571
5.33181256535E-3     -1255.8749515
3.99885942402E-3     177.2179603
4.16369398215E-3     11.8480798
4.17461116346E-3     .7831243
4.7532118644E-3       -606.2697042
4.1753575833E-3       .0260875
4.46428472385E-3     -298.0104305
4.17538287345E-3     .0004364
4.31983379865E-3     -147.7508024
4.1753833001E-3       .0000038
4.24760854938E-3     -73.5647316
4.17538330383E-3     0
4.17538330383E-3

[Joe Horn]

(12-20-2019 04:14 PM)Albert Chan Wrote:  For term (1+i)ⁿ - 1, it use the more accurate expm1(n * log1p(i)) = FNE(n*FNL(i))

Wouldn't using the HP-71's built-in EXPM1 and LOGP1 functions be even more accurate?

[Albert Chan]

(12-20-2019 08:36 PM)Joe Horn Wrote:  Wouldn't using the HP-71's built-in EXPM1 and LOGP1 functions be even more accurate?

Thanks! I don't know HP71B had both expm1 and logp1 built-in. The revised code is much shorter.

Code:

10 INPUT "n, fv, pv, pmt ? ";N,F,P,M
20 DEF FNF(I) @ Y=F+P+(M/I+P)*EXPM1(N*LOGP1(I)) @ DISP I,Y @ FNF=Y @ END DEF
30 I1=-(M*N+P+F)/(N*((N-1)/2*M+P))
40 I2=.9*I1
50 DISP FNROOT(I1,I2,FNF(FVAR))

>RUN
n, fv, pv, pmt ? 168, 10925.76, 0, -45
4.79863130882E-3     -655.6829061
5.33181256535E-3     -1255.8749515
4.21615468581E-3     -41.4514876
4.17684617474E-3     -1.4838737
4.17538677149E-3     -.0035171
4.17538330417E-3     -.0000003
4.17538330387E-3     0
4.17538330387E-3

FYI, 1st column is interest rate, 2nd column is error in FV (thus function named FNF)
Example, I=0.479863%, calculated FV=11581.44, under-estimated true FV by -$655.68

Playing emulator interactively, this is unexpected ...

>FNROOT(I1, .5*I1, FNF(FVAR))
MATH ERR:Kybd FN in FNROOT/INTEGRAL

[J-F Garnier]

(12-20-2019 10:10 PM)Albert Chan Wrote:  Playing emulator interactively, this is unexpected ...
>FNROOT(I1, .5*I1, FNF(FVAR))
MATH ERR:Kybd FN in FNROOT/INTEGRAL

Doesn't come from the emulator, actually.

As the error message implies, it is not allowed to use a User-Defined function in FNROOT/INTEGRAL from the keyboard.

J-F

[Albert Chan]

(12-21-2019 03:44 PM)J-F Garnier Wrote:  

(12-20-2019 10:10 PM)Albert Chan Wrote:  Playing emulator interactively, this is unexpected ...
>FNROOT(I1, .5*I1, FNF(FVAR))
MATH ERR:Kybd FN in FNROOT/INTEGRAL

Doesn't come from the emulator, actually.

As the error message implies, it is not allowed to use a User-Defined function in FNROOT/INTEGRAL from the keyboard.

Thanks J-F Garnier.
A trick to get around this limitation is to wrap FNROOT inside an user function.

>PURGE
>10 INPUT "N, FV, PV, PMT ? ";N,F,P,M
>20 DEF FNF(I) @ Y=F+P+(M/I+P)*EXPM1(N*LOGP1(I)) @ DISP I,Y @ FNF=Y @ END DEF
>30 I1=-(M*N+P+F)/(N*((N-1)/2*M+P))
>40 I2=.9*I1
>50 DEF FNS(A,B)=FNROOT(A,B,FNF(FVAR))
>
>RUN
N, FV, PV, PMT ? 168, 10925.76, 0, -45
>I=FNS(I1,.5*I1)
2.66590628268E-3     1404.76181566
5.33181256535E-3     -1255.8749515
3.99885942402E-3     177.2179603
4.16369398215E-3     11.8480797
4.17461116337E-3     .7831244
4.75321186436E-3     -606.2697042
4.1753575833E-3       .0260875
4.46428472383E-3     -298.0104306
4.17538287345E-3     .0004365
4.31983379864E-3     -147.7508024
4.1753833002E-3       .0000036
4.24760854942E-3     -73.5647316
4.17538330374E-3     .0000001
4.21149592658E-3     -36.7050556
4.17538330384E-3     0

Update: this may improve rate guesses closer to true rate.

y = (1+I)^N-1 ≈ NI / (1 - ½ (N-1)I), Pade[1,1] approximation

Using this rough estimate of y for equation F + P + y(P + M/I) = 0
Solve for I, we have I = (F+P+M*N) / ((N-1)/2*(F-P) - P)

Above should be better than I1 estimate, which assumed no compounding.
Note: I1 is solution of F + P(1+NI) + m Σ(1+kI, k=0 to N-1) = 0

To make I2 still over-estimated true rate, take the arithmetic mean.
For the same reasoning, I1 can be reduced the same way.

>I1=(F+P+M*N) / (N*((1-N)/2*M - P))  → 5.33181256535E-3
>I2=(F+P+M*N) / ((N-1)/2*(F-P) - P)    → 3.68930884387E-3
>I3=(I1+I2)/2                                          → 4.51056070461E-3
>I4=(I1+I3)/2                                          → 4.92118663498E-3
>I=FNS(I3,I4)
4.51056070461E-3     -346.6853844
4.92118663498E-3     -790.3342398
4.18968078622E-3     -14.5135902
4.17566057222E-3     -.2812298
4.17538353441E-3     -.0002338
4.1753833039E-3       0

[Albert Chan]

This version use 3-points fit to interpolate for the root, instead of secant line.
Note: this assumed x = f(y), where f is a parabola, and interpolate for x as f(0)
Note: initial 2 rate guesses are arithmetic mean and harmonic mean of I1, I2.

10 INPUT "N, FV, PV, PMT ? ";N,F,P,M
20 DEF FNI(X1,Y1,X2,Y2)=X1-Y1*(X2-X1)/(Y2-Y1)
30 DEF FNF(I)=F+P+(P+M/I)*EXPM1(N*LOGP1(I))
40 I1=(F+P+M*N)/(N*((1-N)/2*M-P))
50 I2=(F+P+M*N)/((N-1)/2*(F-P)-P)
60 X1=(I1+I2)/2 @ Y1=FNF(X1) @ DISP X1,Y1
70 X2=I1*I2/X1 @ Y2=FNF(X2) @ DISP X2,Y2
80 X3=FNI(X1,Y1,X2,Y2)
90 Y3=FNF(X3) @ DISP X3,Y3
100 IF Y3=0 THEN END
110 E3=-Y3*FNI((X3-X2)/(Y3-Y2),Y2,(X3-X1)/(Y3-Y1),Y1)
120 X1=X2 @ Y1=Y2 @ X2=X3 @ Y2=Y3 @ X3=X3+E3
130 IF .00001*E3+X3-X3 THEN 90
140 DISP X3

>RUN
N, FV, PV, PMT ? 168, 10925.76, 0, -45
4.51056070461E-3     -346.6853844
4.36103281596E-3     -190.3488201
4.17897395274E-3     -3.6426516
4.17538425315E-3     -.0009629
4.1753833038E-3       0

>RUN
N, FV, PV, PMT ? 600, 1E6, 0, -250
1.08792431831E-2     -14144337.7743
4.93576250272E-3     78873.414844
4.96872148862E-3     65799.655312
5.13383450348E-3     -2862.19318
5.12721962668E-3     -6.6198
5.12720426664E-3     .00052
5.12720426785E-3

[Albert Chan]

This version extrapolate for rate, using Steffensen's method
I used HP71B quirky RES behavior on line-50, see HP71B RES bug?

Code:

10 INPUT "N, FV, PV, PMT ? ";N,F,P,M
20 DEF FNA(A,B,C)=C-(C-B)^2/(C-B-(B-A))
30 DEF FNI(I)=EXPM1(LOGP1(-I*(F+P)/(M+I*P))/N)
40 DISP 4*(F+P+M*N)/((N-1)*(F-M*N)-(3*N+1)*P)
50 DISP FNA(RES,FNI(RES),FNI(RES))
60 GOTO 50

To keep code short, it run until extrapolation hit by divide-by-0 (or other) error.

>DEFAULT OFF     ! error if overflow/underflow, divide-by-0 ...
>RUN
N, FV, PV, PMT ? 168, 10925.76, 0, -45
4.36103281596E-3
4.17891608926E-3
4.17538469191E-3
4.17538330384E-3
ERR L20:0/0

>RUN
N, FV, PV, PMT ? 600, 1e6, 0, -250
4.93576250272E-3
5.12771388025E-3
5.12720427126E-3
5.12720426785E-3
ERR L20:0/0

Note: above code require very good guess for rate.
Going for FV errors (see previous posts) is much better.

Example: tried this with default guess failed.
>RUN
N, FV, PV, PMT ? 40, 0, 25000, -1000
4.09556313993E-2
ERR L30:LOG(neg)

We need better guess (for this case, 2.1% ≤ guess ≤ 3.1%)
>.025 @ RUN 50
.025
.025252440621
.025243859111
2.52438486204E-2
ERR L20:0/0

[Albert Chan]

Some improvements over my 3-points fit rate finder (post #6)

• FNF may be more accurate, using Fused Multiply Add.
  FNF = FMA(EXPM1(N*LOGP1(I))/I, P*I+M, F+P)
• FNF(0) singularity removed, now FNF(0) = F+P + M*N
• We use harmonic mean of I1,I2 for X1.
• We use (0, FNF(0)) and (X1,Y1) to estimated X2.
  But, point "0" might be too far, giving bad interpolated rate.
  Instead of interpolation for Y=0 for X2, we cut the corrections in half (thus the factor 0.5)

Quote:10 INPUT "n, fv, pv, pmt ? ";N,F,P,M
20 DEF FNS(X1,Y1,X2,Y2)=X1-Y1*(X2-X1)/(Y2-Y1)
30 DEF FNF(I)
40 A=N @ IF I THEN A=EXPM1(N*LOGP1(I))/I
50 B=P*I+M @ A0=A*B
60 A1=A*1000001 @ A1=A+A1-A1 @ A=A-A1
70 B1=B*1000001 @ B1=B+B1-B1 @ B=B-B1
80 FNF=(F+P+A0)-(A0-A1*B1-A1*B-B1*A-A*B)
90 END DEF
100 X1=4*(F+P+M*N)/((N-1)*(F-M*N)-(3*N+1)*P)
110 Y1=FNF(X1) @ DISP X1,Y1 @ IF Y1=0 THEN END
120 X2=X1+.5*X1*Y1/(F+P+M*N-Y1)
130 Y2=FNF(X2) @ DISP X2,Y2 @ IF Y2=0 THEN END
140 X3=FNS(X1,Y1,X2,Y2)
150 Y3=FNF(X3) @ DISP X3,Y3 @ IF Y3=0 THEN END
160 E3=-Y3*FNS((X3-X2)/(Y3-Y2),Y2,(X3-X1)/(Y3-Y1),Y1)
170 X1=X2 @ Y1=Y2 @ X2=X3 @ Y2=Y3 @ X3=X3+E3
180 IF .00001*E3+X3-X3 THEN 150
190 DISP X3

>RUN
n, fv, pv, pmt ? 168, 10925.76, 0, -45
4.36103281596E-3     -190.348820115
4.24431571126E-3     -70.197349875
4.17612485668E-3     -.75216795
4.17538334438E-3     -.000041145
4.17538330381E-3
>RUN
n, fv, pv, pmt ? 600, 1e6, 0, -250
4.93576250272E-3     78873.4148425
5.18818567872E-3     -26668.4502875
5.12440307537E-3     1206.531985
5.12720670893E-3     -1.0521275
5.12720426781E-3     .0000175
5.12720426785E-3
>RUN
n, fv, pv, pmt ? 30,50000,1000,-1000
3.76850605653E-2     -930.964898075
3.68851990742E-2     -289.785021642
3.65236969622E-2     -3.12429013971
3.65197437154E-2     -1.49842044802E-4
3.65197435258E-2

[Albert Chan]

Again, basing the code on my 3-points fit rate finder (post #6), but with different guess.

\( \large {1\over I} ≈
{\binom{N}{3}M + \binom{N}{2}P \over \binom{N}{2}M + \binom{N}{1}P} -
{\binom{N}{2}M + \binom{N}{1}P \over F + P + M N}\)

10 INPUT "N, FV, PV, PMT ? ";N,F,P,M
20 DEF FNS(X1,Y1,X2,Y2)=X1-Y1*(X2-X1)/(Y2-Y1)
30 DEF FNF(I)
40 A=N @ IF I THEN A=EXPM1(N*LOGP1(I))/I
50 FNF=F+P+A*(P*I+M)
60 END DEF
70 C2=N*(N-1)/2 @ C3=C2*(N-2)/3
80 C3=C3*M+C2*P @ C2=C2*M+N*P @ C1=F+P+M*N @ C3=C3/C2
100 X1=1/(C3-C2/C1)
110 Y1=FNF(X1) @ DISP X1,Y1 @ IF Y1=0 THEN END
120 X2=1/(C3-C2/(C1+Y1))
130 Y2=FNF(X2) @ DISP X2,Y2 @ IF Y2=0 THEN END
140 X3=FNS(X1,Y1,X2,Y2)
150 Y3=FNF(X3) @ DISP X3,Y3 @ IF Y3=0 THEN END
160 E3=-Y3*FNS((X3-X2)/(Y3-Y2),Y2,(X3-X1)/(Y3-Y1),Y1)
170 X1=X2 @ Y1=Y2 @ X2=X3 @ Y2=Y3 @ X3=X3+E3
180 IF .00001*E3+X3-X3 THEN 150
190 DISP X3

>RUN
N, FV, PV, PMT ? 168, 10925.76, 0, -45
4.1171440611E-3       58.8707146
4.17253086536E-3     2.8926652
4.17539297797E-3     -.0098122
4.17538330384E-3     0
>RUN
N, FV, PV, PMT ? 600, 1e6, 0, -250
3.9653228628E-3      385578.34249
4.24281053755E-3     311571.332995
5.41104048363E-3     -130995.27126
5.15294183714E-3     -11161.26946
5.1267015499E-3      216.64945
5.12720400635E-3     .112712
5.12720426786E-3
>RUN
N, FV, PV, PMT ? 30,50000,1000,-1000
3.54153653967E-2     863.8049367
3.63975416852E-2     96.4591994
.03652100616            -.997799
3.65197435458E-2     -.0000155
3.65197435262E-2