TVM solve for interest rate, revisited

[Albert Chan]

From https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=234439

Dieter Wrote:I just had a look at the TVM solver that comes with the 34s library software. In most cases there is no closed form solution for the interest rate, so the iterative solver is used. However, the two required initial guesses are simply -100% and (at least) +1000%. While it is true that the result usually will be somewhere in this interval, I think we can do better ...

To "bracket" rate, it may be better to transform TVM equation, with tanh. abs(tanh(x)) ≤ 1

Solve TVM for rate x, using NPMT=0 form:

f = NPMT/n = x/expm1(n*log1p(x))*(pv+fv) + pv*x + pmt

Let C = x*n / (1-1/K), where K = (1+x)^n
Let Ce = C - n*x/2
Let A = (pv+fv)/2
Let B = (pv−fv)/2

f = (Ce−n*x/2)*(pv+fv)/n + pv*x + pmt  = (2A)*(Ce/n) + B*x + pmt

Ce is even function of n
= (C(n=n) + C(n=-n)) / 2
= (x*n/(1-1/K) + x*(-n)/(1-K)) / 2
= x*n/2 * (1+K)/(1-K)
= x*n/2 * (√K + 1/√K) / (√K - 1/√K)
= x*n/2 / tanh(ln(√K))
= x*n/2 / tanh(n/2*log1p(x))

Let t = tanh(n/2*log1p(x)), put Ce back to NPMT:

f = A*x/t + B*x + pmt

f = 0 → x = -pmt / (A/t + B)

t = -1 .. 1      → x = pmt/fv .. -pmt/pv

This assumed pmt ≠ 0
If pmt = 0, we can directly solve for rate, fv = K*pv

If either fv or pv = 0, we only have 1 "edge", but that is OK,
With x guaranteed 1 solution, we just start from the finite edge.

Comment x ranges assumed its edges have opposite sign.
→ (pv, fv) have opposite sign
→ |A| > |B|
→ Denominator (A/t + B) never hit zero, since t = [-1,1]

[Albert Chan]

Example:

lua> n,pv,pmt,fv = 10,50,-30,100
lua> pmt/fv, -pmt/pv
-0.3      0.6

Note that rate cannot be below -1 (otherwise, log1p(x) = NaN)
Both rates are valid (exclusive) edges.

lua> lo = loan_rate(n,pv,pmt,fv,-0.3)
lua> hi = loan_rate(n,pv,pmt,fv, 0.6)
lua> for i=1,4 do print(i, lo(), hi()) end
1      -0.28463423878569005      0.5821069833315318
2      -0.28443603391367706      0.5820382979602389
3      -0.28443599888025706      0.5820382968834661
4      -0.28443599888025595      0.5820382968834662

loan_rate() use Newton's method to solve NPMT=0, for x
Both true rates are bound within its edges, as expected.

[Albert Chan]

(05-13-2022 09:31 PM)Albert Chan Wrote:  Let t = tanh(n/2*log1p(x)), put Ce back to NPMT:

f = A*x/t + B*x + pmt

f = 0 → x = -pmt / (A/t + B)

We can also solve f=0, for t:

tanh(n/2*log1p(x)) = -A*x/(pmt+B*x)

There is a trivial solution, when x=0. We divide both side by x, to remove it.

XCAS> t1 := tanh(n/2*log1p(x))/x
XCAS> t2 := -A/(pmt+B*x)

We approximate t1 by pade approximation.
To keep things simple, we wanted degree(numerator)=1, degree(denominator)=2.
This allow solving for guess rate x, with simple quadratic.

XCAS> p1 := pade(t1,x,3,2)

(6*n+3*n*x) / (12+12*x+2*x^2+n^2*x^2)

Cross multiply, to get quadratic coefs.

XCAS> P := e2r(numer(p1)*denom(t2) - denom(p1)*numer(t2))

[2*A+3*B*n+A*n^2, 12*A+6*B*n+3*n*pmt, 12*A+6*n*pmt]

Let's try previous post example

XCAS> P2 := P(A=(pv+fv)/2, B=(pv-fv)/2)
XCAS> proot(P2(n=10,pv=50,pmt=-30,fv=100))

[-0.268464164628, 0.485855468976]

This is not much better than using edges for rate guesses.
The example is hard. Most cases, rate approximation is good.

Car lease examples, from Fun math algorithms

proot(P2(n=36,pv=30000,pmt=-550,fv=-15000)) .* 12

[-4.89392664913, 0.0696603112801]      // True APR = 6.96608738330 %

Just as previously done with guess_i(), we avoid using square roots.
(also, we only keep the "small" root)

Code:

function guess_i2(n, pv, pmt, fv)
    pmt, fv = pmt or -1, fv or 0
    local a = (pv+fv)/n
    local b = fv-pv-pmt-2*a
    local c = (a+pmt)/b
    b = ((n*n+2)*a+3*(pv-fv))*c/b   -- I coef = [b/c/3, -2, 4c]
    return c * (b-3)/(b-1.5)        -- pade(4c/(1+sqrt(1-4b/3)),b,2,2)
end

lua> guess_i2(10,50,-30,100) -- the "hard" example
-0.3460122699386503

lua> guess_i2(36,30000,-550,-15000) * 12 -- car lease example
0.06966051503474308

[Albert Chan]

(05-13-2022 09:31 PM)Albert Chan Wrote:  Let t = tanh(n/2*log1p(x)), put Ce back to NPMT:

f = A*x/t + B*x + pmt

f = 0 → x = -pmt / (A/t + B)

t = -1 .. 1      → x = pmt/fv .. -pmt/pv

I made an error here.
If t = -A/B is within ± 1, We will be hit by divide by 0

Previous car lease example:

lua> n,pv,pmt,fv = 36,30000,-550,-15000
lua> -(pv+fv)/(pv-fv) -- within +/- 1
-0.3333333333333333
lua> pmt/fv, -pmt/pv -- t=-1 and t=1 edges, both same sign.
0.03666666666666667      0.018333333333333333

At t=0, x = -pmt / (A/0 + B) = -pmt/∞ = 0

Thus, we have 2 ranges: (x > pmt/fv) || (x < -pmt/pv)

From tanh info, rate gap, pmt/fv .. -pmt/pv, is not possible, *opposite* of OP claim.

[Albert Chan]

(04-11-2022 03:11 PM)Albert Chan Wrote:  This may be how I₀ = 1/P - P/N² comes from

Update: Perhaps formula designed also to match asymptote, when C is big.
When compounding factor C is big, literally all payments goes to paying interest.

C = N/P = I*N/(1-(1+I)^-N) ≈ N*I      → I = 1/P

Rate formula matched this behavior: I₀ = 1/P - P/N^2 = (1 - 1/C^2)/P ≈ 1/P

Another way to show rate edges, without complexity of doing tanh transformation (*)
Above quote assumed FV=0. But, we can fix it ...

I = loan_rate(N, PV, PMT, FV) = loan_rate(N, PV+FV, PMT-FV*I, 0)

From RHS, P = (PV+FV) / -(PMT-FV*I)

I = 1/P
PV*I + FV*I = -PMT + FV*I      → I = -PMT/PV

For asymptote edge I, value of FV does not matter.

To get the other edge, we use time symmetry, travelling "backward in time"
(N,PMT) sign flipped, (PV,FV) get swapped.

I = loan_rate(N, PV, PMT, FV) = loan_rate(-N, FV, -PMT, PV)

The other asymptote edge: I = PMT/FV

--

(*) tanh version, with A=(PV+FV)/2, B=(PV-FV)/2:

I = loan_rate(N, PV, PMT, FV) = loan_rate(N, A, PMT+B*I, A)

Numerically to show both forms equivalent:

XCAS> C := I*N/(1-(1+I)^-N); // compounding factor
XCAS> NPMT := C*PV + C(N=-N)*FV + N*PMT

XCAS> solve(NPMT(N=10, PV=50, PMT=-30, FV=100) = 0, I)      → [-0.28443599888, 0.582038296883]
XCAS> A, B := (50+100)/2, (50-100)/2
XCAS> solve(NPMT(N=10, PV=A, PMT=-30+B*I, FV=A)=0, I)      → [-0.28443599888, 0.582038296883]

[Albert Chan]

Instead of cheap tricks to get good rate guess, why not start from the edge, and iterate with Halley's method ?

For NPMT=0 with 2 rate solutions, this picked the "smaller" edge.
Note: NaN edge moved up front, for rate search early termination.

Code:

function edge_i(n,pv,pmt,fv)
    local a, b = pmt/fv, -pmt/pv
    if abs(a)>abs(b) and b>-1 or b~=b then a,b = b,a end
    return a, b  -- a = small edge or NaN
end

We solve NPMT=0, with PV=0, for rate x

x = loan_rate(n, pv, pmt, fv) = loan_rate(n, 0, pv*x+pmt, pv+fv)

f = (pv+fv)/n * g + pv * x + pmt = 0      // same formula used for Plus42

g = n*x / ((1+x)^n-1)
g'/g = -(g-1+(n-1)*x) / (x+x^2) = -num / den
g''/g = (g+n*x)*(2*(g-1)+(n-1)*x) / (x+x^2)^2 = (num+x+1)*(num+(g-1)) / den^2

Let k = (pv+fv)/n * g

f = k + pv*x + pmt
f' = k*(g'/g) + pv
f'' = k*(g''/g)

Using first derivative to get to second, Halley's correction is cheap.
To make code more robust, we special cased with limits when rate goes 0

Code:

function loan_rate2(n, pv, pmt, fv, i)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n, pv, pmt, fv)
    return function()
        if i==0 then
            local a, b = (pv+fv)/n, pv-fv
            local f, fp = (a+pmt), (a+b)*0.5
            i = -f*fp/(fp*fp-f*(n*n-1)/12*a)
            return i, i
        end
        local x = i/expm1(n*log1p(i))
        local k, y = (pv+fv)*x, n*x-1
        local f = k + pv*i + pmt
        local num, den = y+(n-1)*i, i+i*i
        x = k*num - pv*den      -- f' * -den
        y = k*(num+i+1)*(num+y) -- f'' * den^2
        x = f*x*den / (x*x-.5*f*y)
        i = i + x   -- Halley's method
        return i, x
    end
end

Previous post example:

lua> n, pv, pmt, fv = 10, 50, -30, 100
lua> g = loan_rate2(n,pv,pmt,fv)
lua> g()
-0.2844446981069045        0.015555301893095532
lua> g()
-0.28443599888025756      8.699226646944535e-006
lua> g()
-0.28443599888025595      1.6279605736411871e-015

Fun Math Algorithms, car lease APR example.

lua> n, pv, pmt, fv = 36, 30000, -550, -15000
lua> edge_i(n,pv,pmt,fv) -- Note: 2nd edge can be removed
0.018333333333333333       0.03666666666666667
lua>
lua> g = loan_rate2(n,pv,pmt,fv)
lua> g()
0.005817405851432355      -0.012515927481900979
lua> g()
0.005805072819430545      -1.2333032001809594e-005
lua> g()
0.0058050728194201295    -1.0415537616188019e-014

Note sign of errors.
Iteration does not over-shoot to the other side.
Next iteration guaranteed better estimate (cubic convergence)

Bonus: with one-sided convergence pattern, we can detect for no solution.
If error sign changes is not due to simple rounding errors, no solution.

lua> g = loan_rate2(10,50,-10,100)
lua> g()
0.3175682256333693      0.4175682256333693
lua> g()
-0.0740994732881643     -0.3916676989215336 --> no solution

[Albert Chan]

(05-17-2022 01:31 PM)Albert Chan Wrote:  f = (pv+fv)/n * g + pv * x + pmt = 0      // same formula used for Plus42

g = n*x / ((1+x)^n-1)
g'/g = -(g-1+(n-1)*x) / (x+x^2) = -num / den
g''/g = (g+n*x)*(2*(g-1)+(n-1)*x) / (x+x^2)^2 = (num+x+1)*(num+(g-1)) / den^2

Let k = (pv+fv)/n * g

f = k + pv*x + pmt
f' = k*(g'/g) + pv
f'' = k*(g''/g)

Using first derivative to get to second, Halley's correction is cheap.
To make code more robust, we special cased with limits when rate goes 0

If solution exist, Halley's method seems not to over-shoot, but not yet proved.
However, starting from edge, we can proof no over-shoot for Newton's method

Proof: f is either concave-up, or concave-down (f'' have same sign)

g = n*x / ((1+x)^n-1) = 1 - (n-1)/2*x + (n²-1)/12*x² - (n²-1)/24*x³ + ...

For integer n, g is decaying if n>0, growing if n<0 (= compounding factor C)
Note: if n=1, g=1; if n=-1, g=1+x ⇒ For n=±1, f is a straight line.

g ≥ 1 - (n-1)/2*x            // x=0 tangent line, but inequality holds for other x.
g + (n-1)/2*x - 1 ≥ 0      // if x=0, g=1 ⇒ LHS=0

g''/g = (g+n*x)*(2*(g-1)+(n-1)*x) / (x+x^2)^2
g'' = g^2 * (1+x)^n * 2 * (g+(n-1)/2*x-1) / (x+x^2)^2

Assume n>1:

For -1 < x < ∞, x≠0, RHS all factors are postive ⇒ g'' > 0
Take the limit, g''(x=0) = (n²-1)/6 > 0

g'' is positive thoughtout, except on the edge of x domain:
Take the limit, g''(x=-1) = g''(x=∞) = 0

f'' = (pv+fv)/n * g''      → sign(f'') = sign(pv+fv)

Update:
lua code to automatic rate search, using Newton's 1-sided convergence property.

[Albert Chan]

(06-10-2022 08:25 PM)Albert Chan Wrote:  If solution exist, Halley's method seems not to over-shoot, but not yet proved.

1-sided convergence observation only apply to Halley's rational formula.
Halley's Irrational formula does overshoot, even if guesses not bracketed by 2 roots.

Irrational Halley's modifed slope = f' * (1 + √(1 - 2*f*f''/(f')^2)) / 2

(1 + √(1-2x))/2 = 1 - x/2 - x^2/4 - x^3/4 - ...

Assume correction x is small, cutting off O(x^2) terms made slope correction factor bigger.
Bigger slope made Halley's rational formula correction less aggressive.

Rational Halley's modified slope = f' * (1 - (f*f'')/(f')^2/2) = f' - (f''/2)*(f/f')

In exchange for slightly slower convergence, it is more stable.

(01-07-2020 06:13 PM)Albert Chan Wrote:  Example, for N=30, PV=6500, PMT=-1000, FV=50000, we have I = 8.96% or 11.10%

Does anyone knows what returns should be reported for above investments ?

This is the example I found that Halley's Irrational formula iterations over-shoot.

Code:

Rational    Irrational
0.075       0.075
0.0878885   0.0906481
0.0895961   0.0896056
0.0896059   0.0896059

0.08        0.08 
0.0888735   0.0898506
0.0896050   0.0896059
0.0896059   0.0896059

[Albert Chan]

Examples taken from previous post, where we have 2 rate solutions.

XCAS> tmv := ((pv+fv)/((1+x)^n-1) + pv)*x + pmt
XCAS> f0 := tmv(n=30, pv=6500., pmt=-1000., fv=50000); f1 := f0' :;
XCAS> [a, b] := [limit(x-f0/f1,x=-1), limit(x-f0/f1,x=inf)]

[-1/50, 2/13]            // = 2 edges, [pmt/fv, -pmt/pv]

XCAS> m := fsolve(f1, x=a); f0(x=m)

(0.100019883274, -6.52249630654)              // f extremum

XCAS> [fsolve(f0, x=a), fsolve(f0, x=b)]

[0.0896058562612, 0.111003286405]            // f roots

Extremum(f) = f(x=m) < 0, and f have 2 roots ⇒ f has a concaved-up V shape.

Halley's method, derived from g = f / sqrt(|f'|), Newton correction for g ≡ Halley correction for f

XCAS> g0 := f0 / sqrt(abs(f1)):; g1 := g0':; g2 := g1':;

g have the same roots as f (by its definition)
Interestingly, g'' also have same roots as f

XCAS> fsolve(g2, x=a..b)

[0.0896058562612, 0.111003286405]            // g'' roots = f roots

XCAS> g2(x=a), g2(x=m), g2(x=b)

(796.773155573, -2.55834266547e+38, 371.946352558)

g is concave-down, for intervals between 2 roots, concaved-up otherwise.
⇒ Newton's method on g have 1-sided convergence, for whatever rate guess.
⇒ Halley's method on f have 1-sided convergence, for whatever rate guess.

---

Paradoxically, f and g'' roots matched guaranteed 1-sided convergence is not true.
(at least, not when f has only 1 root)

From the extreme edge (-1,∞), f'' and g'' must have same sign.
We had proved f'' have same sign throughout; g'' signs must be +/-/+, or -/+/-

But, what happens if f has only 1 root ? g'' must still have 2 roots !

Disprove of Halley's 1-sided convergence, by counter-example

We flip above exmple fv sign; f now has only 1 root. (overshoot rate in bold)

lua> iter = loan_rate2(30,6500,-1000,-50000)
lua> for i=1,5 do print(i, iter()) end
1      0.1511929109248802     0.1311929109248802
2      0.16524185831684263    0.014048947391962436
3      0.1652186592665784     -2.319905026425022e-005
4      0.16521865926667817    9.978337255768468e-014
5      0.1652186592666782     2.8509535016486494e-017


This counter-example, g'' roots = [0.0789877153657, 0.165218659267]
g'' second root matched f (only) root.

[Albert Chan]

(06-10-2022 08:25 PM)Albert Chan Wrote:  lua code to automatic rate search, using Newton's 1-sided convergence property.

On 2nd thought, I re-post code here, in case Plus42 repository turned private.

Quote:Newton's method 1-sided convergence is proven. I used it here.

Code:

function edge_i(n,pv,pmt,fv)
    local a, b = pmt/fv, -pmt/pv
    if abs(a)>abs(b) and b>-1 or b~=b then a,b = b,a end
    return a, b  -- a = small edge or NaN
end

function iter_i(n, pv, pmt, fv, i)
    pmt, fv = pmt or -1, fv or 0
    i = i or edge_i(n, pv, pmt, fv)
    return function()
        if i==0 then
            local a, b = (pv+fv)/n, pv-fv
            local f, fp = (a+pmt), (a+b)*0.5
            i = -f/fp
            return i, i, f
        end
        local x = i/expm1(n*log1p(i))
        local k, y = (pv+fv)*x, n*x-1
        local f = k + pv*i + pmt
        local num, den = y+(n-1)*i, i+i*i
        x = f*den / (k*num-pv*den)
        i = i + x   -- Newton's method
        return i, x, f
    end
end

function find_rate(n,pv,pmt,fv,BEGIN,i0)
    if BEGIN then pv, fv = pv+pmt, fv-pmt end
    local g = iter_i(n,pv,pmt,fv,i0)
    if i0 then g() end  -- throw away iteration
    local i,eps,f = g() -- 1-sided convergence
    repeat
        local f0 = f
        i,eps,f = g()
        if f*f0 <= 0 then return i,eps end
    until not (abs(f0) > abs(f))
    return nil
end

Note: NaN edge moved up front, for rate search early termination.
All examples shown in this issue.

Code:

lua> find_rate(10,50,-30,80)
-0.3689336987417774     0
lua> find_rate(10,50,-30,100)
-0.284435998880256      -4.284106772739964e-017
lua> find_rate(10,50,-30,100,true) -- BEGIN mode
-0.20399537076838428    0
lua> find_rate(36,-10000,500,-6000)
0.018066231336911785    1.454230867769741e-017
lua> find_rate(36,-10000,500,-9000)
nil

Update: find_rate() now allow user to enter rate guess.
This is rarely needed, unless we wanted solution (if any) from the other edge.
To maintain 1-sided convergence, first possibly overshooted iteration throw away.

Example below used guesses between 2 roots.
Note: the safe guess is again start from the edge.
Newton's method may fail, if guess get too near where f'=0.

Code:

lua> find_rate(30,6500,-1000,50000,nil,0.09)
0.08960585626123246     0
lua> find_rate(30,6500,-1000,50000,nil,0.11)
0.11100328640549173     -0

[Albert Chan]

Testing for f = 0 (or changed sign), with *guaranteed* root, is good to have.
However, it is possible eps too small to change iterated rate.
This cause f unable to improve (toward 0), giving false negative result.

lua> find_rate(300,-1000,0,10002) -- ???
nil
lua> g = iter_i(300,-1000,0,10002) -- to see what's under the hood
lua> for i=1,9 do print(g()) end

Code:

0.005469683722009796    0.005469683722009796    30.006666666666668
0.007458966325604198    0.001989282603594402    6.432723957167589
0.007702303993544976    0.00024333766794077768  0.636567520936282
0.007705485339429362    3.1813458843856375e-006 0.008111326572127986
0.007705485872366111    5.329367491485611e-010  1.3583484603785223e-006
0.007705485872366127    1.5332655972388135e-017 3.907985046680551e-014
0.007705485872366127    3.484694539179127e-019  8.881784197001252e-016
0.007705485872366127    3.484694539179127e-019  8.881784197001252e-016

(i + eps) rounding errors caused iteration not to advance, f unable to change sign.

For above example, relative to i, 1/2 ULP = 0x1p-61 ≈ 4.337e-19 ≈ 1.25 eps
(that's in a perfect world; i uncertainty is likely > 1/2 ULP)

It is not a guarantee, but if rate converged, we should also consider solution found.

Revised find_rate

Code:

function find_rate(n,pv,pmt,fv,BEGIN,i0)
    if BEGIN then pv, fv = pv+pmt, fv-pmt end
    local g = iter_i(n,pv,pmt,fv,i0)
    if i0 then g() end  -- throw away iteration
    local i,eps,f = g() -- 1-sided convergence
    repeat
        local f0 = f
        i,eps,f = g()
        if f*f0 <= 0 then return i,eps end
    until not (abs(f0) > abs(f))
    return i-eps, (i == i-eps*0.01)
end

Code assumed about 15 digits matched "converged".
Last iteration, with |f| not improving, throw away.

lua> find_rate(300,-1000,0,10002) -- revised version
0.007705485872366127      true

[Albert Chan]

(06-12-2022 05:07 PM)Albert Chan Wrote:  Halley's method, derived from g = f / sqrt(|f'|), Newton correction for g ≡ Halley correction for f

XCAS> g0 := f0 / sqrt(abs(f1)):; g1 := g0':; g2 := g1':;

g have the same roots as f (by its definition)
Interestingly, g'' also have same roots as f

To simplify, assume f1^p really mean |f1|^p
https://www.mathsisfun.com/calculus/product-rule.html

g0 = f0 * f1^(-1/2)

g1 = f0 * (-1/2*f1^(-3/2)*f2) + f1 * f1^(-1/2)
     = (f1^2 - f0*f2/2) * f1^(-3/2)

g2 = (f1^2 - f0*f2/2) * (-3/2*f1^(-5/2)*f2) + (2*f1*f2 - (f1*f2+f0*f3)/2) * f1^(-3/2)
     = f0 * (f2^2 - 2/3*f1*f3) * (3/4*f1^(-5/2))

Note: g'' may have more than f roots (2nd factor may also have root)

BTW, we had also proved that Newton correction for g ≡ Halley's correction for f

-g0/g1 = -(f0 * f1^(-1/2)) * f1^(3/2) / (f1^2 - f0*f2/2) = -f0*f1 / (f1^2 - f0*f2/2)

[jonmoore]

Hi Albert, I've been attempting to follow the thread and am presuming that it serves as a communication channel with Thomas ref your ongoing optimisation of the TVM solver in Plus42? The reason I ask is that I noticed you're prototyping with both XCAS and bespoke Lua code. I know KhiCAS on the Nspire can work with the Nspire's Lua interface, so I wondered if that's what you're doing, or is Thomas using Lua as part of the toolset in Plus42?

No real purpose to my enquiry, I'm just interested in the flow of your problem-solving.

[Albert Chan]

This thread is to answer late Dieter's question (10 years ago), about solving TVM=0, for rate.
[WP34s et al.] Solving the TVM equation for the interest rate, have questions but no answer.

f = NPMT/n = ((pv+fv) / ((1+x)^n-1) + pv)*x + pmt      // solve f=0, for x

Question: What is a robust guess ? (if solution exist, guess will leads to it)

From edges, [pmt/fv, -pmt/pv]
For 1 guess, edge_i(n,pv,pmt,fv) will pick the right one.

Question: What kind of iteration scheme to get true rate ?

Newtons method (guaranteed 1-sided convergence)
In other words, starting from 2 edges, we get 2 roots (if roots exist)

---

Plus42 had arbitrary fixed guess, i = 1e-4
Any guess that is close to f' = 0, Newton's method may fail.
I placed an issue to Plus42, suggesting a more robust guess.

Issue: TVM rate guess, was after this thread were created.
(First update to use edge rate guess on June 6, 2022)

---

Regariding Lua, it is just my preferred language, for posting on-line.
It is not white-space sensitive; I can cut/paste the code, and it will just run.

I like Lua's simplicity, a language I can wrap my head to it.
I had patched the language (LuaJIT 1.1.8) to my liking.

I also like Python, but it use indentation to indicate code block.
Some website comment system removes all indentations.
Many times it require editing to make cut/paste code work.

[jonmoore]

Thanks for the info, Albert.

I've dabbled with Lua as it was once the scripting language of choice in Foundry's visual compositor Nuke. These days the Python API is far better than Lua as the Lua API is no longer maintained, so I don't come into contact with it so often.

A lot of content creation packages have moved away from Lua, which is a pity because Lua is nice and fast, which is super important when you've got 30,000 frames to render in a very short timeframe!

[Albert Chan]

Hi, jonmoore

I've translated Lua code (Newton's method version) to Python.

Code:

from math import *

def expm1(x):           # delete if math library have expm1
    y = exp(x)
    return y-1 - y*(log(y)-x)

def edge_i(n, pv, pmt, fv):
    a = 1/pv if pv else float('inf')
    b = 1/fv if fv else float('inf')
    a, b = pmt*-a, pmt*b
    if abs(a)>abs(b) and b>-1 or b!=b: a,b = b,a
    return a, b         # small edge or NaN

def iter_i(n, pv, pmt=-1, fv=0, i=None, verbose=False):
    if i == None: i = edge_i(n, pv, pmt, fv)[0]
    while 1:
        try:
            if i==0:
                a, b = (pv+fv)/n, pv-fv
                f, fp = (a+pmt), (a+b)*0.5
                i = x = -f/fp
            else:
                x = i/expm1(n*log1p(i))
                k, y = (pv+fv)*x, n*x-1
                f = k + pv*i + pmt
                num, den = y+(n-1)*i, i+i*i
                x = f*den / (k*num-pv*den)
                i = i + x   # Newton's method
        except (ArithmeticError, ValueError):
            i = x = f = float('nan')
        if verbose: print repr(i), x, f
        yield i, x, f

def find_rate(n, pv, pmt, fv, i=None, verbose=False):
    g = iter_i(n, pv, pmt, fv, i, verbose).next
    if i != None: g()   # throw away iteration
    i,eps,f = g()       # 1-sided convergence
    while 1:
        f0 = f
        i,eps,f = g()
        if f*f0 <= 0: return i,eps
        if not(abs(f0) > abs(f)): return i-eps, (i == i-eps*0.01)

>>> find_rate(36,-10000,500,-6000)
(0.018066231336911796, 0.0)
>>> find_rate(36,-10000,500,-9000) # no solution
(-0.25728394143342476, False)

Update:

Python, unlike Lua, may not return NaN, for "bad" calculations.
I've wrapped iter_i() with a try ... except clause, to catch them all (I hope)

>>> n,pv,pmt,fv = 10,50,-30,80
>>> find_rate(n,pv,pmt,fv, i=1e-4) # Plus42 original fixed guess, failed search
(nan, False)

Udpate2: added verbose option

>>> edge_i(n,pv,pmt,fv)
(-0.375, 0.59999999999999998)

>>> find_rate(n,pv,pmt,fv, verbose=1) # default i = small edge
-0.36895081273778957 0.00604918726221 0.447448184647
-0.368933698883462 1.71138543276e-05 0.00125870008025
-0.36893369874177745 1.41684553346e-10 1.04205284401e-08
-0.36893369874177739 4.8305098317e-17 3.5527136788e-15
-0.36893369874177739 0.0 0.0
(-0.36893369874177739, 0.0)

>>> find_rate(n,pv,pmt,fv, verbose=1, i=0.6) # the other edge
0.58466252739896962 -0.015337472601 0.715917095435
0.58461955301798252 -4.29743809871e-05 0.00199461485041
0.58461955266220977 -3.55772689766e-10 1.65125761953e-08
0.58461955266220988 7.65452032864e-17 -3.5527136788e-15
(0.58461955266220988, 7.6545203286364044e-17)

Udpate3: changed i default to None, instead of False
Python False == 0 == 0.0, but only None equal None.

[Thomas Okken]

Albert, thank you very much for these algorithms!

I just released Plus42 1.0.8 which incorporates these latest improvements.

[Albert Chan]

(06-10-2022 08:25 PM)Albert Chan Wrote:  g = n*x / ((1+x)^n-1) = 1 - (n-1)/2*x + (n²-1)/12*x² - (n²-1)/24*x³ + ...

For integer n, g is decaying if n>0, growing if n<0 (= compounding factor C)
Note: if n=1, g=1; if n=-1, g=1+x ⇒ For n=±1, f is a straight line.

g ≥ 1 - (n-1)/2*x            // x=0 tangent line, but inequality holds for other x.

We relied on above identity to prove g'' > 0 ⇒ f'' has same sign throughout.

It is easy to show with plots (pick any n>1), but harder to proof.
(Perhaps I missed something simple ...)

I had placed this question of proving it here (response welcome)

I managed to proof it, using XCAS.

Let R = 1+x > 0, to prove:

g = n*(R-1)/(R^n-1) ≥ 1 - (n-1)/2*(R-1)

→ Or, n ≥ (1 - (n-1)/2*(R-1)) * (R^n-1)/(R-1)      // note: multiplied factor > 0

(RHS, R)' = (1 - n*(n+1)/2*R^(n-1) + (n²-1)*R^n - n*(n-1)/2*R^(n+1)) / (R-1)^2

If n > 1, numerator have 3 sign changes, 1 or 3 positive roots (Descartes rules of signs)
Taking taylor series, we showed that it had triple equal roots, at x=0

XCAS> M := (1 - n*(n+1)/2*R^(n-1) + (n²-1)*R^n - n*(n-1)/2*R^(n+1))
XCAS> series(M(R=1+x), x, 0, 4, polynom)

(n/6-n^3/6)*x^3 + (-n/4+n^2/8+n^3/4-n^4/8)*x^4

2 roots get cancelled by denominator, x^2, we have only 1 extremum, at x=0
Taking limit at x=0, we have

LHS = g(0) = 1
RHS = 1 - (n-1)/2*0 = 1         → LHS = RHS

For x≠0, two sides don't touch. Example, pick x=∞

LHS = g(∞) = 0 (decay to nothing)
RHS = 1 - (n-1)/2*∞ = -∞      → LHS > RHS

For x > -1, LHS ≥ RHS      QED

[Albert Chan]

(06-19-2022 06:57 PM)Albert Chan Wrote:  Let R = 1+x > 0, to prove:

g = n*(R-1)/(R^n-1) ≥ 1 - (n-1)/2*(R-1)

→ Or, n ≥ (1 - (n-1)/2*(R-1)) * (R^n-1)/(R-1)      // note: multiplied factor > 0

(RHS, R)' = (1 - n*(n+1)/2*R^(n-1) + (n²-1)*R^n - n*(n-1)/2*R^(n+1)) / (R-1)^2

If n > 1, numerator have 3 sign changes, 1 or 3 positive roots (Descartes rules of signs)
Taking taylor series, we showed that it had triple equal roots, at x=0 ...

I was being stupid.

Multiplied factor is a polynomial = (R^(n-1) + R^(n-2) + ... + R + 1)
In other words, RHS is just a polynomial. All its derivatives also polynomial.

(RHS,R)' roots = (1 or 3) - 2 = (-1 or 1) = 1 positive root --> 1 extremum (*)

No CAS needed for the proof

(*) we can say extremum must be where curve and tangent line touched.
In other words, numerator M have triple equal roots, at R = 1

Or, we do the math, and confirm M's 3 roots:

M = 1 - n*(n+1)/2*R^(n-1) + (n^2-1)*R^n - n*(n-1)/2*R^(n+1)

M(R=1) = 1 - n*(n+1)/2 + (n^2-1) - (n^2-n)/2 = 0      // M first root at R=1
M' = -(n+1)*n*(n-1)/2 * R^(n-2) * (R-1)^2                  // other 2 roots, also R=1

[Albert Chan]

(06-18-2022 11:03 PM)Thomas Okken Wrote:  Albert, thank you very much for these algorithms!

I just released Plus42 1.0.8 which incorporates these latest improvements.

TVM improvement nicely remove the need for iterator, using counter to ensure 2 f's before rate search.

Here is an updated Python find_rate(), using sentinel values, instead of counter.

Code:

def find_rate(n, pv, pmt, fv, i=None, verbose=False):
    i0, f0 = i, float('nan')
    for (i,eps,f) in iter_i(n, pv, pmt, fv, i, verbose):
        if i0!=None: i0=None; continue  # skip 1 iteration
        if f != f:    return i, False   # bad f -> bad i
        if f0*f <= 0: return i, True    # locked-in rate
        if abs(f) < abs(f0)): f0=f; continue
        return i-eps, (i == i-eps*0.01)

Above relied on IEEE-754 NaN special propety:
(NaN != x) returns True; all other kinds of NaN comparisons returns False.

We could also remove NaN's ahead of time, using function math.isnan()

Code:

from math import isnan
        
def find_rate(n, pv, pmt, fv, i=None, verbose=False):
    i0, f0 = i, float('nan')
    for (i,eps,f) in iter_i(n, pv, pmt, fv, i, verbose):
        if i0!=None: i0=None; continue  # skip 1 iteration
        if isnan(f):  return i, False   # bad f -> bad i
        if isnan(f0): f0=f; continue    # ensure 2 f's
        if f0*f <= 0: return i, True    # locked-in rate
        if abs(f) < abs(f0): f0=f; continue
        return i-eps, (i == i-eps*0.01)

Udpate: using counter is shorter and simpler

Code:

def find_rate(n, pv, pmt, fv, i=None, verbose=False, c=2):
    for (i,eps,f) in iter_i(n, pv, pmt, fv, i, verbose):
        if c: c-=1; f0=f; continue      # skip iterations
        if f0*f <= 0: return i, True    # locked-in rate
        if abs(f) < abs(f0): f0=f; continue
        return i-eps, (i == i-eps*0.01)