Looking for TVM contributions

[Albert Chan]

(06-18-2024 03:27 AM)robve Wrote:  Let me illustrate with the "time reversed" TVM problem, which is essentially the same. The forward npmt() is best for this TVM example (the backward npmt() line is fatter):
./tvm 5 '?' -369494.09 17500 540000
3/0 i%= 11.9999997132 y=-2.05654941965e-08
4/0 i%= 11.9999997132 y= 0
4, 11.9999997132285614%, 0
s=0.567426862982966185 r=-0.43257313701703376

./tvm -5 '?' -540000 17500 369494.09
3/0 i%= 11.9999997132 y=-2.05764081329e-08
4/0 i%= 11.9999997132 y=-7.27595761418e-12
5/0 i%= 11.9999997132 y= 2.91038304567e-11
5, 11.9999997132285596%, 0
s=1.76234166063797959 r=0.762341660637979479

The |s| and |r| values are closer in the forward than in the backward npmt().

I think both time direction are equally good here.
1 ULP difference, may even be caused by scaling by 100 for percentage!

If I have to pick one, I would go for r = expm1(log1p(i)*n) to be positive. (== i*n > 0)
We can map positive to negative without too much error, but not the other way.

lua> i, n = .1, 100
lua> r1, r2 = EFF(i,n), EFF(i,-n)
lua> r1, r2
13779.612339822266      -0.9999274342840985
lua> -r2/(1+r2), -r1/(1+r1)
13779.612339813315      -0.9999274342840985

But, there are many things going on ... this is only 1 factor out of many.
See solve for n example, https://www.hpmuseum.org/forum/thread-21...#pid187533

[Albert Chan]

(06-17-2024 09:03 PM)robve Wrote:  Assuming TYPE is 0 for END and 1 for BGN, then we recognize PV and FV are switched and their signs inverted and N is inverted:
PV*S + PMT*K*R/I + FV = 0
with K=1+B*I and R=EXPM1(-N*LNP1(I)) and S=R+1

Compare this to

FV = (K*PMT*R/I-PV)/S
K*PMT*R/I - PV - FV*S = 0

This gives the same TVM problem.

Time reversed NFV turns into NPV, and vice versa.
We are solving same problem, but curve totally different.

(05-31-2024 10:39 PM)Albert Chan Wrote:  10 B=1 @ N=40 @ P=900 @ M=-40 @ F=-1000
20 DEF FNF(I) @ S=(1+I)^(-N) @ U=(1-S)/I @ K=1+B*I
30 FNF=P+K*M*U+F*S @ END DEF ! NPV
40 V=.01 @ W=FNF(V) @ I=.02 @ Y=FNF(I)
50 H=I-V @ V=I @ I=I-H*Y/(Y-W) @ W=Y @ Y=FNF(I)
60 DISP 100*I,Y @ IF V/I+I/V<>2 THEN 50

>run
 3.55877069046      -223.49575008
 4.3407684864        -68.532154113
 4.68660450765      -10.546204073
 4.74950349013      -.604877167
 4.75333056057      -.005731135
 4.75336716844      -.000003149
 4.75336718857      -.000000001

Let's turn NPV into NFV for rate problem, by time-reversal.

>15 VARSWAP P,F @ M=-M @ N=-N
>run
 11.3570289363       36923.539379
 2.35993812666      -1387.63023153
 2.68581267787      -1287.38125489
 6.87064131886       3585.77746844
 3.79135227821      -769.04731268
 4.33514429757      -379.62865086
 4.86526534953       115.24756927
 4.7418099072        -11.55575603
 4.75306056619      -.30740914
 4.75336803852       .00085219
 4.75336718851      -.00000007

Very unstable with multiple overshoot. NFV (n>0) formula for rate not recommended

[Albert Chan]

Hi, robve

I was looking at your test samples pdf. #23 just quit, returned NaN

Fix is simple, have Secant's initial 2 points guaranteed one-sided convergence.
Bonus: this allowed for simple no solution test, just like Newton's method.

1st point = edge (X = -PV or FV)  --> edge rate = PMT/X

I recently showed edge rate is really result of Newton's from *extreme* rates, -1 or ∞
Turns out, whatever edge X we pick, slope (from extreme) is -X !

2nd point = PMT/X + f(PMT/X) / X      // not as good as Newton on PMT/X, but good enough

Just as Newton's final half-correction, we do the same for Secants.

10 INPUT "B,N,P,M,F? ";B,N,P,M,F
20 IF B THEN P=P+M @ F=F-M ! end mode
30 IF ABS(P)>ABS(F) THEN X=F @ F=-P @ P=-X @ N=-N ! time reversed
40 DEF FNF(I) @ IF N*I*I+1-1 THEN FNF=((P+F)/EXPM1(LOGP1(I)*N)+P)*I+M @ END
50 Z=P+F @ FNF=(Z+N*M+I*((Z+N*(P-F))/2+I*Z*(N*N-1)/12))/N @ END DEF
60 X=F @ IF M/X<=-1 THEN X=-P ! guess = edge X
70 I0=M/X @ Y0=FNF(I0) @ I=I0+Y0/X @ Y=FNF(I) @ DISP 100*I0,Y0
80 DISP 100*I,Y @ Z=I-I0 @ I0=I @ I=I-Y*Z/(Y-Y0) @ Y0=Y @ Y=FNF(I)
90 IF SGN(Y)=SGN(Y0) AND ABS(Y)<ABS(Y0) THEN 80
100 DISP 100*I,Y @ IF Y-Y0 THEN J=I-Y*(I-I0)/(Y-Y0)/2 @ DISP 100*J

Code use small edge (LINE 60), but we can use the other, as long as edge rate is valid.

>run ! sample #23
B,N,P,M,F? 0,10,50,-30,80
-37.5                      .4474481846
-36.9406897692      .0348105026
-36.8935058326      .0000999938
-36.8933699054      .0000000229
-36.8933698743      .0000000002
-36.893369874        0
-36.893369874
>X=-P @ run 70 ! small P, big edge
 60                         .7159170954
 58.5681658092      .0493058621
 58.4622602744      .0001415644
 58.4619553284      .0000000289
 58.4619552661      -.0000000001
 58.4619552662

Although -X slope proof assumed positive N, time-symmetry make it valid for negative N too
Example, if we swap (P, F), code will do time-reversal to keep size of P small.

>run
B,N,P,M,F? 0,10,80,-30,50
37.5                      2.1052759718
34.8684050352      .291631744
34.4452478265      .0037390629
34.4397519889      .0000071623
34.4397414412      0
34.4397414412
>X=-P @ run 70 ! small P, big edge
-60                       .0081797505
-59.983640499      .0000312799
-59.983577699      .0000000002
-59.9835776986     0
-59.9835776986

Lets test tiny rate branch (LINE 50)

>run ! sample #21
B,N,P,M,F? 0,480,100000,-208.333333,0
 .208333333             121.444894064
 .086888438936        46.544537395
 1.14202213657E-2   5.774169845
 7.319118817E-4      .366932959
 6.607422444E-6      .003310945
 3.19310587E-9        1.93321242436E-6
-6.652600197E-10    1.03444729167E-11
-6.65280665633E-10 -0
-6.65280665633E-10

What about no solution case?
Y get worse no sign change, and I not converging.

B,N,P,M,F? 0,10,1000,-287,2000
-14.35                    116.137711015
-8.54311444925      61.534445887
-1.99911961932      20.995382798
 1.39004791475      8.609614417
 3.74593127745      3.290957549
 5.20365121116      1.292347407
 6.14624653622      .507826981
 6.75639679154      .207532994
 7.17807126641      .093719606
 7.52529900771      .056815067
 8.05986109467      .099136777
 6.80767250381

Update: I made a stupid mistake, using X for everything.
Now X refer only to edge, and I use Z for temp variable.

[Albert Chan]

Here is rate search method comparison, Secant vs Newton

Code:

              SECANT                                  NEWTON
>run ! 1
B,N,P,M,F? 1,40,900,-40,-1000
 4.16666666667       -5.1783387869       4.16666666667       -5.1783387869
 4.70607695697       -.4168082133        4.7524512828        -.0080716937
 4.75329510161       -.0006352866        4.75336718645       -.0000000187
 4.75336717995       -.000000076         4.75336718857        0
 4.75336718857        0                  4.75336718857
 4.75336718857

>run ! 2
B,N,P,M,F? 1,40,900,-400,-1000
 66.6666666667       -66.666666756       66.6666666667       -66.666666756
 77.7777777927       -11.111111044       79.9999999644       -.000000183
 79.9999999982       -.000000014         80.000000001         0
 80.000000001         0                  80.000000001
 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   -3.01507228575E-6     6.2977646195E-3
-1.1162507291E-6      2.4765911094E-3   -1.1099216589E-7      7.022277322E-4
 1.1441952785E-7      3.284869711E-4     3.07577857249E-7     1.52692719E-5
 3.02612985476E-7     2.32387146E-5      3.17093002025E-7     7.8839E-9
 3.16940255029E-7     2.527511E-7        3.17097920006E-7     0
 3.17097796093E-7     1.986E-10          3.17097920006E-7
 3.17097919979E-7     0
 3.17097919979E-7
>res * n                                >res * n
 10.0000000045                           10.0000000053

>run ! 4
B,N,P,M,F? 0,365*24*60,0,-0.01,5527.7829015
-1.80904354932E-4     6.2977549796E-3   -1.80904354932E-4     6.2977549796E-3
-6.6975224722E-5      2.4765920415E-3   -6.659604136E-6       7.022284099E-4
 6.8651170126E-6      3.284878032E-4     1.84546671641E-5     1.52693549E-5
 1.81567718641E-5     2.32388617E-5      1.90255801068E-5     7.8841E-9
 1.90164151578E-5     2.527541E-7        1.90258751938E-5    -6.E-14
 1.90258677529E-5     1.988E-10          1.90258751927E-5
 1.90258751936E-5    -7.E-14
 1.90258751923E-5
>res * n                                >res * n
 10.0000000011                           10.0000000013

>run ! 5
B,N,P,M,F? 0,365*24*60,0,-0.01,5260.382426
-1.90100247286E-4     .005827435363     -1.90100247286E-4     .005827435363
-7.9320556931E-5      2.2391630652E-3   -2.2550930229E-5      6.131837704E-4
-1.01915247867E-5     2.787870632E-4    -1.232943057E-7       1.15811690742E-5
-3.6061488282E-7      1.78257882E-5      3.16928059687E-7     4.4652E-9
 3.10916622228E-7     1.624902E-7        3.17097921452E-7    -8.E-14
 3.17094250469E-7     9.65E-11           3.1709791993E-7
 3.17097921431E-7    -1.E-14
 3.17097921241E-7
>res * 60*n                             >res * 60*n
 10.0000000443                           10.0000000029

>run ! 6
B,N,P,M,F? 0,32,-999999,0,1000000
 0                    .03125             0                    .03125
 .000003125           1.611328125E-8     3.12500161133E-6     0
 3.12500161133E-6     0                  3.12500161133E-6
 3.12500161133E-6

>run ! 7
B,N,P,M,F? 0,360,0,-1000,1000000
-.1                   2306.3389698      -.1                   2306.3389698
 .13063389698         1177.39713526      .313028391463        504.47068744
 .37116682905         327.89051248       .471270710313        62.96905732
 .464007148654        80.61001117        .497374290776        1.51994784
 .49427180422         8.66562343         .49803616688         .0009567
 .49791714914         .27395873          .498036584009        0
 .498036157001        .00097936          .498036584009
 .498036583962        .0000001
 .498036584006        .00000001
 .498036584011       -.000000003
 .49803658401

[robve]

(06-18-2024 06:15 PM)Albert Chan Wrote:  I was looking at your test samples pdf. #23 just quit, returned NaN

That's because guess_i() returns i%=101.94 that npmt cannot handle. If I go back to a small initial I% it works fine. I just left it in for now, knowing that iguess() should be changed later after investigating before deciding. By the way, this iguess() function is your suggested initial best guess.

It is the same test as on the contributions list of benchmarks #6.

- Rob

EDIT: I realize it may be a bit difficult to see the TVM parameters in the PDF. They are n, i% (to solve with ?), pv, pmt, fv below (the last one is a repeat which I neglected to take out):

./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
./tvm 328 ? 35000 -324 0
./tvm 58 ? -775 -49.56 4000
./tvm 32 ? -6000 0 10000
./tvm 48 ? 19198.60 -450 0
./tvm 5 ? -369494.09 17500 540000
./tvm 348 ? 243400 -1363.29 0
./tvm 30 ? -3200 -717.44 60000
./tvm 60 ? 243400 -1363.29 -222975.98
./tvm 24 ? 0 -50.26 1281.34
./tvm 6 ? -32000 0 28346.96
./tvm 60 ? 8000 -150.97 0
./tvm 456 ? 270000 -1215.3333333333 0
./tvm 480 ? 100000 -208.333333333333343 0
./tvm 10 ? 50 -30 400
./tvm 10 ? 50 -30 80
./tvm 10 ? -100 10 1e-10
./tvm 32 ? -999999 0 1e6
./tvm 525600 ? 0 -0.01 5260.38242600032663
./tvm 40 ? 900 -400 -1000 1

EDIT 2: the nan problem with iguess() returning a rate larger than 100% can be fixed with a test such as:

Code:

  i = iguess();
  if (i == 0.0)
    i = 1e-5;
  else if (fabs(i) >= 100)
    i = copysign(50, i);

Then we get an answer:
evals=7, i%=-36.8933698741777363%, FVerror=1.4e-14

But is +/-50% a reasonable new starting point to always avoid this? Don't think so. It looks like your idea to use edge is nice.

[robve]

(06-18-2024 11:54 PM)Albert Chan Wrote:  Here is rate search method comparison, Secant vs. Newton.

I'm no longer just comparing Secant vs. Newton, because we know Secant takes more steps. I'm comparing Hybrid to Newton.

I'm doing these important things:

- optionally allow one or more retries after overshooting, this is parameter g. When g=0 we stop. When g=1 we allow one retry to run iterations again after overshooting from the current I%.

- if the final iteration when we stop y=npmt() is nonzero, then we should pick the previous i%, not the last that overshot.

- I'm measuring FVerror to check the TVM model's accuracy with the i% solved, not the i% rate places because some final digits simply don't matter, it's not anomalous.

- I count npmt() evaluations, not steps, because Secant requires two initial i% which we count as 2 evaluations.

- for termination test I use for the Secant, Newton and Hybrid implementations:

Code:

    } while (y != 0.0 && fabs(y) < fabs(w))
  } while (y != 0.0 && y != w && g-- > 0);

My improved Newton implementation:

Code:

  i = iguess();
  if (i == 0.0)
    i = 1e-5;
  g = 1; // pick g=0 or g=1 (see PDF)
  do {
    do {
      w = y;
      y = npmt();
      if (y == 0)
        break;
      d = b*pmt - (pv + fv)*(1.0 + n*s/(r + r/j))/r - fv;
      h = 100*y/d;
      rr = r; // only to compute FV error later
      v = i;
      i -= h;
      if (isnan(i))
        return;
    } while (y != 0.0 && fabs(y) < fabs(w));
  } while (y != 0.0 && y != w && g-- > 0);
  if (y != 0) {
    i = v; // pick previous i% after overshooting
    r = rr; // only to compute FV error
    s = r + 1; // only to compute FV error
  }
  j = 0.01*i;
  printf("%.2g\n", fabs(fv - (k*pmt*r/j-pv)/s);

As you can see, it iterates until a root was found or when overshooting. When overshooting and g=1 then we continue iterating again until a root or overshooting. After overshooting we pick the previous i% which is closer to the root.

The results Newton with g=0 and g=1 versus Hybrid with g=0. Hybrid easily beats Newton:

  tvmperfcomp-newton-hybrid.pdf (Size: 34.65 KB / Downloads: 5)

I'm curious what your implementation does differently. Or if I'm doing something wrong. Perhaps the difference is that our convergence tests aren't the same? But from the first test in your list, it looks to me that it has the same number of evaluations: 1 initial npmt and 4 steps npmt = 5 evals. These evals are all nptm + derivative, so more expensive than the Hybrid method that only uses the derivative once after the initial guess.

PS. As I had mentioned before, Newton's derivate expression has subtle evaluation errors that get worse as we get closer to the root. Newton requires both npmt and derivative to be very precise (like having both a fat line to evaluate f(i%) and a sloppy slope to evaluate f'(i%)). Secant just takes the slope of two npmt points. Secant only uses npmt. In other words, it's less risky as we get closer to the root. The Newton convergence slows down anyway, no longer quadratic or better.

- Rob

[Albert Chan]

Hi, robve

Here is plot of f(x) = ((pv+fv)/((1+x)^n-1) + pv)*x + pmt, for test sample #23
I also showed asymptote line, f(x → ∞) = pv*x + pmt

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

And for the other root, also with asymptote line, f(x → -1) = -fv*x + pmt

plot {(130/((1+x)^10-1)+50)*x-30 , -80x-30}, x = -0.5 .. 0

Goal is to iterate rate from the outside, going in.
Edge rates are where asymptote line intersect to x-axis, thus always outside.
If we iterate from the inside, Newton will definitely overshoot ... perhaps by a lot.

Since asymptote line slope is more extreme (in size) than at any other point,
we can use it for Secant's 2nd point, and have a second guess, also from outside.

(06-19-2024 02:01 AM)robve Wrote:  - optionally allow one or more retries after overshooting, this is parameter g. When g=0 we stop. When g=1 we allow one retry to run iterations again after overshooting from the current I%.

I assumed when you say overshoot, it means f changes sign.
With correct guess for f (concave up or down), there is no need to retry after overshoot.

With wrong guess, Newton's next guess should overshoot back to the outside, so g=1
It is not as simple for Secant's method. It may take a while before *both* points are outside.

Quote:- if the final iteration when we stop y=npmt() is nonzero, then we should pick the previous i%, not the last that overshot.

If previous i% from the outside, overshoot is not really an overshoot, we should keep it.

If previous i% from the inside, we don't really know which is better.
All we know is true rate is somewhere between the two.

Quote:- I'm measuring FVerror to check the TVM model's accuracy with the i% solved, not the i% rate places because some final digits simply don't matter, it's not anomalous.

If i% getting very accurate (error of few ULP's), FVerror probably not that useful
Example, here we compare, for a given rate, calculated PMT vs FV

lua> n,i,pv,pmt,fv = 10,nil,50,-30,80 -- sample #23
lua> i = tvm(n,i,pv,pmt,fv)
lua> for k=0,5 do print(i, tvm(n,i,pv,nil,fv), tvm(n,i,pv,pmt,nil)); i=nextafter(i) end

Code:

-0.36893369874177734    -29.999999999999996     80.00000000000001
-0.3689336987417774     -30.000000000000004     80.00000000000001
-0.36893369874177745    -30.000000000000007     80
-0.3689336987417775     -30.00000000000001      79.99999999999999
-0.36893369874177756    -30.000000000000014     79.99999999999999
-0.3689336987417776     -30.00000000000002      79.99999999999997

Since rate search is solving npmt=0, it look like our solved i is good, so does (i+1ULP)
Using Plus42 for reference, true rate is (i+0.87ULP) = (i+1ULP) - 0.13 ULP

But from FV numbers, true rate = (i+2ULP)

I am not suggesting PMT_error is better test than FV_error.
It is just that solved i is too accurate for these tests to work.

From what I can tell, all your methods look equally good.

Quote:- for termination test I use for the Secant, Newton and Hybrid implementations:

Code:

    } while (y != 0.0 && fabs(y) < fabs(w))
  } while (y != 0.0 && y != w && g-- > 0);

From wrong guess, we don't know if y is always improving (toward 0).
Unless your guesses are from the outside, (fabs(y) < fabs(w)) is probably wrong.

Quote:... Newton's derivate expression has subtle evaluation errors that get worse as we get closer to the root.

Since both Newton and Secant are self-correcting, we can simply look at last 2 points.
With 2 points close together, Secant's slope is going to be much worse than Newton's.

With both methods use the *same* f, and Newton get better f', I would think Newton is better.
Of course, final f is tiny, and not very accurate. Perhaps more accurate f' make no difference.

[Albert Chan]

(06-19-2024 12:20 PM)Albert Chan Wrote:  asymptote line slope is more extreme (in size) than at any other point,

Here is a trick that work with any rate guess, and have one-sided convergence.

lua> n,pv,pmt,fv = 10,50,-30,80 -- sample #23
lua> f = fn'x: npmt(n,x,pv,pmt,fv)'
lua> i = 0.01 -- any guess OK

lua> X = -pv -- go for big root
lua> for k=1,5 do eps=f(i)/X; i=i+eps; print(i, eps) end
0.3514866009669545      0.3414866009669545
0.5527206785990931      0.20123407763213863
0.5821388175291544      0.029418138930061274
0.5844404930347874      0.002301675505633014
0.5846067021783057      0.0001662091435182589

We can speed up convergence, with smaller slope.
What is the ideal slope? We assume constant rate of convergence.

lua> r = 0.0001662091435182589 / 0.002301675505633014
lua> X = X * (1-r)
lua> for k=1,5 do eps=f(i)/X; i=i+eps; print(i, eps) end
0.5846195592477864      1.2857069480758773e-05
0.5846195526588186      -6.588967862434568e-09
0.5846195526622116      3.3930050292746237e-12
0.5846195526622099      -1.6848616523122454e-15
0.5846195526622099      -0

[robve]

I realize dumping a set of results with code isn't helping much to further this discussion. I'm happy to comment on questions and to avoid confusion.

(06-19-2024 12:20 PM)Albert Chan Wrote:  

(06-19-2024 02:01 AM)robve Wrote:  - optionally allow one or more retries after overshooting, this is parameter g. When g=0 we stop. When g=1 we allow one retry to run iterations again after overshooting from the current I%.

I assumed when you say overshoot, it means f changes sign.

No. See code of the termination test. It's when the magnitude of npmt increases rather than monotonically drops in magnitude. Sign is completely ignored. Sign changes are fine and can happen as an overshoot, but that can be an overshoot that is not necessarily detrimental. It can be auto-corrected and often will.

I don't understand why sign change should be treated differently. In fact, I initially tried your suggestion to take sign changes into account in the termination test. But the experimental results were worse for several tests in C double precision and on a 10 BCD calculator. The npmt sign can flip due to overshooting, but it is not necessarily bad and we should never stop right after the sign flips as was shown in the HP-71B code. It can be auto-corrected. So we want to keep iterating. This put me on the wrong path and I eventually dropped it from the termination tests.

(06-19-2024 12:20 PM)Albert Chan Wrote:  

Quote:- if the final iteration when we stop y=npmt() is nonzero, then we should pick the previous i%, not the last that overshot.

If previous i% from the outside, overshoot is not really an overshoot, we should keep it.

The term "overshoot" is confusing perhaps. An overshoot is not necessarily bad. I'm talking about non-monotonicity, i.e. when moving away from the root when the npmt magnitude increases. If we move away and we stop when that happens, then experiments show it is better to use the previous i% since we made a bad overshoot. A bad overshoot can be many orders of magnitude away from the root! There are plenty of examples, such as this TVM problem:

n=480 pv=100000 pmt=-208.333333333333343 fv=0
evals=4, i%=-1.39404957469160446e-17%, FVerror=5.4e+05

(06-19-2024 12:20 PM)Albert Chan Wrote:  If previous i% from the inside, we don't really know which is better.
All we know is true rate is somewhere between the two.

Yes. For example, I tested the use of a weighted average of i% after stopping with a sign change to judge where the root might be. But it didn't help to improve the result. Experimenting with this showed that taking the previous i% appears to reduce errors in i%. There is no theoretical argument why this should be bad. After all, it all is part of the solving sequence of points. But taking the last value when things go haywire is sure to be terrible. See the example above.

(06-19-2024 12:20 PM)Albert Chan Wrote:  

Quote:- I'm measuring FVerror to check the TVM model's accuracy with the i% solved, not the i% rate places because some final digits simply don't matter, it's not anomalous.

If i% getting very accurate (error of few ULP's), FVerror probably not that useful
Example, here we compare, for a given rate, calculated PMT vs FV

lua> n,i,pv,pmt,fv = 10,nil,50,-30,80 -- sample #23
lua> i = tvm(n,i,pv,pmt,fv)
lua> for k=0,5 do print(i, tvm(n,i,pv,nil,fv), tvm(n,i,pv,pmt,nil)); i=nextafter(i) end

Code:

-0.36893369874177734    -29.999999999999996     80.00000000000001
-0.3689336987417774     -30.000000000000004     80.00000000000001
-0.36893369874177745    -30.000000000000007     80
-0.3689336987417775     -30.00000000000001      79.99999999999999
-0.36893369874177756    -30.000000000000014     79.99999999999999
-0.3689336987417776     -30.00000000000002      79.99999999999997

Since rate search is solving npmt=0, it look like our solved i is good, so does (i+1ULP)
Using Plus42 for reference, true rate is (i+0.87ULP) = (i+1ULP) - 0.13 ULP

But from FV numbers, true rate = (i+2ULP)

I am not suggesting PMT_error is better test than FV_error.
It is just that solved i is too accurate for these tests to work.

Perhaps FV error, or PV error (also when "time reversed") or PMT error are all OK. When testing dozens of TVM cases I saw that FV error is more sensitive to an i% deviation.

By contrast, comparing i% digits to a higher precision i% solver is not meaningful when finding "better" values with a lower precision solver is impossible both theoretically and practically. A solver should just do its job well enough to limit the TVM model error to ensure TVM computations are reversible in most cases. Getting 1 or 2 ULP error in the other TVM parameters after solving to is attainable and should suffice to maintain TVM model equilibrium.

(06-19-2024 12:20 PM)Albert Chan Wrote:  

Quote:- for termination test I use for the Secant, Newton and Hybrid implementations:

Code:

    } while (y != 0.0 && fabs(y) < fabs(w))
  } while (y != 0.0 && y != w && g-- > 0);

From wrong guess, we don't know if y is always improving (toward 0).
Unless your guesses are from the outside, (fabs(y) < fabs(w)) is probably wrong.

There is no sign change test to stop at, so starting inside/outside shouldn't matter as much. We can't control the error in the slope to avoid a sign change in the first place. It will happen. The fabs(y) < fabs(w) was also in your code, so...

(06-19-2024 12:20 PM)Albert Chan Wrote:  

Quote:... Newton's derivate expression has subtle evaluation errors that get worse as we get closer to the root.

Since both Newton and Secant are self-correcting, we can simply look at last 2 points.
With 2 points close together, Secant's slope is going to be much worse than Newton's.

How? I don't see it in the experiments. You're assuming the two point npmt values have errors in different directions. Then it is worse. But expression evaluation errors in two close points are highly correlated, not independent. Claiming Secant is going to be much worse in general, i.e. assuming these errors are independent, is a stretch.

Furthermore, if the two points are very close so that expression evaluation for the two points return the same npmt value, then we should stop because the two points are very close to the root numerically and we can pick either point. Let me also add that self-correcting can also mean over-correcting when the correction suffers from inaccuracies, which the Newton derivative is not immune to, far from it

(06-19-2024 12:20 PM)Albert Chan Wrote:  With both methods use the *same* f, and Newton get better f', I would think Newton is better.
Of course, final f is tiny, and not very accurate. Perhaps more accurate f' make no difference.

I never said Secant is better than Newton. All I'm saying is that the Hybrid appears to be a promising approach (see code and results in the PDF). I'm not sure why you seem so skeptical when the evidence is there using a set of TVM problems I took from the HP forums (i.e. unbiased). Try the code if you don't believe the numbers.

I'm not completely happy about the iguess() (it's your guess_i() code by the way), which works fine in END mode, but is not as good in BEGIN mode. It also can return an overflow >100% causing NaN. I can replicate your symbolic derivation and understand it, but inserting a BEGIN mode adjustment that requires i% itself that we want to guess is something to figure out.

EDIT: I don't use "time reversal" when testing TVM rate solving, because I want to test and compare the solvers on whatever inputs are given, not on modified input that suits the solver. After all, the TVM problems are meant to test TVM rate solving under strain.

- Rob

EDIT: fixed typo

[Albert Chan]

(06-19-2024 03:47 PM)robve Wrote:  n=480 pv=100000 pmt=-208.333333333333343 fv=0
evals=4, i%=-1.39404957469160446e-17%, FVerror=5.4e+05

(-pmt) > (pv/n), by a smudge. I would expect positive rate (0 is good too)
Still, -1.4e-17% is not too bad (using fma, true rate = +1.89e-17%)

i ≈ 1/p - p/n² = (n²-p²)/(p*n²) = (n+p)*(n-p) / (p*n²)      -- where p=-pv/pmt

lua> fma(n,pmt,-pv) * fma(n,pmt,pv) / (-pv*pmt*n*n)      -- scale by pmt²/pmt²
1.8947806286936e-19

More concerning is huge FVerror, perhaps bug in FV code?

n=480 pv=100000 pmt=-208.333333333333343 fv=-5.4e+5 --> solve rate = 0.4598% ?

(06-19-2024 12:20 PM)Albert Chan Wrote:  Since both Newton and Secant are self-correcting, we can simply look at last 2 points.
With 2 points close together, Secant's slope is going to be much worse than Newton's.

robve Wrote:How? I don't see it in the experiments. You're assuming the two point npmt values have errors in different directions. Then it is worse. But expression evaluation errors in two close points are highly correlated, not independent. Claiming Secant is going to be much worse in general, i.e. assuming these errors are independent, is a stretch.

I said Secant's slope is worse, not Secant's method.
Here is an example, Newton's slope with 15 good digits, Secant's slope only 1

Assuming last point from Newton and Secant reach about the same place.

Newton final step: x - y / y'
Secant final step: x - y / ((y-y0)/(x-x0))

Also, you missed my conclusion. Bad slope might not matter.

(06-19-2024 12:20 PM)Albert Chan Wrote:  With both methods use the *same* f, and Newton get better f', I would think Newton is better.
Of course, final f is tiny, and not very accurate. Perhaps more accurate f' make no difference.

robve Wrote:I'm not completely happy about the iguess() (it's your guess_i() code by the way), which works fine in END mode, but is not as good in BEGIN mode. It also can return an overflow >100% causing NaN.

guess_i() is retired from use. Edge rate is now preferred.

Anyway, BEGIN mode {n,pv,pmt,fv} == END mode {n,pv+pmt,pmt,fv-pmt}

test sample #1 Wrote:./tvm 40 ? 900 -40 -1000 1

lua> n,pv,pmt,fv = 40,900,-40,-1000 -- BEGIN mode
lua> pv,fv = pv+pmt, fv-pmt -- END mode
lua> n,pv,pmt,fv
40      860      -40      -960
lua> pmt/fv, pmt/-pv
0.041666666666666664      0.046511627906976744

1 sign changes --> 1 solution, either edge is good.

[Albert Chan]

robve Wrote:The npmt sign can flip due to overshooting, but it is not necessarily bad and we should never stop right after the sign flips as was shown in the HP-71B code. It can be auto-corrected. So we want to keep iterating...

I think we are solving rate problem with different approach.

Your code seems to treat f(x) like a black-box, and solve it with Secant/Hybrid method.
Code can be easily adapted for solving other f(x)=0 problem, say, change npmt to NPV
Rate guess can be changed, even user supplied. But it is not guaranteed to work.

Example, test sample #23 were designed to fail Plus42 rate solver
Plu42 had hard coded guess=0.0001, next iteration give an impossible rate of -200%

Instead of general solver, my approach is specialized rate solver.
I wanted rate solver to always work, even if more iterations required.
Even harder, it had to figured out no solution case, and let user know.

Before start coding, I proved npmt(x) is either concave up (or down).
With this out of the way, edge guess from the outside, I get one-sided convergence.



Technically this setup should never encounter overshoot (f sign flip)
If we do see overshoot, it is due to rounding error (f = 0), rate is found.
No more auto-correction needed. It cannot get any better.

If f not improving (toward 0), and i not converging --> no solution.

[robve]

(06-19-2024 02:01 AM)robve Wrote:  My improved Newton implementation:
...

Previously in that post I showed the results for 64 bit float (double). I also have results for 128 bit floats (long double) and updated to the most recent implementations:

  tvmperfcomp-newton-hybrid-128bit.pdf (Size: 34.47 KB / Downloads: 3)
Hybrid has lower FVerror overall. It is hard to argue that the difference in accuracy is due to roundoff alone.

The implementations aren't using any tricks or special prep. They take the same initial guess, then Hybrid performs one Newton step to get the next point and iterates Secant until convergence without retrying (g=0). Two case for Newton are shown (g=0) and with one retry (g=1) after overshooting to try to find the root.

- Rob

[robve]

(06-19-2024 09:06 PM)Albert Chan Wrote:  

(06-19-2024 03:47 PM)robve Wrote:  n=480 pv=100000 pmt=-208.333333333333343 fv=0
evals=4, i%=-1.39404957469160446e-17%, FVerror=5.4e+05

(-pmt) > (pv/n), by a smudge. I would expect positive rate (0 is good too)
Still, -1.4e-17% is not too bad (using fma, true rate = +1.89e-17%)

i ≈ 1/p - p/n² = (n²-p²)/(p*n²) = (n+p)*(n-p) / (p*n²)      -- where p=-pv/pmt

lua> fma(n,pmt,-pv) * fma(n,pmt,pv) / (-pv*pmt*n*n)      -- scale by pmt²/pmt²
1.8947806286936e-19

More concerning is huge FVerror, perhaps bug in FV code?

No. The error is in the i% result. As a consequence the TVM model's parameter space suffers large errors.

With previous i% after overshooting the TVM model parameter space isn't badly affected
./tvm 480 '?' 100000 -208.333333333333343 0
1/1 i%= 1e-05 y= 0.00501045666647
2/1 i%= 7.98321085531e-11 y= 3.999917908e-08
3/1 i%= 1.03199085624e-16 y= 2.84217094304e-14
4/1 i%= 6.09287401526e-17 y= 2.84217094304e-14
4, i%=6.0928740152605195e-17%
FVerror=1.5e-11
PVerr=1.5e-11
PMTerr=2.8e-14

With last i% after overshooting the TVM model parameters space error is bad:
./tvm 480 '?' 100000 -208.333333333333343 0
1/1 i%= 1e-05 y= 0.00501045666647
2/1 i%= 7.98321085531e-11 y= 3.999917908e-08
3/1 i%= 1.03199085624e-16 y= 2.84217094304e-14
4/1 i%= 6.09287401526e-17 y= 2.84217094304e-14
4, i%=-1.39404957469160446e-17%
FVerr=5.4e+05
PVerr=5.4e+05
PMTerr=2.6e+02

Which result should we trust? I am inclined to accept the former over the latter.

(06-19-2024 09:06 PM)Albert Chan Wrote:  

robve Wrote:I'm not completely happy about the iguess() (it's your guess_i() code by the way), which works fine in END mode, but is not as good in BEGIN mode. It also can return an overflow >100% causing NaN.

guess_i() is retired from use. Edge rate is now preferred.

Anyway, BEGIN mode {n,pv,pmt,fv} == END mode {n,pv+pmt,pmt,fv-pmt}

Yes, but for iguess() adjusted for pmt this produces bad initial i% guess for some of the more extreme cases and that worries me. I don't understand why and how an adjustment can be safely made. But I digress since guess_i() is retired (it's quite good actually, except for BEGIN mode and a few other cases). For example, I noticed that adjusting pmt for BGN mode in this case n=40 pv=900 pmt=-400 fv=-1000 gives i%=371.725. Now, this is a hard example so it's not unexpected but still.

(06-19-2024 09:06 PM)Albert Chan Wrote:  

test sample #1 Wrote:./tvm 40 ? 900 -40 -1000 1

lua> n,pv,pmt,fv = 40,900,-40,-1000 -- BEGIN mode
lua> pv,fv = pv+pmt, fv-pmt -- END mode
lua> n,pv,pmt,fv
40      860      -40      -960
lua> pmt/fv, pmt/-pv
0.041666666666666664      0.046511627906976744

1 sign changes --> 1 solution, either edge is good.

With BEGIN mode and the Newton solver:
1/1 i%= 4.55232470519 y= -1.77265084526
2/1 i%= 4.75326326796 y=-0.000915829625612
3/1 i%= 4.75336718854 y=-2.40206077251e-10
4/1 i%= 4.75336718857 y= -7.1054273576e-15
i%=4.75336718857102181%
FVerr=1.1e-13
PVerr=0
PMTerr=7.1e-15

After adjusting PV and FV for BEGIN mode there is a sign change because we start at "the wrong side" of the curve but that's no problem:
1/1 i%= 4.75977785174 y= 0.0564948432948
2/1 i%= 4.75336708491 y=-9.13549698112e-07
3/1 i%= 4.75336718857 y= 7.1054273576e-15
4/1 i%= 4.75336718857 y= -7.1054273576e-15
i%=4.75336718857102092%
FVerr=2.3e-13
PVerr=0
PMTerr=7.1e-15

Both are equally good.

When we have to do one more iteration to get on the right side of the curve, then that isn't a big deal.

About the npmt curvature: adjusting npmt for positive f'' is just a sign change of npmt. But root finders don't need to care as we know, because they flip to the other side of the curve anyway (as shown above). Theoretically then they never need to flip back again unless f'' changes sign. For TVM f'' won't change sign I believe, but numerically there is a nonzero possibility that the sign of npmt may flip close to root due to inaccuracies, especially with lower precision solvers. That worries me.

- Rob

[Albert Chan]

(06-20-2024 04:03 AM)robve Wrote:  With last i% after overshooting the TVM model parameters space error is bad:
./tvm 480 '?' 100000 -208.333333333333343 0
1/1 i%= 1e-05 y= 0.00501045666647
2/1 i%= 7.98321085531e-11 y= 3.999917908e-08
3/1 i%= 1.03199085624e-16 y= 2.84217094304e-14
4/1 i%= 6.09287401526e-17 y= 2.84217094304e-14
4, i%=-1.39404957469160446e-17%
FVerr=5.4e+05
PVerr=5.4e+05
PMTerr=2.6e+02

Because f(ε) was not implemented, y3 and y4 almost the same.
Bad y/y' cause extrapolated i% go crazy, even its sign is wrong.

f(ε) part is a 1-liner. I think it is worth adding.

40 DEF FNF(I) @ IF N*I*I+1-1 THEN FNF=((P+F)/EXPM1(LOGP1(I)*N)+P)*I+M @ END
50 Z=P+F @ FNF=(Z+N*M+I*((Z+N*(P-F))/2+I*Z*(N*N-1)/12))/N @ END DEF

But bad i does not explain huge errors. How are they calculated?

lua> i = -1.3940495746916044e-19
lua> NFV(480, i, 1e5, -208.333333333333343, 0)
0
lua> NPV(480, i, 1e5, -208.333333333333343, 0)
0
lua> npmt(480, i, 1e5, -208.333333333333343, 0)
0

robve Wrote:iguess() adjusted for pmt this produces bad initial i% guess for some of the more extreme cases and that worries me ...
BGN mode in this case n=40 pv=900 pmt=-400 fv=-1000 gives i%=371.725

guess_i() assumed (n*i) is small, and NPMT is quadratic of rate.
Above example, if we solve actual quadratic, we have imaginary numbers.
Pade approx was used to return real guess, even if quadratic have no real root.

BTW, why is i% > 100 a problem?
I just tried 1 billion %, and it run OK (next iteration ≈ edge rate, as expected)

lua> tvm_begin(40, 1e7, 900, -400, -1000, true)

Code:

1e+07                   4999999600              500
0.8000000007450581      3.6760650345968315e-07  500.00000010323674
0.8000000000098451      0                       500.00000010323674
0.8000000000098451

robve Wrote:About the npmt curvature: adjusting npmt for positive f'' is just a sign change of npmt. But root finders don't need to care as we know, because they flip to the other side of the curve anyway (as shown above). Theoretically then they never need to flip back again unless f'' changes sign. For TVM f'' won't change sign I believe, but numerically there is a nonzero possibility that the sign of npmt may flip close to root due to inaccuracies, especially with lower precision solvers. That worries me.

Problem is in Secant's method. We are estimating slope with points, which may not be both on the outside.
Depends on where they are, we don't know how many iterations until *both points* moved outside.

Hybrid work better perhaps is due to this. It moved next point to the outside, on the first try.
Switch back to Secant with both points outside, give us nice one-sided convergence.

(04-09-2022 05:47 PM)Albert Chan Wrote:  \(\displaystyle I ≈ \frac{1}{P} - \frac{P}{N^2}\)

This is even better, and work well with big N.

find_rate() thus skip over user-supplied guess, and consider next one as initial guess.
Below example, the better guess is on the inside ... Newton pushed it out.

lua> n,pv,pmt,fv = 365*24*60*60,0,-0.01,331667.0067 -- problem 3
lua> i = pmt/fv
lua> i -- edge rate
-3.015072285753529e-08

lua> find_rate(n, i-1/i/n^2, pv, pmt, fv, true) -- perhaps better guess

Code:

3.1987564370323693e-09  -4.4522291350098525e-06 -160259.18297182902
3.1709750078970886e-09  6.720436082630066e-10   -160307.5638939386
3.1709792001110547e-09  1.5612511283791264e-17  -160307.55659320104
3.170979200111152e-09   0                       -160307.55659320077
3.170979200111152e-09

[robve]

(06-19-2024 09:08 PM)Albert Chan Wrote:  Technically this setup should never encounter overshoot (f sign flip)
If we do see overshoot, it is due to rounding error (f = 0), rate is found.
No more auto-correction needed. It cannot get any better.

If f not improving (toward 0), and i not converging --> no solution.

Well, try this example with Newton. The curve is parabolic with f''>0. We stop at sign change and take the last i%:
./tvm 10 '?' 50 -12.8277219567 80
1/0 i%= 3.08778325056 y= 0.0103695628259
2/0 i%= 3.59460137697 y= 0.00259881409493
3/0 i%= 3.8488202702 y= 0.00064807442555
4/0 i%= 3.97517959151 y= 0.000159382918897
5/0 i%= 4.03630943599 y= 3.72170552208e-05
6/0 i%= 4.06308588579 y= 7.13292596899e-06
7/0 i%= 4.07140713916 y= 6.88560398032e-07
8/0 i%= 4.07240259392 y= 9.85261117137e-09
9/0 i%= 4.07241725752 y= 2.14050999148e-12
10/0 i%= 4.07241726071 y= -3.5527136788e-15
10, i%=4.07241726070281285%, FVerr=2e-10

This is not a winning strategy IMHO. Newton stopping at non-monotonicity is better with previous i% when the last step wasn't monotonically decreasing:

./tvm 10 '?' 50 -12.8277219567 80
1/0 i%= 3.08778325056 y= 0.0103695628259
2/0 i%= 3.59460137697 y= 0.00259881409493
3/0 i%= 3.8488202702 y= 0.00064807442555
4/0 i%= 3.97517959151 y= 0.000159382918897
5/0 i%= 4.03630943599 y= 3.72170552208e-05
6/0 i%= 4.06308588579 y= 7.13292596899e-06
7/0 i%= 4.07140713916 y= 6.88560398032e-07
8/0 i%= 4.07240259392 y= 9.85261117137e-09
9/0 i%= 4.07241725752 y= 2.14050999148e-12
10/0 i%= 4.07241726071 y= -3.5527136788e-15
11/0 i%= 4.0724172607 y= 1.7763568394e-15
12, i%=4.07241726070545784%, FVerr=2.8e-14

With the Hybrid method we have no sign change and we hit the root even closer:
./tvm 10 '?' 50 -12.8277219567 80
3/0 i%= 3.7640993412 y= 0.00115435116452
4/0 i%= 3.89955465008 y= 0.000413574938126
5/0 i%= 3.97517927056 y= 0.000159383755689
6/0 i%= 4.02259765825 y= 5.81463133074e-05
7/0 i%= 4.04983268476 y= 2.02404140044e-05
8/0 i%= 4.06437522885 y= 6.04424116091e-06
9/0 i%= 4.07056694296 y= 1.27674678296e-06
10/0 i%= 4.0722250993 y= 1.29426121376e-07
11/0 i%= 4.07241215141 y= 3.43175443618e-09
12/0 i%= 4.07241724621 y= 9.73088276623e-12
13/0 i%= 4.0724172607 y= 1.7763568394e-15
14/0 i%= 4.0724172607 y= 1.7763568394e-15
14, i%=4.07241726070100896%, FVerr=1.4e-14

Note how many trailing digits of i% are actually irrelevant. This is an example why comparing decimal places of i% is just wrong in general, the point I made earlier. The reason for this is the flat curvature, hence more iterations are needed than usual to get to the root for this example. It depends on the TVM problem how many places are relevant. The relationship is nonlinear.

When we increase PMT a little bit then there is one root numerically (two roots very very close numerically indistinguishable) on the x-axis:
./tvm 10 '?' 50 -12.8277106125 80
For this specific case we can stop at a sign change without losing precision.

IMPORTANT: all such examples are anecdotal. Anecdotal examples are fine as counter examples to prove claims are incorrect (counter to a proof), but shall never be solely used to offer a general proof.

- Rob

[robve]

(06-20-2024 03:07 PM)Albert Chan Wrote:  

(06-20-2024 04:03 AM)robve Wrote:  ...
FVerr=5.4e+05
PVerr=5.4e+05
PMTerr=2.6e+02

But bad i does not explain huge errors. How are they calculated?

lua> i = -1.3940495746916044e-19
lua> NFV(480, i, 1e5, -208.333333333333343, 0)
0
lua> NPV(480, i, 1e5, -208.333333333333343, 0)
0
lua> npmt(480, i, 1e5, -208.333333333333343, 0)
0

I'm literally using the same parameters to find the root at npmt()=0 and to compute the errors:

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;
}

After root finding with npmt()=0 all global parameters are set accordingly. I only need to make sure to use the previous j, r, s when I use the previous i% rate that is more precise when the last step overshot. Then fabs(fv - (k*pmt*r/j-pv)/s) is the FV error. It's just the TVM FV formula.

The npmt() is zero or close of course, because we solve using npmt. Also nfv(0 and npv probably are also zero or close. But the TVM reverse computation of PV, PMT and FV should not have large errors like that. I really want to maintain TVM model equilibrium, i.e. inputs and outputs are reversible, at least with minimal drift, in the presence of roundoff.

- Rob

[robve]

(06-19-2024 03:47 PM)robve Wrote:  I'm not completely happy about the iguess() (it's your guess_i() code by the way), which works fine in END mode, but is not as good in BEGIN mode. It also can return an overflow >100% causing NaN. I can replicate your symbolic derivation and understand it, but inserting a BEGIN mode adjustment that requires i% itself that we want to guess is something to figure out.

After some more experimentation with TVM problems for which the initial guess isn't that good, such as
n=10 pv=-100 pmt=10 fv=<small> for <small>=1e-1 to 1e-99
n=40 pv=900 pmt=-400 fv=-1000 BGN
n=480 pv=100000 pmt=-208.333333333333343 fv=0
and others, the following updated version is a nice improvement:

Code:

double iguess()
{
  double P = pv + b*pmt;
  double F = fv - b*pmt;
  double A = (F + P)/n;
  double B = (F - P) - A;
  double C = (pmt + A)/B;
  double D = (n*n - 1.0)*A*C/B;
  double I = 100.0*C*(D - 3.0)/(D - 1.5);
  if (fabs(I) < 1e-16)
    I = copysign(1e-8, I);
  else if (fabs(I) > 99)
    I = copysign(99, I);
  return I;
}

It is based on Albert's guess_i(). It may not start at a point to approach the root from a side that does not flip the npmt sign (it doesn't affect these root finders). It does appear to offer a fairly accurate starting point. In the Hybrid method the next point we pick after the guess is often quadratically closer to the root using the Newton derivative. That's all we need for some TVM problems as it hits the root exactly with only 1 to 3 npmt evaluations.

Results for Hybrid with the new guess compared to the previous implementation looks good, saving 5 iteration steps total over all tests: 98 instead of 103, not counting the NaN case that is now fixed. Final accuracy is the same.

  tvmperfcomp-hybrid-newguess.pdf (Size: 33.22 KB / Downloads: 2)

- Rob

EDIT: minor edit to correct typo and I=1e-8 that works for Newton when I=1e-9 fails with NaN in Newton.

[Albert Chan]

(06-20-2024 04:30 PM)robve Wrote:  Well, try this example with Newton. The curve is parabolic with f''>0. We stop at sign change and take the last i%:
./tvm 10 '?' 50 -12.8277219567 80
...
10, i%=4.07241726070281285%, FVerr=2e-10

This is not a winning strategy IMHO. Newton stopping at non-monotonicity is better with previous i% when the last step wasn't monotonically decreasing:

12, i%=4.07241726070545784%, FVerr=2.8e-14

With the Hybrid method we have no sign change and we hit the root even closer:

14, i%=4.07241726070100896%, FVerr=1.4e-14

Solve with mpmath, mp.dps = 34, rounded back to float
rate1(%) = 4.14000005643961 -- ignored for error analysis
rate2(%) = 4.0724172607043148

Newton error =  1.502e-12 %
Secant error = -1.143e-12 %
Hybrid error =  3.306e-12 %

[robve]

I've expanded the comparisons by including another Newton method. The following results compare Newton and Hybrid to another Newton solver that was mentioned on the HP Forums before. To ensure fairness, I used the exact same initial i% guess that performs very well thanks to Albert and I've used the same termination condition for all four methods compared. To emulate BEGIN mode I added PV+=PMT and FV-=PMT to solve rate i% with this other Newton solver (PV and FV updates are undone to compute FVerror.)

  tvmperfcomp-voidware-newton-hybrid.pdf (Size: 36.22 KB / Downloads: 4)

This other Newton solver performs objectively worse than the Newton and Hybrid methods used and discussed in this thread. There is one minor benefit of this method. When we disregard the problematic case that requires as many as 53 evaluations to converge (!), this Newton method requires slightly fewer evaluations on average. However, the results are less accurate and some are wrong. Perhaps interesting to note is that additional iterations (with g=1 to continue with one retry) may improve the results in some cases, but mostly harms the result. Interestingly, this is the opposite of the Newton method used and discussed in this thread. I suppose this aligns with the claims that Newton should not need continuation theoretically after the termination condition was met, but depending on the choice of f and f' expressions may benefit from additional iterations.

- Rob

[robve]

The Newton versus Hybrid method performance on a SHARP PC-1403H with 10 digits BCD. It took some effort to implement, double check and run the numbers by hand.

  tvmperfcomp-pc1403h-newton-hybrid.pdf (Size: 33.56 KB / Downloads: 5)

I've used g=1 for both methods to continue one more phase of iterations after one time "overshooting" (an overshoot means that the last step was too large, i.e. it is an "overstep" that is likely moving away from the root). Otherwise, FVerror will be larger for some cases and setting g=1 does improve these cases. However, I should mention that overshooting does not happen a lot. It happens in a couple of cases such as the 6th and 21st test case. But not 22nd, even though it is similar to 21st and also has a nonzero FVerror. After the first overshoot, we iterate again. Then, if an overshoot happens again, we take the previous i% and stop, because this point is closer to the root.

- Hybrid evals are cheaper, i.e. no derivative to compute

- Newton required 93 evals versus 103 for Hybrid

- The 3rd case FVerror=100 looks really bad, but it means that we computed FV=-900 instead of FV=-1000 which was discussed in this thread before why this happens. It is very hard (impossible?) to improve this result with only 10 digits BCD precision as the least significant digit in i% causes FV to sway.

- Sometimes Newton is more accurate and sometimes Hybrid is more accurate as indicated with comparatively larger FVerror shown in red. However, Newton method errors appear worse, such as in the 6th, 21st and 22nd and last TVM test cases. Hybrid method errors are not as large.

My takeaway from these experiments with 10 digits BCD, IEEE 754 double and long double precision is that I personally prefer the Hybrid method, because Newton behaves erratic when rates approach zero and does not appear to maintain TVM equilibrium as well as Hybrid does (i.e. lower FV, PV and PMT errors in reversed computations).

This is all experimental using unbiased test cases collected on this forum. If something changes significantly in the future when implementing and testing these methods, then I will be happy to change my mind.

The TVM test cases are (I've removed the last three that were duplicates left after merging test sets):
n pv pmt fv (BGN mode = 1)
40 900 -40 -1000 1
40 900 1000 -1000 1
40 900 -400 -1000 1
31536000 0 -0.01 331667.0067
525600 0 -0.01 5282.36776893622391
525600 0 -0.01 5260.38242600032663
32 -999999 0 1000000
360 0 -1000 1000000
328 35000 -324 0
58 -775 -49.56 4000
32 -6000 0 10000
48 19198.60 -450 0
5 -369494.09 17500 540000
348 243400 -1363.29 0
30 -3200 -717.44 60000
60 243400 -1363.29 -222975.98
24 0 -50.26 1281.34
6 -32000 0 28346.96
60 8000 -150.97 0
456 270000 -1215.3333333333 0
480 100000 -208.333333333333343 0
10 50 -30 400
10 50 -30 80
10 -100 10 1e-10

The PC-1403H program:

Code:

' 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

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 25
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
25 G=1,C=0 : IF P+F=-N*M LET I=0 : RETURN
' solve Y=NPMT=0 for rate I% pick initial nonzero best guess and perform newton step
26 V=P+B*M,W=F-B*M
27 Y=(V+W)/N,W=W-V-Y,V=(M+Y)/W,W=(N*N-1)*Y*V/W,I=100*V*(W-3)/(W-1.5)
28 IF ABS I<1E-9 LET I=1E-5*(SGN I+(I=0))
29 IF ABS I>99 LET I=99*SGN I
30 GOSUB 38 : IF Y=0 RETURN

' USE THIS PART FOR HYBRID:
31 V=I,I=I-(K*M-(P+F*S)*J/R)*100/(B*M-(P+F)*(1+N*S/(R+R/J))/R-F)
32 GOSUB 38 : IF Y=0 OR Y=W RETURN
' secant loop
33 T=I,I=I-(I-V)*Y/(Y-W),V=T
34 GOSUB 38 : IF Y<>0 IF ABS Y<ABS W GOTO 33
35 IF Y<>0 IF Y<>W IF G LET G=G-1 : GOTO 33

' USE THIS PART FOR NEWTON:
31 V=I,I=I-(K*M-(P+F*S)*J/R)*100/(B*M-(P+F)*(1+N*S/(R+R/J))/R-F)
34 GOSUB 38 : IF Y<>0 IF ABS Y<ABS W GOTO 31
35 IF Y<>0 IF Y<>W IF G LET G=G-1 : GOTO 31

' pick previous I% after overstepping
36 IF Y<>0 LET I=V
' round very small I% to zero
37 I=(ABS I>1E-9)*I : RETURN
' compute Y=NPMT and save previous Y in W
38 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

To count evals and compute the FVerror, add Z=Z+1,Q=R at the start of line 38 then set Z=0 at line 25. Change line 36 to 36 IF Y<>0 LET I=V,J=.01*I,R=Q,S=R+1 to pick the previous rate and restore the previous TVM parameters. Compute FVerror with ABS(F-(K*M*R/J-P)/S)

- Rob

EDIT: minor update to the PDF so all significant i% digits are now visible
EDIT 2: minor update to the PDF with a minor update to the program to test P+F=-N*M to return zero rate, which does not apply to the last test case