Bossy Lobster

A blog by Danny Hermes; musing on tech, mathematics, etc.

Edit on GitHub

Reverse Calculating An Interest Rate

I was recently playing around with some loan data and only happened to have the term (or length, or duration) of the loan, the amount of the recurring payment (in this case monthly) and the remaining principal owed on the loan. I figured there was an easy way to get at the interest rate, but wasn't sure how. After some badgering from my coworker +Paul, I searched the web and found a tool from CALCAmo (a site just for calculating amortizations).

Problem solved, right? Wrong. I wanted to know why; I had to go deeper. So I did a bit of math and a bit of programming and I was where I needed to be. I'll break the following down into parts before going on full steam.

  • Break down the amortization schedule in terms of the variables we have and the one we want
  • Determine a function we want to find zeros of
  • Write some code to implement the Newton-Raphson method
  • Utilize the Newton-Raphson code to find an interest rate
  • Bonus: Analyze the function to make sure we are right

Step I: Break Down the Amortization Schedule

We can do this using the series {Pi}i\left\{P_i\right\}_i of principal owed, which varies over time and will go to zero once paid off. In this series, P0P_0 is the principal owed currently and PiP_i is the principal owed after ii payments have been made. (Assuming monthly payments, this will be after ii months.) If the term is TT periods, then we have PT=0P_T = 0.

We have already introduced the term TT; we also need the value of the recurring (again, usually monthly) payment R,R, the interest rate rr and the initial principal owed P0=PP_0 = P.

Time-Relationship between Principal Values

If after ii periods, PiP_i is owed, then after one period has elapsed, we will owe PimP_i \cdot m where m=m(r)m = m(r) is some multiplier based on the length of the term. For example if each period is one month, then we divide our rate by 1212 for the interest and add 11 to note that we are adding to existing principal:

m(r)=1+r12.m(r) = 1 + \frac{r}{12}.

In addition to the interest, we will have paid off RR hence

Pi+1=PimR.P_{i + 1} = P_i \cdot m - R.

Formula for PiP_i

Using this, we can actually determine PiP_i strictly in terms of m,Rm, R and PP. First, note that

P2=P1mR=(P0mR)mR=Pm2R(m+1)\begin{aligned} P_2 = P_1 \cdot m - R &= (P_0 \cdot m - R) \cdot m - R \\ &= P \cdot m^2 - R(m + 1) \end{aligned}

since P0=PP_0 = P. We can show inductively that

Pi=PmiRj=0i1mj.\displaystyle P_i = P \cdot m^i - R \cdot \sum_{j = 0}^{i - 1} m^j.

We already have the base case i=1,i = 1, by definition. Assuming it holds for i,i, we see that

Pi+1=PimR=m(PmiRj=0i1mj)R=Pmi+1Rj=1imjR,\begin{aligned} P_{i + 1} = P_i \cdot m - R &= m \cdot \left(P \cdot m^i - R \cdot \sum_{j = 0}^{i - 1} m^j\right) - R \\ &= P \cdot m^{i + 1} - R \cdot \sum_{j = 1}^{i} m^j - R, \end{aligned}

and our induction is complete. (We bump the index jj since we are multiplying each mjm^j by mm.) Each term in the series is related to the previous one (except P0,P_0, since time can't be negative in this case).

Step II: Determine a Function we want to find Zeros of

Since we know PT=0P_T = 0 and PT=PmTRj=0T1mj,\displaystyle P_T = P \cdot m^T - R \cdot \sum_{j = 0}^{T - 1} m^j, we actually have a polynomial in place that will let us solve for mm and in so doing, solve for rr.

To make our lives a tad easier, we'll do some rearranging. First, note that

j=0T1mj=mT1++m+1=mT1m1.\displaystyle \sum_{j = 0}^{T - 1} m^j =m^{T - 1} + \cdots + m + 1 = \frac{m^T - 1}{m - 1}.

We calculate this sum of a geometric series here, but I'll just refer you to the Wikipedia page instead. With this reduction we want to solve

0=PmTRmT1m1PmT(m1)=R(mT1).\begin{aligned} & 0 = P \cdot m^T - R \cdot \frac{m^T - 1}{m - 1} \\ & \Longleftrightarrow P \cdot m^T \cdot (m - 1) =R \cdot(m^T - 1). \end{aligned}

With that, we have accomplished Step II, we have found a function (parameterized by P,TP, T and RR which we can use zeros from to find our interest rate:

fP,T,R(m)=PmT(m1)R(mT1)=PmT+1(P+R)mT+R.\begin{aligned} f_{P, T, R}(m) &= P \cdot m^T \cdot (m - 1) -R \cdot(m^T - 1) \\ &= P \cdot m^{T + 1} - (P + R) \cdot m^T + R. \end{aligned}

Step III: Write some code to implement the Newton-Raphson method

We use the Newton-Raphson method to get super-duper-close to a zero of the function.For in-depth coverage, see the Wikipedia page on the Newton-Raphson method, but I'll give some cursory coverage below. The methods used to show that a fixed point is found are not necessary for the intuition behind the method.

Intuition behind the method

For the intuition, assume we know (and can compute) a function f,f, its derivative ff' at a value xx. Assume there is some zero yy nearby xx. Since they are close, we can approximate the slope of the line between the points (x,f(x))(x, f(x)) and (y,f(y))(y, f(y)) with the derivative nearby. Since we know x,x, we use f(x)f'(x) and intuit that

f(x)=slope=f(y)f(x)yxyx=f(y)f(x)f(x).f'(x) = \text{slope} = \frac{f(y) - f(x)}{y - x} \Rightarrow y - x = \frac{f(y) - f(x)}{f'(x)}.

But, since we know that yy is a zero, f(y)f(x)=f(x)f(y) - f(x) = -f(x) hence

yx=f(x)f(x)y=xf(x)f(x).y - x = \frac{-f(x)}{f'(x)} \Rightarrow y = x - \frac{f(x)}{f'(x)}.

Using this method, one can start with a given value x0=xx_0 = x and compute better and better approximations of a zero via the iteration above that determines yy. We use a sequence to do so:

xi+1=xif(xi)f(xi)x_{i + 1} = x_i - \frac{f(x_i)}{f'(x_i)}

and stop calculating the xix_i either after f(xi)f(x_i) is below a preset threshold or after the fineness of the approximation xixi+1\left|x_i - x_{i + 1}\right| goes below a (likely different) preset threshold. Again, there is much that can be said about these approximations, but we are trying to accomplish things today, not theorize.

Programming Newton-Raphson

To perform Newton-Raphson, we'll implement a Python function that takes the initial guess x0x_0 and the functions ff and ff'. We'll also (arbitrarily) stop after the value f(xi)f(x_i) drops below 10810^{-8} in absolute value.

def newton_raphson_method(guess, f, f_prime):
    def next_value(value):
        return value - f(value)*1.0/f_prime(value)

    current = guess
    while abs(f(current)) > 10**(-8):
        current = next_value(current)

    return current

As you can see, once we have f and f_prime, everything else is easy because all the work in calculating the next value (via next_value) is done by the functions.

Step IV: Utilize the Newton-Raphson code to find an Interest Rate

We first need to implement fP,T,R(m)=PmT+1(P+R)mT+Rf_{P, T, R}(m) = P \cdot m^{T + 1} - (P + R) \cdot m^T + R and fP,T,Rf'_{P, T, R} in Python. Before doing so, we do a simple derivative calculation:

fP,T,R(m)=P(T+1)mT(P+R)TmT1.f_{P, T, R}'(m) = P \cdot (T + 1) \cdot m^T - (P + R) \cdot T \cdot m^{T - 1}.

With these formulae in hand, we write a function which will spit out the corresponding f and f_prime given the parameters PP (principal), TT (term) and RR (payment):

def generate_polynomials(principal, term, payment):
    def f(m):
        return (principal*(m**(term + 1)) - (principal + payment)*(m**term) +

    def f_prime(m):
        return (principal*(term + 1)*(m**term) -
                (principal + payment)*term*(m**(term - 1)))

    return (f, f_prime)

Note that these functions only take a single argument (m), but we are able to use the other parameters from the parent scope beyond the life of the call to generate_polynomials due to closure in Python.

In order to solve, we need an initial guess, but we need to know the relationship between mm and rr before we can determine what sort ofguessmakes sense. In addition, once a value for mm is returned from Newton-Raphson, we need to be able to turn it into an rr value so functions m and m_inverse should be implemented. For our dummy case here, we'll assume monthly payments (and compounding):

def m(r):
    return 1 + r/12.0

def m_inverse(m_value):
    return 12.0*(m_value - 1)

Using these, and assuming that an interest rate of 10% is a good guess, we can put all the pieces together:

def solve_for_interest_rate(principal, term, payment, m, m_inverse):
    f, f_prime = generate_polynomials(principal, term, payment)

    guess_m = m(0.10)  # ten percent as a decimal
    m_value = newton_raphson_method(guess_m, f, f_prime)
    return m_inverse(m_value)

To check that this makes sense, let's plug in some values. Using the loan calculator, if we have a 30-year loan (with 1230=36012 \cdot 30 = 360 months of payments) of $100,000 with an interest rate of 7%, the monthly payment would be $665.30. Plugging this into our pipeline:

>>> principal = 100000
>>> term = 360
>>> payment = 665.30
>>> solve_for_interest_rate(principal, term, payment, m, m_inverse)

And we see the rate of 7% is approximated quite well!

Bonus: Analyze the function to make sure we are right

Coming soon. We will analyze the derivative and concavity to make sure that our guess yield the correct (and unique) zero.