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:

NPMT = 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

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

Ce is also 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:

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

NPMT = 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:

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

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

We can also solve NPMT=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:  NPMT = A*x/t + B*x + pmt

NPMT = 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 x ranges: (x > pmt/fv) || (x < -pmt/pv)

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

Update: this example, we can remove (x > pmt/fv)

t = tanh(n/2*log1p(x))

Let rate edge X = x(t=±1) = -pmt / (A/±1 + B)

Any finite x, same sign as X, |x| ≥ |X|, will end up with |t| < 1
Next iteration, x will always get "smaller" (x will never converge to itself)

Rule of thumb, if both edges have same sign, pick the smaller sized one.
Above example, based on tanh info only, x = (-1, pmt/pv) ≈ (-1, 0.01833)

[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]