Request TVM formula guidance

[iMatt04bit]

I am writing code for an HP12C calculator sim as a personal project for learning Java.
All gui components are complete and working, as are all scientific, math, statistics, dates and all buttons except the TVM buttons which are in process.

Now its down to the code for solving compound interest for amortized loans and I have had good results, including the capability for solving 2 roots. My only concern is that I have actually written a good compliment of f(x) where f'(x) is the derivative.

The results appear precise and a match to a variety of real and simulated calculators.
I just need to know if I have a correct derivative. ...if only I had taken Calculus...

Could someone take a look, if they like, and let me know if I am on the right track.

This is psuedo-Java-code of what I am using, residing within a loop.

K = 0 (month end interest computation, else 1 for beginning)
PYR = number of annual compounding periods x years, ie: 12
I = initial guess followed by storing computed guesses (display format 5 for 0.05)
R = I / 100 / PYR
TVM values supplied: N,PMT,PV,FV,PYR,K
Solve for I using Newton's method in iterative code.

// formula from the hp12c manual, solves to 0 for matched I converted to R, used below
f(x) = PV + (1+R*K) * PMT * ((1-(1+R)^-N)/R) + FV * (1+R)^-N

// the derivative as best I know
f'(x) = PV + (1+R*K) * PMT * [1/N * (1-(1+R)^-(N-1)) / (1/R)] + FV * [1/N * (1+R)^-(N-1)]

// compute for +I
I = I - f(x)/f'(x)

// or for -I
I = I + f(x)/f'(x)

// breakpoint on difference of prior computed I and current I is <= 5e-12 for fractional I | I < 1 percent
// or computed I <= 5e-15 for I | I >= 1 percent

Thanks

[Thomas Klemm]

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // formula from the hp12c manual, solves to 0 for matched I converted to R, used below
f(x) = PV + (1+R*K) * PMT * ((1-(1+R)^-N)/R) + FV * (1+R)^-N

D[P+(1+R*K)*T*((1-(1+R)^(-N))/R)+F*(1+R)^(-N),R]

I just had to use single variable names P, T and F instead of PV, PMT and FV and use parentheses for the negative exponents.

HTH
Thomas

[Ángel Martin]

From the SandMath manual, TVM$ section:

f ’ (i) = (PMT / i^2 ) * [ (1+i)^(-n) - 1 ] + n * [PMT (1 + ip)/i – FV ] * (1+i)^[-(n+1)]

Cheers,
ÁM.

[Dieter]

(12-02-2014 05:33 AM)iMatt04bit Wrote:  The results appear precise and a match to a variety of real and simulated calculators.

Even with very small values of i ?-)

Some HP calculators use special functions like ln(1+x) or exp(x)-1 to ensure exact results even under these circumstances. Have you implemented something similar in Java? What's the working precision? Standard double precision, i.e. approx. 15 decimal digits?

(12-02-2014 05:33 AM)iMatt04bit Wrote:  I just need to know if I have a correct derivative. ...if only I had taken Calculus...

If you didn't there is always Wolfram Alpha. ;-)

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // formula from the hp12c manual, solves to 0 for matched I converted to R, used below
f(x) = PV + (1+R*K) * PMT * ((1-(1+R)^-N)/R) + FV * (1+R)^-N

I assume you mean f(R), not f(x).

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // the derivative as best I know
f'(x) = PV + (1+R*K) * PMT * [1/N * (1-(1+R)^-(N-1)) / (1/R)] + FV * [1/N * (1+R)^-(N-1)]

Hmmm... not quite. ;-)

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // compute for +I
I = I - f(x)/f'(x)

// or for -I
I = I + f(x)/f'(x)

Could you explain this a bit more in detail?
And what about cases with, say, two different positive solutions?
What is the initial guess for i, what value does the iteration start with?
How do you handle cases where no solution exists?

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // breakpoint on difference of prior computed I and current I is <= 5e-12 for fractional I | I < 1 percent
// or computed I <= 5e-15 for I | I >= 1 percent

Sorry, I do not quite understand what you want to say here.

The usual exit condition for the iteration loop is a relative (!) error near the working precision, i.e. here 1E-15. Due to numeric pitfalls this may not occur under some circumstances, so I often use another idea: Since approach Newton's method converges roughly quadratically, the idea is to quit as soon as a relative error as small as, say, 1E-10 is reached. Assuming quadratic convergence, this means that the next iteration would have an error as low as 1E-20 and thus the approximation of i would not change — so you're done and the result is fine. Just an idea...

Finally, here are two links that may be helpful:

First you may take a look at the way the TVM equation is solved in the HP41 standard application pac. The code can be found in the HP41 software section. Hint: R07 holds 1+i%/100 and R09 stores i%/100.

And finally, as you might have expected, the TVM issue has been discussed numerous time before. Maybe a look at this thread in the old forum may provide some useful information.

Dieter

[iMatt04bit]

(12-03-2014 09:17 PM)Dieter Wrote:  [quote='iMatt04bit' pid='22217' dateline='1417498436']
The results appear precise and a match to a variety of real and simulated calculators.

Quote:Even with very small values of i ?-)

It works with an interest value entry of 0.005/12, and up to 999/12.
The very small value is not solving as well as I want, though still accurate, but normal values do solve really well. I still get results that match the hp12c and hp20b for all values, so I feel okay with the code so far. It should solve to the maximum of display characters as input - still have to test it out to see how well that works.

I still have more scenarios to run through the code to be certain that it is covering all the conditions. You should see the spreadsheet that I have with all the TVM scenarios and much of where I worked out the i guesser code formulas.

Quote:Some HP calculators use special functions like ln(1+x) or exp(x)-1 to ensure exact results even under these circumstances. Have you implemented something similar in Java? What's the working precision? Standard double precision, i.e. approx. 15 decimal digits?

The precision is exact in all tested cases to the hp12c emulator I am using, and 1e-10 almost always vs the hp20b emulator, except a rare occasional 1e-9 difference.
Java seems to be producing 1e-15 precision, but I really don't have a method to test it fully. I do have logarithms, Eulers, all kinds of good intrinsic math formulas in Java. I am very interested in applying the two functions you mention.

(12-02-2014 05:33 AM)iMatt04bit Wrote:  I just need to know if I have a correct derivative. ...if only I had taken Calculus...

Quote:If you didn't there is always Wolfram Alpha. ;-)

I checked that out (thanks!), its really nice site.
I also found one that solves and gives a tutorial at calculus-calculator.com

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // formula from the hp12c manual, solves to 0 for matched I converted to R, used below
f(x) = PV + (1+R*K) * PMT * ((1-(1+R)^-N)/R) + FV * (1+R)^-N

Quote:I assume you mean f(R), not f(x).

Yes, as in R being the mutating value of i/100, fed into the formulas for fx and fx'

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // the derivative as best I know
f'(x) = PV + (1+R*K) * PMT * [1/N * (1-(1+R)^-(N-1)) / (1/R)] + FV * [1/N * (1+R)^-(N-1)]

Hmmm... not quite. ;-)

I found the real one and its working well.
Its really odd how that version and similar others were producing reasonable results.
At times I think it had better results than the real thing.
But it seemed to be confined to particular ranges and I wanted more.

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // compute for +I
I = I - f(x)/f'(x)

// or for -I
I = I + f(x)/f'(x)

Quote:Could you explain this a bit more in detail?
And what about cases with, say, two different positive solutions?
What is the initial guess for i, what value does the iteration start with?
How do you handle cases where no solution exists?

I have since moved to just using one of those formulas.
I was having a bad sign day...

These functions are inside a loop to form the complete function of the Newton method.

The fx function is the HP TVM formula that solves to zero, is in a variable named fxA.
The fx' function which has this massively long derivative formula, is in a variable named fxB.

Both functions receive the input of R, which at first is the seed value (aka guess of i), in the form R = x/100, where x is the i value in the loop.

Next these function results are computed into the evolving version of i, in variable x.
So next: x = x - fxA/fxB;

When the loop cycles again the x value feeds into R, which in turn feeds into fxA,fxB.
As in: R = x/100;

I really hope this is the right way, it seemed so from examples I have reviewed.
I still think it takes far too many iterations to solve.
The partial derivatives that I was using solved much faster.
But still the same format of x = x - fxA/fxB was in use.
Do you know how fast this should solve?
A rough estimate of expected iterations?

For an initial value/guess of i, there is a set of conditional code statements which calculates an estimate depending on variations of the TVM inputs of PV,PMT,FV and if all are non-zero, or if one of them is zero, and also per signs. The estimates are within 0.5 of the solution at times, and aimed slightly above the expected value of i.
I also use a rough estimater that simply walks up the powers of 10 from 1 to 9e+10 and tests for a good match using the formula assigned to variable fxA. This involves multiple stepped loops which do not iterate through every digit until refining the value.

I am solving both pos/neg double roots - found an example in the WP34S archives.
I look for a condition where a double root is expected, like (PV > 0.0) && (PMT < 0.0) && (FV > 0.0), and then I compare multiple guess values, where if there are different signs, then both pos/neg values are queued for the Newton solver. It also catches conditions where a lower value of i will also solve. So both types of roots are covered.

When no solution is found, a value of zero is returned.
Best idea for an error code that I have at the moment since i != 0 is expected.

(12-02-2014 05:33 AM)iMatt04bit Wrote:  // breakpoint on difference of prior computed I and current I is <= 5e-12 for fractional I | I < 1 percent
// or computed I <= 5e-15 for I | I >= 1 percent

Quote:Sorry, I do not quite understand what you want to say here.

I save the prior value of computed interest, payment, and the result of the TVM formula in fxA. The differences of these values are compared on each iteration. When little or no difference it is time to break out of the loop.

I print every iteration to get a visual on how well the solver is working.
The values adf, idf, pdf : are differences of prior and current values.
adf : prior fxA and current fxA
idf : prior/cur interest
pdf : prior/cur pmt
dfPMT : computed pmt vs actual PMT
x : permutating interest
fxA : TVM formula which should resolve to 0.00 under ideal conditions
ct : certain near exit conditions engage a counter as a backup loop exit method

Code:


Using: N: 10  PV: 50.00  PMT: -30.00  FV: 400.00
Using i guess of: 60.200000000000000
   1)  fxA: 42.065379721231366  fxB: 14.581248976719598  x: 57.315104488758620  pmt: -32.555287548788040  PMT: -30.00
         ct:  0   adf: 55.796304578028526  idf: 2.884895511241382  pdf: 32.555287548788040  dfPMT: 2.555287548788037
   2)  fxA: 25.304006313140775  fxB: 14.691976764111903  x: 55.592803486013310  pmt: -31.466094515010040  PMT: -30.00
         ct:  0   adf: 16.761373408090590  idf: 1.722301002745311  pdf: 1.089193033777995  dfPMT: 1.466094515010042
   3)  fxA: 14.953452283311517  fxB: 14.723078734743643  x: 54.577156381714110  pmt: -30.841422907263926  PMT: -30.00
         ct:  0   adf: 10.350554029829258  idf: 1.015647104299198  pdf: 0.624671607746116  dfPMT: 0.841422907263926
// snipped
  56)  fxA: -0.000000000001847  fxB: -41.025090567119340  x: 14.435871328079948  pmt: -29.999999999999964  PMT: -30.00
         ct:  5   adf: 0.000000000002274  idf: 0.000000000000044  pdf: 0.000000000000050  dfPMT: 0.000000000000036
  57)  fxA: -0.000000000000426  fxB: -41.025090567119560  x: 14.435871328079937  pmt: -29.999999999999993  PMT: -30.00
         ct:  5   adf: 0.000000000001421  idf: 0.000000000000011  pdf: 0.000000000000028  dfPMT: 0.000000000000007
Break countdown (5) on abs(fxA) < 1e-12 : -0.000000000000426
I is: 53.172213268384730,14.435871328079937

Quote:The usual exit condition for the iteration loop is a relative (!) error near the working precision, i.e. here 1E-15. Due to numeric pitfalls this may not occur under some circumstances, so I often use another idea: Since approach Newton's method converges roughly quadratically, the idea is to quit as soon as a relative error as small as, say, 1E-10 is reached. Assuming quadratic convergence, this means that the next iteration would have an error as low as 1E-20 and thus the approximation of i would not change — so you're done and the result is fine. Just an idea...

Quote:Finally, here are two links that may be helpful:

I have a set of variables where the difference of the prior value and current value is compared. This is done for computed values of pmt, i, fxA, and a difference of pmt vs PMT, where pmt is computed in the loop using the permutating i, and PMT is the TVM value received by the solver function. I have been able to use fairly extreme cutoffs of 1e-15 for many, but have some backup solutions for where they may fail. Still testing things, but I can break out on the solve of 0.005/12, so they seem to be ok.

Quote:First you may take a look at the way the TVM equation is solved in the HP41 standard application pac. The code can be found in the HP41 software section. Hint: R07 holds 1+i%/100 and R09 stores i%/100.

I'll have to look at that to see how it is used, it looks like good reading.

Quote:And finally, as you might have expected, the TVM issue has been discussed numerous time before. Maybe a look at this thread in the old forum may provide some useful information.

I think I have been everywhere and now have a nice collection of documents and manuals. Did a lot of reading online about algorithms for solving for root. There is not enough info regarding 5th degree polynomials though. ; (

Quote:Dieter

Thanks for the info, much appreciated, Matt

[Dieter]

(12-11-2014 10:02 AM)iMatt04bit Wrote:  

(12-03-2014 09:17 PM)Dieter Wrote:  Even with very small values of i ?-)

It works with an interest value entry of 0.005/12, and up to 999/12.

That's not "small". ;-) Try something like 1E-14/12. Or anything that rounds to 1 if you evaluate 1+R.

Quote:I am very interested in applying the two functions you mention.

These functions, ln1+x and e^x-1 have been discussed several times here. The classic method by W. Kahan for ln1+x can be found in the HP15C advanced functions handbook, others have been discussed in this forum, for instance here.

Quote:I found the real one and its working well.
Its really odd how that version and similar others were producing reasonable results. At times I think it had better results than the real thing.

If you consider the way Newton's method works, it's well possible that a slightly wrong derivative may work better (or worse) than the correct one.

Quote:I really hope this is the right way, it seemed so from examples I have reviewed.
I still think it takes far too many iterations to solve.
The partial derivatives that I was using solved much faster.
But still the same format of x = x - fxA/fxB was in use.
Do you know how fast this should solve?
A rough estimate of expected iterations?

Let's take a look at your example.

Quote:For an initial value/guess of i, there is a set of conditional code statements which calculates an estimate depending on variations of the TVM inputs of PV,PMT,FV and if all are non-zero, or if one of them is zero, and also per signs. The estimates are within 0.5 of the solution at times, and aimed slightly above the expected value of i.

Does this mean the estimate is off by 50% or by 0.5%?
The latter would be sensational. ;-)

Quote:I save the prior value of computed interest, payment, and the result of the TVM formula in fxA. The differences of these values are compared on each iteration. When little or no difference it is time to break out of the loop.

Please see my previous post on this question.

Now let's take a look at your example case.
I have removed most of the data and left only x, the current approximation for the interest rate.

Code:


Using: N: 10  PV: 50.00  PMT: -30.00  FV: 400.00

Using i guess of: 60.200000000000000

   1)  x: 57.315104488758620
   2)  x: 55.592803486013310
   3)  x: 54.577156381714110 
// snipped
  56)  x: 14.435871328079948
  57)  x: 14.435871328079937
Break countdown (5) on abs(fxA) < 1e-12 : -0.000000000000426
I is:
   53.172213268384730%
   14.435871328079937%

WHAT? FIFTY-SEVEN ITERATIONS? #-)
You must be kidding. There's definitely something going wrong here.

I tried the quick-and-dirty spreadsheet I once used to estimate an initial guess for the interest rate. My first guesses for this case were 9,836% and 60%. Starting with these, the iteration converges as follows. The third column shows the results for your estimate 60,2% which works equally well.

Code:


 0)    9,836 06557 37705      60,000 00000 00000      60,200 00000 00000
 1)   13,745 64880 79411      53,563 10589 39787      53,582 35595 39331
 2)   14,416 85886 71670      53,174 01364 11121      53,174 19311 90885
 3)   14,435 85632 97618      53,172 21330 74545      53,172 21331 56320
 4)   14,435 87132 80707      53,172 21326 83847      53,172 21326 83847
 5)   14,435 87132 80799      53,172 21326 83847      53,172 21326 83847
 6)   14,435 87132 80799      53,172 21326 83847      53,172 21326 83847

As you can see, the approximation converges after five resp. four (!) iterations.

So it looks like there is something wrong with your iteration loop. Your initial guess of 60,2% is fine and converges within four iterations as well, so there must be something wrong within the loop.

Convergence should develop roughly quadratic, and this is true for the example above. Take a look at the right column starting at 60,2%. The first approximation has two correct digits, the next one four, then eight and finally all 15 digits are correct. The left column shows the same pattern: the first approximation has 1–2 valid digits, and then three, six, twelve and finally all digits are correct.

You may check whether the derivative is really correct by comparing its results with a numeric approximation: calculate f(r) and f(r+0,0001) and divide the difference of these two by 0,0001.

Addendum:
It looks like both your f(r) and f'(r) are not correct. Starting with 60,2%, the TVM equation (fxA) should return a value of 4,2065 and not 42,065 as stated in your results. Also the first derivative at this point (fxB) is roughly 56,81 and not 14,581. The correct values lead to a first approximation of 0,602 – 4,2065/56,81 = 0,52795... or 52,8%.

The second iteration then should give fxA = –0,234 and fxB = 62,97 which leads to a new x of 53,17% – which is already exact to four significant digits.

Dieter

[Dieter]

(12-11-2014 10:02 AM)iMatt04bit Wrote:  // formula from the hp12c manual, solves to 0 for matched I converted to R, used below
f(x) = PV + (1+R*K) * PMT * ((1-(1+R)^-N)/R) + FV * (1+R)^-N
...
// the derivative as best I know
...

Since there seems to be something wrong with your iterative calculation, here's a first suggestion on how to handle the TVM function, its derivative and the iterative approximation. You may want to try it and see what you get.

Here r is the (decimal) interest rate (e.g. 0.05 for 5%) and k is either 0 or 1, depending on END or BEGIN mode.

Since the equations will try to divide by zero if r=0, it might be a good idea to provide a first guess for r which equals the approximation that would follow an estimate of r=0. I used the following formula:

Code:

  r = 2 * (pv + n * pmt + fv)  /  ((fv * (2*k-1 + n) + pv * (2*k-1 - n)))

– If r = 0 we're already done: r=0 is the exact solution.
– If r > 0 we have a first approximation that will be improved during the following iteration loop.
– If r < 0 using half that (r := r/2) as the initial guess seems to work better, especially for large n.

Now the iteration may start. The TVM function and its derivate may be evaluated quite efficiently this way:

Code:

  q = 1 + r
  a = q ^ (-n)
  b = pmt / r * (1 - a)

  TVMfunction   =  pv + b * (1 + k * r) + fv * a
  TVMderivative =  n * a / q * (pmt * (1 / r + k) - fv) - b / r

Finally r is adjusted:

Code:

  correction = TVMfunction / TVMderivative
  r = r - correction

Repeat this until abs(correction) is sufficiently small, e.g. ≤ abs(1E-12 * r), or the max. number of iterations is ecceded (which means no convergence and should throw an error).

For your example case this yiels the following results:

Code:

 #  Approximation for r    TVM(r)   /    TVM'(r)    =   correction
-------------------------------------------------------------------
 0  0,09836 06557 37705    20,8922      -624,4903      -3,3455 E-2
 1  0,13181 54769 78560     4,3458      -381,4222      -1,1394 E-2
 2  0,14320 92370 54887     0,36344     -319,0893      -1,1390 E-3
 3  0,14434 82151 43331     0,003289    -313,3273      -1,0497 E-5
 4  0,14435 87123 97370     2,7676 E-7  -313,2746      -8,8343 E-10
 5  0,14435 87132 80800    -9,948 E-14  -313,2746       3,1754 E-16
 6  0,14435 87132 80799

BTW the TVM equation may be written in numerous different ways. Some can be evaluated more easily, some less so, some provide faster convergence, others converge slower. So the way the TVM equation is shown in the 12C manual may not be the final word. ;-)

Dieter

[iMatt04bit]

I found the problem in the solver code and its resolved now.
There was a misplaced sign for the derivative and the rate needed a correction.
I am using R in the formulas, where R = x/100, so after x = (fxA/fxB * 100) it works well now with results matching up exactly to the hp12c emu and usually exact or 1e-10 difference to the hp20b emu. I tested it with all the TVM scenarios in the hp12c manual, where I turned the scenarios to solve for the interest rate.

I am using the following formulas, where x is the interest rate in the form displayed on the calculator. Where fxA is the TVM solution to zero, and fxB the derivative.

Code:

R = x/100;
fxA = PV + (1+R*K) * PMT * ((1-Math.pow(1+R,-N))/R) + FV * Math.pow(1+R,-N);
fxB = (1/Math.pow(R,2)) * Math.pow(1+R,(-N-1)) * (-FV * N * Math.pow(R,2) + K * N * PMT * Math.pow(R,2) + N * PMT * R - PMT * R * Math.pow(1+R,N) + PMT * R - PMT * Math.pow(1+R,N) + PMT);

I feel like the Newton method loses resolution for extremely small numbers, like 0.005/100 for an interest rate or smaller. And it has a breaking point on the larger scale (more than 99999999) for the interest rate (can be solved on the hp20b).
The reason for the problem with large numbers is the size of the value in fxB.

So for an alternative I am experimenting with the Secant method (I have it working). Wiki is such a great resource.
It should work well for extremely large numbers since it only requires the TVM solution formula (fxA) and 2 points, which are easy, just use p1=guess + guess/2 and p2=guess - guess/2.

Before the loop:

Code:

x0 = i + i/2;
x1 = i - i/2;

This is the code for the Secant method (in the loop):

Code:

R = x0/100;
fx0 = PV + (1+R*K) * PMT * ((1-Math.pow(1+R,-N))/R) + FV * Math.pow(1+R,-N);
R = x1/100;
fx1 = PV + (1+R*K) * PMT * ((1-Math.pow(1+R,-N))/R) + FV * Math.pow(1+R,-N);

x1 = x1 - fx1 * ((x1 - x0)/(fx1 - fx0));

There are an impressive amount of root finder methods out there, so its hard to know what others may prove useful for an alternative. I think a primary concern is the difficulty in solving small numbers that approach zero.

I want to thank everyone for the help - its much appreciated.

[Dieter]

(12-17-2014 01:56 AM)iMatt04bit Wrote:  I am using the following formulas, where x is the interest rate in the form displayed on the calculator.

I have not checked this yet.
But how fast does it converge (# iterations) and which initial guess do you use?

(12-17-2014 01:56 AM)iMatt04bit Wrote:  I feel like the Newton method loses resolution for extremely small numbers, like 0.005/100 for an interest rate or smaller.

That's not the fault of Newton's method. Expressions like (1+i)ⁿ – 1 cause severe roundoff errors if implemented the way you did, i.e. with standard math functions like pow or exp.

(12-17-2014 01:56 AM)iMatt04bit Wrote:  So for an alternative I am experimenting with the Secant method (I have it working). Wiki is such a great resource.

Methods like the secant method or Regula Falsi work fine if there are two initial approximations that bracket the true root with different signs of the function to solve, i.e. f(x1) * f(x2) < 0. How do you make sure this condition is met, and what are the two guesses you start with?

(12-17-2014 01:56 AM)iMatt04bit Wrote:  There are an impressive amount of root finder methods out there, so its hard to know what others may prove useful for an alternative. I think a primary concern is the difficulty in solving small numbers that approach zero.

Again, that problem that can be addressed with a careful implementation. I will give an example. For the moment let's assume 15 digits working precision.

Code:


  let   i  =  1E-10%
        n  =  10 years with 12 payments/yr
  then
              r  =  8,3333 33333 33333 E-14
            1+r  =  1,0000 00000 00008
   (1+r)^120     =  1,0000 00000 00960
   (1+r)^120 - 1 =  9,6000 00000 00000 E-12
 
compare this with the true result:
 
   (1+r)^120 - 1 =  1,0000 00000 00496 E-11
 

Do the same calculation with 12 digits and 1+r will even numerically equal 1, leading to a plain zero.

These problems can be avoided by using the already mentioned special functions ln(1+x) and e^x–1. These are also available in Java as log1p and expm1. So you better do the above calculation this way:

Code:


              r  =  8,3333 33333 33333 E-14
       log1p(r)  =  8,3333 33333 33298 E-14
       120 * "   =  9,9999 99999 99958 E-12
       expm1(")  =  1,0000 00000 00496 E-11
 

Of course this implementation is strongly recommended for any solving method.

Dieter

Edited to correct some typos

[iMatt04bit]

I found an improvement for the Secant method where it will solve much faster:

Code:

x0 = x0 - fx0 * ((x0 - x1)/(fx0 - fx1));

This should follow the other line with similar code.
Now the secant method decrements from both x0 and x1 points as it solves.

I have run some preliminary tests using a variety of TVM scenarios with good results.
Min/max number of iterations to solve for 21 scenarios tested are:
Newton: 2,6 and Secant: 5,9
The Newton method appears to offer the best accuracy.

I left out the stats for the double root scenarios because I am still looking for a formulaic method for determining:
1) a double root exists 2) an approximation of the 2nd root value - with an indication of sign.

I used the Euler methods (thanks Dieter) which have shown improved results.
A code sample:

Code:

f1 = Math.exp(N * Math.log1p(R));  // (1+R)^N
f2 = Math.exp(-N * Math.log1p(R));         // (1+R)^-N
f3 = Math.exp( (-N-1) * Math.log1p(R));  // (1+R)^(-N-1)

I attached 2 files for anyone interested.
The pdf file shows a comparison of values between emu-calculators and both solver methods. And the text file shows the output giving an indication of how well the methods are solving.

[Dieter]

(12-20-2014 05:47 AM)iMatt04bit Wrote:  I found an improvement for the Secant method where it will solve much faster:

Code:

x0 = x0 - fx0 * ((x0 - x1)/(fx0 - fx1));

This should follow the other line with similar code.

Hmm... looks like the usual formula for the Secant method to me. ;-)

(12-20-2014 05:47 AM)iMatt04bit Wrote:  Now the secant method decrements from both x0 and x1 points as it solves.

Determining if (and how fast) the iteration converges is a quite tricky subject. Using different forms of the TVM equation may drastically change this behaviour, so the version given in the 12C manual is not neccessarily the best choise. I also would reccommend a look at the Kahan paper linked below.

(12-20-2014 05:47 AM)iMatt04bit Wrote:  Min/max number of iterations to solve for 21 scenarios tested are:
Newton: 2,6 and Secant: 5,9

So Newton's method converges faster?

Looking at "results.txt" the initial guesses look very good. How do you evaluate these?

Here's a paper by W. Kahan (who contributed much of the mathematics used in HP's calculators since the mid-Seventies and most of the later iEEE754 floating point standard): Mathematics written in sand addresses several pitfalls and problems using floating point arithmetics on computers and calculators. Especially pages 14 ff. may be interesting here as they deal with financial calculators and problems at solving for the interest rate. The "A Penny for your Thoughts" example also demonstrates why using log1p instead of pow is a good idea. You may also want to read p. 17–22 concerning the convergence of the Newton resp. Secant method.

Dieter

[iMatt04bit]

(12-21-2014 06:14 PM)Dieter Wrote:  

(12-20-2014 05:47 AM)iMatt04bit Wrote:  I found an improvement for the Secant method where it will solve much faster:

Code:

x0 = x0 - fx0 * ((x0 - x1)/(fx0 - fx1));

This should follow the other line with similar code.

Quote:Hmm... looks like the usual formula for the Secant method to me. ;-)

Wiki shows something like this:
x2 = x1 - f(x1) * [ (x1-x0)/(f(x1)-f(x0) ]
x3 = x2 - f(x2) * [ (x2-x1)/(f(x2)-f(x1) ]
x4 ...etc
It looks like single variable substitution returned to the input of the next iterative function. But maybe I have been staring at it for too long. Before I tried a dual solving method - the Wiki version of the function set only converged from one endpoint.

This is the code I used where I added a second convergence process as x0 = formula:

Code:

// initial guess converted to 2 points expected to provide a reasonable closed set
x0 = iguess + iguess/2;
x1 = iguess - iguess/2;

for(;;){

  // the TVM formula from the perspective of each point
  R = x0/100;
  fx0 = PV + (1+R*K) * PMT * ((1-Math.pow(1+R,-N))/R) + FV * Math.pow(1+R,-N);
  R = x1/100;
  fx1 = PV + (1+R*K) * PMT * ((1-Math.pow(1+R,-N))/R) + FV * Math.pow(1+R,-N);

  // double convergence via dual variable substitution
  x1 = x1 - fx1 * ((x1 - x0)/(fx1 - fx0));
  x0 = x0 - fx0 * ((x0 - x1)/(fx0 - fx1));

  // breakpoint code for when a solution is reached goes here
  // something like { abs(x0 - x1) < 1e-15 } and/or other methods
}

(12-20-2014 05:47 AM)iMatt04bit Wrote:  Now the secant method decrements from both x0 and x1 points as it solves.

Quote:Determining if (and how fast) the iteration converges is a quite tricky subject. Using different forms of the TVM equation may drastically change this behaviour, so the version given in the 12C manual is not neccessarily the best choise. I also would reccommend a look at the Kahan paper linked below.

I reformed the formulas to make use of the logarithmic functions and to provide for large scale numbers. Note that in the code snippet I am showing the original format of the formulas for clarity.

Which brings a question to mind: is there a min/max limit to the size of the interest rate value that the ln/log function can handle?

(12-20-2014 05:47 AM)iMatt04bit Wrote:  Min/max number of iterations to solve for 21 scenarios tested are:
Newton: 2,6 and Secant: 5,9

Quote:So Newton's method converges faster?

Yes, it also looks a bit more accurate at times, by the appearance of the fxA value vs the fx1 value.

Quote:Looking at "results.txt" the initial guesses look very good. How do you evaluate these?

I have a set of simple calculations + the one you showed me is in there.
There are 10 calculations, and depending on the TVM values some will provide a reasonable guess under certain conditions. I evaluate each result using the TVM formula, and then cheat a little by running the best one through a single pass of the Newton formula.

A last evaluation of the guess before and after the Newton is evaluated and used as the final estimate of i.

I still am experimenting with these estimating formulas.

Quote:Here's a paper by W. Kahan (who contributed much of the mathematics used in HP's calculators since the mid-Seventies and most of the later iEEE754 floating point standard): Mathematics written in sand addresses several pitfalls and problems using floating point arithmetics on computers and calculators. Especially pages 14 ff. may be interesting here as they deal with financial calculators and problems at solving for the interest rate. The "A Penny for your Thoughts" example also demonstrates why using log1p instead of pow is a good idea. You may also want to read p. 17–22 concerning the convergence of the Newton resp. Secant method.

I am already reading it - thanks very much!
This is the stuff that makes math an adventure.

Dieter

[iMatt04bit]

Ran a test using the scenarios listed on page 17 of MathSand and results look good.
Using the full compliment of logarithms.
I cut the number of pre-guess selections down to 6 of the best ones.
Different TVM scenarios will call for the best choice, which is usually the Dieter formula.

Dieter I have to ask: how did you determine that formula? Its awesome.
2 * (PV + N * PMT + FV) / (FV * (2*K-1.0+N) + PV * (2*K-1.0-N))

I'll post results.

-matt

[Dieter]

(12-24-2014 05:17 AM)iMatt04bit Wrote:  Dieter I have to ask: how did you determine that formula? Its awesome.

It's trivial. ;-)

That's simply the next approximation you get with Newton's method if you start with r=0, using the TVM formula as stated in the documentation for the WP34s TVM solver. So it's nothing but

r   ≈   0 – TVM(0) / TVM'(0)

You can also do this with the TVM formula as it is given in the 12C manual. In this case the result is

r   ≈   –2 * (pv + n*pmt + fv) / (n * (pmt*(n+2*k-1) + 2*pv))

In other words: you could just as well start with r=0 and you'd get these approximations after the first iteration step. ;-)

However, this would require special handling of the (removable) singularity at r=0 where the TVM function and its derivative are not defined (both approach a limit as r→0, and these limits are used for the first guess). If the mentioned approximation is used as the first guess and the next iteration step still yields r=0, the solution is r=0. Else we can expect the following iterations to lead away from r=0, thus avoiding problems with the singularity at that point.

In many cases these two approximations bracket the true interest rate, or their mean is an even better estimate that comes quite close. Using your original example (n=10, PV=50, PMT=–30, FV=400, k=0) the two approximations are 9,836% and 17,647% (mean = 13,74%) while the exact rate is 14,436%.

However, this (one approximation a bit low, the other one a bit high) is not the case under all circumstances. ;-)

Addendum (Dec 26): Obviously we ain't seen nothing yet. I just read HP Journal 10/1977 on the new HP92 financial calculator. The article mentioned the "new" TVM formula (with signed parameters) and the problems when solving for the interest rate. Most of these problems look quite familar, and the author stresses the importance of a decent initial guess. Quote: "...a strategy was developed that produces inital guesses accurate to five decimal places". Well, that's awesome. :-)

BTW there also is an interesting read on the way the good old HP80 tries to solve for the interest rate (at least in a standard case like FV=0) which is described in HP Journal 5/1973 (cf. Appendix A following the article). Here the author even describes how the formula for the first guess was developed. Please note that at that time all variables were entered unsigned, i.e. usually positive.

We now can see whether we can develop a strategy that produces decent guesses for various scenarios. We can tell quite easily whether a solution exists and if there are multiple or no solutions at all. Upper limits for the interest rate can be given for some standard cases, e.g. FV=0 or PV=0, where the interest rate approaches a limit as n approaches infinity (which helps solve these ciritcal cases for large n where convergence may become extremely slow – which is also mentioned in the HP Journal article). Designing an algorithm that estimates the interest rate for different scenarios will be a tricky task, but I think it can be done. However, compared to HP's, I assume my skills are too limited to get a result with five valid digits in the initial guess. On the other hand, for large n and FV=0 or PV=0 this can be accomplished quite easily.

Your turn now. :-)

Dieter

[iMatt04bit]

(12-25-2014 12:42 PM)Dieter Wrote:  

(12-24-2014 05:17 AM)iMatt04bit Wrote:  Dieter I have to ask: how did you determine that formula? Its awesome.

Quote:It's trivial. ;-)

Sure, for you... I was staring at it in amazement. :0

That's simply the next approximation you get with Newton's method if you start with r=0, using the TVM formula as stated in the documentation for the WP34s TVM solver. So it's nothing but

r   ≈   0 – TVM(0) / TVM'(0)

You can also do this with the TVM formula as it is given in the 12C manual. In this case the result is

r   ≈   –2 * (pv + n*pmt + fv) / (n * (pmt*(n+2*k-1) + 2*pv))

In other words: you could just as well start with r=0 and you'd get these approximations after the first iteration step. ;-)

However, this would require special handling of the (removable) singularity at r=0 where the TVM function and its derivative are not defined (both approach a limit as r→0, and these limits are used for the first guess). If the mentioned approximation is used as the first guess and the next iteration step still yields r=0, the solution is r=0. Else we can expect the following iterations to lead away from r=0, thus avoiding problems with the singularity at that point.

In many cases these two approximations bracket the true interest rate, or their mean is an even better estimate that comes quite close. Using your original example (n=10, PV=50, PMT=–30, FV=400, k=0) the two approximations are 9,836% and 17,647% (mean = 13,74%) while the exact rate is 14,436%.

However, this (one approximation a bit low, the other one a bit high) is not the case under all circumstances. ;-)

Quote:Addendum (Dec 26): Obviously we ain't seen nothing yet. I just read HP Journal 10/1977 on the new HP92 financial calculator. The article mentioned the "new" TVM formula (with signed parameters) and the problems when solving for the interest rate. Most of these problems look quite familar, and the author stresses the importance of a decent initial guess. Quote: "...a strategy was developed that produces inital guesses accurate to five decimal places". Well, that's awesome. :-)

I am still trying to figure out how that is done - 5 decimal places on the 1st guess!
Sure sounds like rocket science to me. ; )

Quote:BTW there also is an interesting read on the way the good old HP80 tries to solve for the interest rate (at least in a standard case like FV=0) which is described in HP Journal 5/1973 (cf. Appendix A following the article). Here the author even describes how the formula for the first guess was developed. Please note that at that time all variables were entered unsigned, i.e. usually positive.

I am so impressed with the minds behind those HP innovations.
Maybe I can deduce how to plot some of this... my graphing skills are a bit rusty.

We now can see whether we can develop a strategy that produces decent guesses for various scenarios. We can tell quite easily whether a solution exists and if there are multiple or no solutions at all. Upper limits for the interest rate can be given for some standard cases, e.g. FV=0 or PV=0, where the interest rate approaches a limit as n approaches infinity (which helps solve these ciritcal cases for large n where convergence may become extremely slow – which is also mentioned in the HP Journal article). Designing an algorithm that estimates the interest rate for different scenarios will be a tricky task, but I think it can be done. However, compared to HP's, I assume my skills are too limited to get a result with five valid digits in the initial guess. On the other hand, for large n and FV=0 or PV=0 this can be accomplished quite easily.

Your turn now. :-)

Dieter

Thanks for the info, I am reading the journal now.

I revisited the 2nd root prediction and have a solution that seems to be working.
I can get a first approximation in 4-5 iterations, and then feed that to the Newton solver for an exact value.

Essentially I scan in 10ths of the first root between root-1/10th and 0 + 1/10.
I scan to convergence from top-down, bot-up; and for both pos and neg values simultaneously.
If there is either a pos or neg solution, then the loop breaks on convergence.

I am using the TVM formula to evaluate for the best solution (within the loop).

I am working on the hypothesis that the absolute value of the 2nd root is less than the absolute value of the primary root (and greater than zero). Does this seem sane?

This method solves 2 scenarios (the only ones I have at the moment) which have double roots:
N:10 PV:50 PMT:-30 FV:100 with 58.203829688346610, -28.443599888025595
N:10 PV:50 PMT:-30 FV:400 with 53.172213268384720, 14.435871328079967

-matt

[Dieter]

(12-28-2014 06:46 AM)iMatt04bit Wrote:  I am still trying to figure out how that is done - 5 decimal places on the 1st guess!
Sure sounds like rocket science to me. ; )

I can hardly believe it either. Ok, for very large n a limit may exist which can agree with the true root to five or more decimals. But otherwise...

Example – try a standard case: let either FV=0 or PV=0, use a sufficiently large number of periods. Then see how close you get with an estimate of PMT/(FV–PV). ;-)

(12-28-2014 06:46 AM)iMatt04bit Wrote:  I revisited the 2nd root prediction and have a solution that seems to be working.
I can get a first approximation in 4-5 iterations, and then feed that to the Newton solver for an exact value.

Essentially I scan in 10ths of the first root between root-1/10th and 0 + 1/10.
I scan to convergence from top-down, bot-up; and for both pos and neg values simultaneously.

Sounds a bit complicated. ;-)

(12-28-2014 06:46 AM)iMatt04bit Wrote:  I am working on the hypothesis that the absolute value of the 2nd root is less than the absolute value of the primary root (and greater than zero). Does this seem sane?

If two solutions exist, which is the primary root and which is the secondary? In other words: what makes a root primary or secondary? If two roots exist, do they both make sense economically/financially?

(12-28-2014 06:46 AM)iMatt04bit Wrote:  This method solves 2 scenarios (the only ones I have at the moment) which have double roots:

Simply let FV and PV have the same sign (both values positive or negative) and then try different values for PMT. Depending on the latter, either two or no solutions exist.

(12-28-2014 06:46 AM)iMatt04bit Wrote:  N:10 PV:50 PMT:-30 FV:100 with 58.203829688346610, -28.443599888025595
N:10 PV:50 PMT:-30 FV:400 with 53.172213268384720, 14.435871328079967

FTR: the exact 30-digit values are
  –28,443 59988 80255 92557 23584 03219 %
    58,203 82968 83466 11648 87241 82543 %
resp.
    14,435 87132 80799 56974 20470 12989 %
    53,172 21326 83847 24310 46024 13744 %

If your Java version uses standard binary double precision arithmetics you should not trust more than 15 digits. At most. ;-)

Dieter