(10C) Mortgage Loan Interest Rate

[Gamo]

This program was adapted from HP-11C Solutions Handbook on page 152

This program will calculate the "interest rate" on a loan with equal periodic payments.
The user must specify the number of periods, the present value and the payment amount.

Procedure:

n [ENTER] PV [ENTER] PMT [CHS] [R/S]

Example:

You recently obtained a $2500 car loan for 36 months.
If your monthly payment is $86.67, What is the annual percentage rate?

FIX 2
36 [ENTER]
2500 [ENTER]
86.67 [CHS] [R/S] --> 1.25 %Monthly Rate
12 [x] ---> 15.01 %Annual Rate

Program:

Code:


01  ÷
02 STO 3
03 Rv
04 STO 1
05 RCL 3
06 X^2
07 √x
08 ENTER
09 1/x
10 X<>Y
11 RCL 1
12 X^2
13  ÷
14  -
15 STO 2
16 RCL 3
17 X^2
18 √x
19 RCL 2
20  x
21  1
22 RCL 2
23  1
24  +
25 RCL 1
26 CHS
27 Y^X
28 STO 0
29  -
30  -
31 RCL 1
32 RCL 2
33 1/x
34  1
35  +
36  ÷
37  1
38  +
39 RCL 0
40  x
41  1
42  -
43 RCL 2
44  ÷
45  ÷
46 STO+2
47 X^2
48 √x
49 EEX
50  6
51 CHS
52 X≤Y
53 GTO 16
54 RCL 2
55 EEX
56  2
57  x

Gamo

[Albert Chan]

I had thought of this loan estimate before ...

Let compounding factor C = nR/P, so C = 1 if i = 0

C = ni / (1 - (1+i)^-n)
C = 1 + (n+1)/2 i + (n²-1)/12 i² - (n²-1)/24 i³ + ...

if keep only linear terms, i1 = 2*(C-1) / (n+1)
if keep quadratics too, i2 = (2 i1) / (1 + sqrt(1 + (2 i1)*(n-1)/3))

For P=2500, R=86.67, n=36:
C = nR/P = 1.248048
i1 = 0.013408, or annual rate of 16.1%
i2 = 0.012497, or annual rate of 15.0%

For fast calculation, assume n² ~ n²-1, and drop O(i³) terms:

C = i/2 + ((ni + 3)² + 3) / 12

Ignoring i/2 term, ni = √(12 C - 3) - 3 = 0.460719
i = (ni)/n = 0.012798, or annual rate of 12I = 15.4%

Repeat again with i/2 term, solve for ni, annual rate = 15.0%

[Dieter]

(11-01-2018 03:22 AM)Albert Chan Wrote:  For fast calculation, assume n² ~ n²-1, and drop O(i³) terms:

C = i/2 + ((ni + 3)² + 3) / 12

Ignoring i/2 term, ni = √(12 C - 3) - 3 = 0.460719
i = (ni)/n = 0.012798, or annual rate of 12I = 15.4%

Repeat again with i/2 term, solve for ni, annual rate = 15.0%

Albert, this is another case where you may have a good idea but... even after reading your post several times I simply don't understand what you want to say. I neither understand where you get these formulas from nor which method you want to propose.

Could you please give the complete algorithm (for the general case) and a detailled explanation for calculating the interest rate? If the first estimate is 0,0134808, how do you get the exact value 0,01250454600 from this?

Dieter

[Gamo]

Here is the formula to use for the initial guess for [i]
According to this Handbook.

Gamo

[Albert Chan]

(11-01-2018 03:22 AM)Albert Chan Wrote:  Let compounding factor C = nR/P, so C = 1 if i = 0

C = ni / (1 - (1+i)^-n)
C = 1 + (n+1)/2 i + (n²-1)/12 i² - (n²-1)/24 i³ + ...

if keep only linear terms, i1 = 2*(C-1) / (n+1)
if keep quadratics too, i2 = (2 i1) / (1 + sqrt(1 + (2 i1)*(n-1)/3))

For P=2500, R=86.67, n=36:
C = nR/P = 1.248048
i1 = 0.013408, or annual rate of 16.1%
i2 = 0.012497, or annual rate of 15.0%

Hi, Dieter

The formulas were meant for rough estimation, using "normal" calculators (no X^Y)

All formula were derived from C Taylor expansion at i = 0
The quadratic formula were using more accurate form: 2c / (-b - sqrt(b²-4ac))
where, a= (n²-1)/12, b=(n+1)/2, c=1-C, a non-positive number.

This is how the quick formula is derived:
C = 1 + (n+1)/2 i + (n²-1)/12 i² - (n²-1)/24 i³ + ...
~ 1 + (n+1)/2 i + (ni)² / 12
~ i/2 + (3 + 9 + 6(ni) + (ni)²) / 12
= i/2 + ((ni + 3)² + 3) / 12

To get more accurate interest rate, do interpolation (but X^Y is required):
i1 = 0.013408, C1 = n*i1 / (1 - (1+i1)^-n) = 1.267246
i2 = 0.012497, C2 = n*i2 / (1 - (1+i2)^-n) = 1.247888
Interpolate for C = nR/P = 1.248048, we get i = 0.0125045

[Albert Chan]

Hi, Gamo

Thanks for the very good rate guess (post 4).
Rewrite in terms of C, i0 = (C - 1/C)/n

For C = nR/P = (36)(86.67)/2500 = 1.248048, i0 = 0.012411

Using my quick formula, solve for ni, thus i:

i1 = (√(12*(C - i0/2) - 3) - 3) / n = 0.012498, very good estimate

Doing interpolation with above two rates (to match C), we get i = 0.012504549

Edit: doing some plots of C vs i, discovered a even better interpolation scheme.
Interpolate to match √C instead, the line is almost straight, we get i = 0.012504546

[Dieter]

(11-01-2018 02:39 PM)Albert Chan Wrote:  The formulas were meant for rough estimation, using "normal" calculators (no X^Y)

Who does a mortgage calculation, let alone determine the effective interest rate, with such a calculator?

But anyway, if I understand this correctly you simply determine an estimate for the interest rate. The exact result then is determined by (repeated) interpolation, e.g. with the secant method:

(11-01-2018 02:39 PM)Albert Chan Wrote:  To get more accurate interest rate, do interpolation (but X^Y is required):
i1 = 0.013408, C1 = n*i1 / (1 - (1+i1)^-n) = 1.267246
i2 = 0.012497, C2 = n*i2 / (1 - (1+i2)^-n) = 1.247888
Interpolate for C = nR/P = 1.248048, we get i = 0.0125045

But you don't know how accurate this is unless you do some more iterations. ;-)

Re. the HP method from the posted manual page: I remember this has been discussed before, and for this special case it seems to be a quite good estimate.

Dieter

[Albert Chan]

PHP Code:

def solve_r(i, n):  # r = payment / principle
    try: return i/(1.-(1.+i)**-n)
    except ZeroDivisionError: return 1./n 


PHP Code:

def solve_i(n, r, eps=1e-8, verbal=True, f=solve_r):
    i0 = r - 1.0/(r*n*n); f0 = f(i0, n)  # guess i0
    i1 = i0 + 1.5*(r-f0); f1 = f(i1, n)  # guess i1
    if verbal: print '%.15g =>\t%.15g' % (i0, i1)
    while f0 != f1:
        di = (i0-i1)*(r-f1)/(f0-f1); i0 = i1; i1 += di
        if verbal: print '%+.6e => \t%.15g' % (di, i1)
        if abs(di) < eps: break
        f0 = f1; f1 = f(i1, n)
    return i1 


Here is an interest finder by secant's method.
Below example (accurate to 10 places), only 3 f(i) calls are needed.

Code:

>>> solve_i(n=36, r=86.67/2500)
0.0124110212798083 =>   0.0124934261260045
+1.112498e-05 =>        0.0125045511103353
-5.107705e-09 =>        0.0125045460026299
0.012504546002629922

[Albert Chan]

PHP Code:

def nsolve_i(n, r, eps=1e-8, verbal=True):
    i = r - 1.0/(r*n*n)
    if verbal: print 'guess i0 =>\t\t%.15g' % i
    while 1:
        fi = i/(1.-(1.+i)**-n)      # newton's method
        di = (fi - r) / (fi - n*fi*(fi-i)/(1.+i)) * i
        # di = (fi - (fi*r)**0.5) / (fi - n*fi*(fi-i)/(1.+i)) * (i+i)
        i -= di
        if verbal: print '%+.6e => \t%.15g' % (-di, i)
        if abs(di) < eps: break
    return i 


Tried using Newton's method. It is not as fast, unless f(i) is inlined.
Also, if eps is set too small, i will not converge.

Code:

>>> nsolve_i(n=36, r=86.67/2500)
guess i0 =>             0.0124110212798083
+9.356770e-05 =>        0.0125045889822448
-4.297933e-08 =>        0.012504546002918
-9.079922e-15 =>        0.012504546002909
0.012504546002908958

If newton's method apply to √f(i), i converge faster.
Just replace above di, with the commented version.

Code:

>>> nsolve_i(n=36, r=86.67/2500)  # √f(i) version
guess i0 =>             0.0124110212798083
+9.353061e-05 =>        0.0125045518848863
-5.881977e-09 =>        0.0125045460029088
0.012504546002908835

Correct i = 0x1.99bfbc11210d4P-7 ~ 0.01250454600290888
Although √f(i) version does 1 less iteration, it is more accurate (-26 ULP vs +45 ULP)

Correction: Newton's method is actually faster (with similar error)
I had compared the speed with secant's method using same eps = 1e-8
This is not fair to the Newton version, since the result had a few more good digits.

[Dieter]

(11-03-2018 01:51 PM)Albert Chan Wrote:  Tried using Newton's method. It is not as fast, unless f(i) is inlined.
Also, if eps is set too small, i will not converge.

The TVM solvers that I have programmed over the years (well, decades ;-)) all use the Newton method, and usually it converges within a few iterations. In this regard Newton's and the secant method / regula falsi are quite similar. But Newton's method only requires one initial guess instead of two that preferably bracket the solution. How can you realiably ensure this?

If you want to try a secant method I'd recommend implementing a modified algorithm, for instance the Illinois method. This avoids some problems and improves a reasonably fast convergence.

(11-03-2018 01:51 PM)Albert Chan Wrote:  Although √f(i) version does 1 less iteration, it is more accurate (-26 ULP vs +45 ULP)

I think you cannot judge the different methods based on such minor differences. The whole calculation is not that accurate: you can't get 16 valid digits with 16 digit precision. And the implementation already loses several digits when 1+i is calculated, so the last two digits are more or less random. For example, for i=1/30 and n=36 the expression (1+i)ⁿ evaluated with 16 digits will cause an error of 38 ULP. Just for this single calculation, even when evaluated without any additional errors. That's why better implementations of the TVM formula use special functions like e^x–1 and ln(1+x) to avoid such problems, or at least reduce them significantly.

You also mentioned the problem that successive iterations may oscillate in the last digits so that there is no convergence to a unique 16 digit result. The essential reason for this is that neither f(x) nor f'(x) are evaluated absolutely accurate (see above). Since Newton's method and (approximately) the secant method both converge nearly quadratically close to the true root, I would recommend something like this: in each iteration step check the relative error (!), and then compare this not with 1E–15 but 1E–8. Or maybe 1E–11 for a more conservative approach. If the last relative correction term is below this threshold the next iteration should reach an accuracy level close to machine precision, plus or minus some random deviation in the last few digits that cannot be realiably avoided. You may now exit the iteration or, just to be sure, set a flag and do one more final iteration.

Dieter

[Albert Chan]

(11-03-2018 06:11 PM)Dieter Wrote:  The TVM solvers that I have programmed over the years (well, decades ;-)) all use the Newton method, and usually it converges within a few iterations. In this regard Newton's and the secant method / regula falsi are quite similar. But Newton's method only requires one initial guess instead of two that preferably bracket the solution. How can you realiably ensure this?

I started the code from plot of the f(i) vs i. It is relatively straight.
(unlike post 4 formula, which did the flipped version, with i in the denominator)

With my f setup, f(0+) = 1/n ~ 0, f(1) = 1/(1-2^-n) ~ 1. So, for 0 < i <= 1, f'(i) ~ 1
In other words, (r - f(i)) gap is about the same size as i correction.

Since true i is closer to 0% than 100%, I expect flatter f(i).
I wanted i1 closer to true i than i0, so scaled the gap by 1.5 (to compensate)
For small n, that is not quite enough, but big n have slight over-correction.

Quote:I think you cannot judge the different methods based on such minor differences.
The whole calculation is not that accurate: you can't get 16 valid digits with 16 digit precision.

My mistake. The last iteration reached limit of IEEE double precision.
I should compare both version with first iteration, not the last.

Relative error (vs true i), f(i) version = 3.44e-6, √f(i) version = 4.70e-9.
So, starting from the same guess, √f(i) version is 732 times closer after 1 iteration.

Quote:For example, for i=1/30 and n=36 the expression (1+i)ⁿ evaluated with 16 digits will cause an error of 38 ULP

Your really pick a good example, same issue if done with binary math.

float(31/30) = exact(31/30) + 0.467 ULP, almost a half-way case.

lua> (31/30)^36 -- over-estimated 26 ULP
3.255785675929093
lua> 31^36 / 30^36 -- only off 1 ULP, last 82 should be 815
3.255785675929082

26 ULP sounded bad, but not really: 26 ULP = 26/2^51 ~ 1.15e-14
For 16 decimals setup, 38 ULP = 38/10^15 = 3.8e-14, or 3.3 times worse.

Edit: based on above issue, eps is not recommended to set below 1e-12

[Albert Chan]

PHP Code:

def nsolve_i(n, r, eps=1e-8, verbal=True):
    i = r - 1.0/(r*n*n)         # guess
    prev = r
    while 1:
        fi = i/(1.-(1.+i)**-n)  # newton's method
        di = (fi - r) / (fi - n*fi*(fi-i)/(1+i)) * i
        if verbal: print '%.15g\t%+.5e' % (i, -di)
        curr = abs(di)
        if curr >= prev: break  # bad di
        i = i - di
        if curr < eps: break
        prev = curr
    return i 


Code:

>>> nsolve_i(n=36, r=86.67/2500, eps=0)
0.0124110212798083      +9.35677e-05
0.0125045889822448      -4.29793e-08
0.012504546002918       -9.07992e-15
0.012504546002909       -9.44595e-17
0.0125045460029089      +3.54223e-16
0.012504546002908864

Here is a way to avoid Newton's non-convergence due to rounding error.
Stop when adjustment size is bigger than previously.

The last bad adjustment is thrown away, and return previous iteration.
I learn this trick from Kahan's cubic equation solver paper

Edit: above for illustration only.
With eps = 1e-8, we already reached IEEE double limitation. (i = 0.012504546002909)