Repository: numpy/numpy-financial

Issue:
Slow NPV calculations for monte carlo simulation purposes
When using a large number of monte carlo simulations (~1e6 or more) for uncertainty calculations on NPV, the present implementation of Numpy-financial NPV is very slow. We suggest a new implementation of the NPV function , which allows calculation of serveral projects simultaneously, and is approximately 200 times faster than the current version for 1e6 simulations for a cash-flow of length 10.

In cases where the number of entries in the cash-flow is significantly higher than the number of projects to be calculated, the old implementation will be faster. To solve this, we suggest an adaptive approach which chooses method based on the input dimensions of the 'values' array.

Example code:
```
import numpy_financial as npf
import numpy as np
import time

def faster_npv(rate: float, values: list):
    """ A faster way to calculate NPV for several projects based 
    on numpy arrays. """
    discounted_values = []
    for i in range(np.shape(values)[0]):
        discounted_values.append(values[i] / (1 + rate) ** i)
    return np.sum(discounted_values, axis=0)

rate=0.05
no_simulations = int(1e6)

capex = np.random.normal(loc=10e5, scale = 10e3, size=no_simulations)
lifetime = 10 #years

opex = []
for yr in range(lifetime):
    opex.append(np.random.normal(loc=10e4, scale= 10e2,size=no_simulations))

start_time=time.time()
all_npv=np.nan*np.ones_like(capex)
for i in range(no_simulations):
    values = [capex[i]]+[op[i] for op in opex]
    all_npv[i]=npf.npv(rate, values)
print(f"Standard NPF calculations took {time.time()-start_time} s.")

start_time=time.time()
values=[capex]+opex
faster_npv(rate,values )
print(f"Faster NPV calculations took {time.time()-start_time} s.")
```


Discussion:

Patch:
diff --git a/numpy_financial/_financial.py b/numpy_financial/_financial.py
index e268424..070599d 100644
--- a/numpy_financial/_financial.py
+++ b/numpy_financial/_financial.py
@@ -675,7 +675,36 @@ def rate(nper, pmt, pv, fv, when='end', guess=None, tol=None, maxiter=100):
     return rn
 
 
-def irr(values, guess=0.1, tol=1e-12, maxiter=100):
+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):
     """
     Return the Internal Rate of Return (IRR).
 
@@ -692,13 +721,6 @@ def irr(values, guess=0.1, tol=1e-12, maxiter=100):
         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 an initial guess of 0.1 (i.e. 10%) is assumed instead.
-    tol : float, optional
-        Required tolerance to accept solution. Default is 1e-12.
-    maxiter : int, optional
-        Maximum iterations to perform in finding a solution. Default is 100.
 
     Returns
     -------
@@ -708,7 +730,7 @@ def irr(values, guess=0.1, tol=1e-12, maxiter=100):
     Notes
     -----
     The IRR is perhaps best understood through an example (illustrated
-    using np.irr in the Examples section below). Suppose one invests 100
+    using np.irr in the Examples section below).  Suppose one invests 100
     units and then makes the following withdrawals at regular (fixed)
     intervals: 39, 59, 55, 20.  Assuming the ending value is 0, one's 100
     unit investment yields 173 units; however, due to the combination of
@@ -749,36 +771,28 @@ def irr(values, guess=0.1, tol=1e-12, maxiter=100):
     if values.ndim != 1:
         raise ValueError("Cashflows must be a rank-1 array")
 
-    # If all values are of the same sign no solution exists
-    # we don't perform any further calculations and exit early
-    same_sign = np.all(values > 0) if values[0] > 0 else np.all(values < 0)
-    if same_sign:
-        return np.nan
+    # 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]
 
-    # We aim to solve eirr such that NPV is exactly zero. This can be framed as
-    # simply finding the closest root of a polynomial to a given initial guess
-    # as follows:
-    #           V0           V1           V2           V3
-    # NPV = ---------- + ---------- + ---------- + ---------- + ...
-    #       (1+eirr)^0   (1+eirr)^1   (1+eirr)^2   (1+eirr)^3
-    #
-    # by letting x = 1 / (1+eirr), we substitute to get
-    #
-    # NPV = V0 * x^0   + V1 * x^1   +  V2 * x^2  +  V3 * x^3  + ...
-    # 
-    # which we solve using Newton-Raphson and then reverse out the solution 
-    # as eirr = 1/x - 1 (if we are close enough to a solution)
-    npv_ = np.polynomial.Polynomial(values)
-    d_npv = npv_.deriv()
-    x = 1 / (1 + guess)
+    res = _roots(values[::-1])
 
-    for _ in range(maxiter):
-        x_new = x - (npv_(x) / d_npv(x))
-        if abs(x_new - x) < tol:
-            return 1 / x_new - 1
-        x = x_new
+    mask = (res.imag == 0) & (res.real > 0)
+    if not mask.any():
+        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)
 
-    return np.nan
+    rate = rate.item(np.argmin(np.abs(rate)))
+    return rate
 
 
 def npv(rate, values):
@@ -849,8 +863,15 @@ def npv(rate, values):
     3065.22267
 
     """
-    values = np.asarray(values)
-    return (values / (1+rate)**np.arange(0, len(values))).sum(axis=0)
+    values = np.atleast_2d(values)
+    timestep_array = np.arange(0, values.shape[1])
+    npv = (values / (1 + rate) ** timestep_array).sum(axis=1)
+    try:
+        # If size of array is one, return scalar
+        return npv.item()
+    except ValueError:
+        # Otherwise, return entire array
+        return npv
 
 
 def mirr(values, finance_rate, reinvest_rate):


Test patch:
diff --git a/numpy_financial/tests/test_financial.py b/numpy_financial/tests/test_financial.py
index 791f1a7..35ac855 100644
--- a/numpy_financial/tests/test_financial.py
+++ b/numpy_financial/tests/test_financial.py
@@ -148,6 +148,19 @@ class TestNpv:
             npf.npv(Decimal('0.05'), [-15000, 1500, 2500, 3500, 4500, 6000]),
             Decimal('122.894854950942692161628715'))
 
+    def test_npv_broadcast(self):
+        cashflows = [
+            [-15000, 1500, 2500, 3500, 4500, 6000],
+            [-15000, 1500, 2500, 3500, 4500, 6000],
+            [-15000, 1500, 2500, 3500, 4500, 6000],
+            [-15000, 1500, 2500, 3500, 4500, 6000],
+        ]
+        expected_npvs = [
+            122.8948549, 122.8948549, 122.8948549, 122.8948549
+        ]
+        actual_npvs = npf.npv(0.05, cashflows)
+        assert_allclose(actual_npvs, expected_npvs)
+
 
 class TestPmt:
     def test_pmt_simple(self):
@@ -564,7 +577,7 @@ class TestFv:
 class TestIrr:
     def test_npv_irr_congruence(self):
         # IRR is defined as the rate required for the present value of
-        # a series of cashflows to be zero, so we should have
+        # a a series of cashflows to be zero, so we should have
         #
         # NPV(IRR(x), x) = 0.
         cashflows = numpy.array([-40000, 5000, 8000, 12000, 30000])
@@ -593,12 +606,9 @@ class TestIrr:
             decimal=2,
         )
 
-    @pytest.mark.parametrize('v', [
-        (1, 2, 3),
-        (-1, -2, -3),
-    ])
-    def test_numpy_gh_6744(self, v):
+    def test_numpy_gh_6744(self):
         # Test that if there is no solution then npf.irr returns nan.
+        v = [-1, -2, -3]
         assert numpy.isnan(npf.irr(v))
 
     def test_gh_15(self):