ENH: Use Newton's method to calculate IRR

Author: Kai-Striega

numpy_financial/_financial.py

@@ -671,36 +671,7 @@ def rate(nper, pmt, pv, fv, when='end', guess=None, tol=None, maxiter=100):
         return rn
 
 
-def _roots(p):
-    """Modified version of NumPy's roots function.
-
-    NumPy's roots uses the companion matrix method, which divides by
-    p[0]. This can causes overflows/underflows. Instead form a
-    modified companion matrix that is scaled by 2^c * p[0], where the
-    exponent c is chosen to balance the magnitudes of the
-    coefficients. Since scaling the matrix just scales the
-    eigenvalues, we can remove the scaling at the end.
-
-    Scaling by a power of 2 is chosen to avoid rounding errors.
-
-    """
-    _, e = np.frexp(p)
-    # Balance the most extreme exponents e_max and e_min by solving
-    # the equation
-    #
-    # |c + e_max| = |c + e_min|.
-    #
-    # Round the exponent to an integer to avoid rounding errors.
-    c = int(-0.5 * (np.max(e) + np.min(e)))
-    p = np.ldexp(p, c)
-
-    A = np.diag(np.full(p.size - 2, p[0]), k=-1)
-    A[0,:] = -p[1:]
-    eigenvalues = np.linalg.eigvals(A)
-    return eigenvalues / p[0]
-
-
-def irr(values):
+def irr(values, guess=0.1):
     """
     Return the Internal Rate of Return (IRR).
 
@@ -717,6 +688,9 @@ def irr(values):
         are negative and net "withdrawals" are positive.  Thus, for
         example, at least the first element of `values`, which represents
         the initial investment, will typically be negative.
+    guess : float, optional
+        Initial guess of the IRR for the iterative solver. If no guess is
+        given 0.1 is used instead.
 
     Returns
     -------
@@ -767,28 +741,23 @@ def irr(values):
     if values.ndim != 1:
         raise ValueError("Cashflows must be a rank-1 array")
 
-    # Strip leading and trailing zeros. Since we only care about
-    # positive roots we can neglect roots at zero.
-    non_zero = np.nonzero(np.ravel(values))[0]
-    values = values[int(non_zero[0]):int(non_zero[-1])+1]
+    solution_found = False
+    p = np.polynomial.Polynomial(values)
+    pp = p.deriv()
+    x = 1 / (1 + guess)
 
-    res = _roots(values[::-1])
+    for i in range(100):
+        x_new = x - (p(x) / pp(x))
+        if abs(x_new - x) < 1e-12:
+            solution_found = True
+            break
+        x = x_new
 
-    mask = (res.imag == 0) & (res.real > 0)
-    if not mask.any():
+    if solution_found:
+        return 1 / x - 1
+    else:
         return np.nan
-    res = res[mask].real
-    # NPV(rate) = 0 can have more than one solution so we return
-    # only the solution closest to zero.
-    rate = 1/res - 1
-
-    # If there are any positive solutions prefer those over negative
-    # rates.
-    if (rate > 0).any():
-        rate = np.where(rate > 0, rate, np.inf)
-
-    rate = rate.item(np.argmin(np.abs(rate)))
-    return rate
+
 
 
 def npv(rate, values):