User defined XIRR function: exception handling

Question

This is my first attempt at Python. I started with some code from https://github.com/peliot/XIRR-and-XNPV (which aims to mimic the Excel function XIRR) and tried to 'improve' on it (and the Excel version). But I would like to know if in these few lines I already made beginner's mistakes. I'm particularly interested to learn whether there is a more elegant way to handle the exceptions (an 'improvement' of mine), because this looks a bit clunky. (This is not the only change compare to the original, so the other bits may be not completely correct either.)

Note that the function newton itself may also give rise to an immediate exception.

import datetime
from scipy import optimize
def xnpv(rate, cashflows):
 t0 = min(cashflows, key = lambda t: t[0])[0]
 return sum([cf/(1+rate)**((t-t0).days/365.0) for (t,cf) in cashflows])
def xirr(cashflows,guess=0.1):
 try:
 outc = optimize.newton(lambda r: xnpv(r, cashflows), guess, 
 maxiter=100)
 if outc.imag == 0:
 return outc
 else:
 raise
 except (RuntimeError, OverflowError):
 try:
 outc = optimize.newton(lambda r: xnpv(r, cashflows), -guess, 
 maxiter=100)
 if outc.imag == 0:
 return outc
 else:
 raise
 except (RuntimeError, OverflowError):
 return float("NaN")
#cftest = [(datetime.date(2001, 12, 5), -2000), (datetime.date(2007, 12, 5), -10), (datetime.date(2017, 12, 5), 20)]
cftest = [(datetime.date(2001, 12, 5), -2000), (datetime.date(2000, 12, 5), -1000), (datetime.date(2017, 12, 5), 2000)]
#cftest = [(datetime.date(2001, 12, 5), -2000), (datetime.date(2007, 12, 5), -1000), (datetime.date(2017, 12, 5), 20)]
print(cftest)
print(xirr(cftest))

Answer 1

Error handling

Let's start with this, because it was the focus of your question.

Exceptions are your friend. They're Python-idiomatic, and have lots of wonderful features; they -

• tell you where the inner problem happened (via the backtrace)
• tell you what kind of problem happened (via the exception type)
• tell you a description of the problem (via the properties of the exception object)
• allow you to unwind the stack in a predictable manner up to a point where they're deliberately caught or allowed to terminate the program
• generally don't pass unnoticed, nor poison later stages of a numerical process

Your except: return NaN undoes all of those advantages. In an exception-friendly language like Python, we should usually have the freedom to write naive code - without checking value validity at the boundary of an API interface. In other words, the consumer of the x*** functions should not be burdened with the need to check for NaN; it should merrily continue on with the values it was given, or allow its execution context to be unwound if an exception has occurred. With a NaN-signalling mechanism, naive code will continue churning through calculations that may accidentally throw in surprising places or may just produce nonsense output, and debugging will be made more difficult. We should strive to develop systems designed to fail early, instead of potentially plowing past the first point of failure. Python is generally a weak language when it comes to fail-fast properties, but for the purposes of this question we'll stick with it.

For the rest of this answer I will assume the use of idiomatic Python exception handling. A reinterpretation of xirr() that uses simple exceptions but preserves the logic may look like

def xirr(cashflows: typing.Iterable[tuple[datetime.date, float]], guess: float = 0.1) -> float:
 def attempt(sign: int) -> float:
 # May raise RuntimeError since disp=True
 return optimize.newton(
 func=xnpv, args=(cashflows,), x0=sign*guess, maxiter=100,
 )
 try:
 return attempt(1)
 except RuntimeError:
 pass # First attempt is non-fatal
 return attempt(-1) # Second attempt is fatal

Science

Notice that newton() varies its algorithm based on the presence or absence of the first- and second-order derivatives:

The Newton-Raphson method is used if the derivative fprime of func is provided, otherwise the secant method is used. If the second order derivative fprime2 of func is also provided, then Halley’s method is used.

This is crucial, because the secant method on its own does not have very good convergence. With no derivatives (the secant method), the three test cases run like

nfev=44 niter=43
nfev=9 niter=8
nfev=10 niter=9

With one derivative (Newton-Raphson):

nfev=26 niter=13
nfev=12 niter=6
nfev=14 niter=7

With two derivatives (Halley):

nfev=24 niter=8
nfev=12 niter=4
nfev=15 niter=5

via

def dnpv_dr(rate: float, dates: np.ndarray, flows: np.ndarray) -> float:
 years = (dates - dates.min())/FISCAL_YEAR
 return -(
 flows*years*(1 + rate)**(-1 - years)
 ).sum()
def d2npv_dr2(rate: float, dates: np.ndarray, flows: np.ndarray) -> float:
 years = (dates - dates.min())/FISCAL_YEAR
 return (
 flows*years*(years + 1)*(rate + 1)**(-2 - years)
 ).sum()
def xirr(dates: np.ndarray, flows: np.ndarray, guess: float = 0.1) -> float:
 irr, result = optimize.newton(
 func=xnpv, args=(dates, flows), fprime=dnpv_dr, fprime2=d2npv_dr2,
 x0=guess, maxiter=100, full_output=True,
 )
 print(f'nfev={result.function_calls} niter={result.iterations}')

But those are all with "cheating" estimates of -0.15 -> -0.251, -0.1 -> -0.024 and -0.3 -> -0.334 already known to converge. All three of the Newton-based algorithms are very sensitive to the estimate, and replacing the wrapper interface with root_scalar() does not help so long as one of the three methods above is still used.

In this situation, needing an estimate at all is a problem we can eliminate. Examine the table in Scalar functions. Alefeld, Potra, Shi (1995) Algorithm 748, ACM Transactions on Mathematical Software 21(3) - casually toms748 - is a very efficient scalar root algorithm that doesn't require any derivatives or estimates, only a bracket. We could easily assert that all practical IRRs are in the bracket [-99%, 99%]. Give this bracket to toms748, and it has guaranteed convergence to a simple root with a convergence factor of 2.7 per iteration. Again for the three test cases, it takes iteration counts of:

toms748 converged 7
toms748 converged 5
toms748 converged 7

Even though one end of that bracket will always be very, very poorly scaled (function order of magnitude 10^33 in this case), 748 has no problem converging. Even if you reduce xtol (the absolute tolerance) to 1e-15 it typically incurs only one or two more iterations.

Vectorisation

As you've seen in bits and pieces above: vectorise with Numpy. You already had access to that because you'd been using scipy.optimize for which numpy is a dependency. This process is easy because your decision variable is a scalar; only dates and cashflows are vectors.

Other than the obvious - using np.array - you should entirely remove datetime.date and replace it with datetime64.

Observe that the NPV is really a dot product of the cash flow vector and (rate + 1) raised to a negative years vector. This allows for a more efficient root function:

cashflows.dot(rate1**years)

Testing

Write unit and regression tests. I used the three scenarios (including the commented ones) in your code, and cross-checked the results against the two functions you've implemented with Google Sheets.

Prior to vectorisation, this may look like

import datetime
import typing
import numpy as np
from scipy import optimize
def xnpv(rate: float, cashflows: typing.Iterable[tuple[datetime.date, float]]) -> float:
 t0 = min(cashflows, key=lambda t: t[0])[0]
 return sum(
 cf/(1 + rate)**((t - t0).days/365)
 for t, cf in cashflows
 )
def xirr(cashflows: typing.Iterable[tuple[datetime.date, float]], guess: float = 0.1):
 return optimize.newton(
 func=xnpv, args=(cashflows,), x0=guess, maxiter=100,
 )
def test() -> None:
 # OP scenario 1
 cftest = [
 (datetime.date(2001, 12, 5), -2000),
 (datetime.date(2007, 12, 5), -10),
 (datetime.date(2017, 12, 5), 20),
 ]
 assert np.isclose(-1998.30, xnpv(0.05, cftest), atol=0, rtol=1e-4)
 assert np.isclose(-0.251_269, xirr(cftest, -0.15), atol=0, rtol=1e-4)
 # OP scenario 2
 cftest = [
 (datetime.date(2001, 12, 5), -2000),
 (datetime.date(2000, 12, 5), -1000),
 (datetime.date(2017, 12, 5), 2000),
 ]
 assert np.isclose(-2032.634, xnpv(0.05, cftest), atol=0, rtol=1e-4)
 assert np.isclose(-0.024_507, xirr(cftest, -0.1), atol=0, rtol=1e-4)
 # OP scenario 3
 cftest = [
 (datetime.date(2001, 12, 5), -2000),
 (datetime.date(2007, 12, 5), -1000),
 (datetime.date(2017, 12, 5), 20),
 ]
 assert np.isclose(-2736.958, xnpv(0.05, cftest), atol=0, rtol=1e-4)
 assert np.isclose(-0.334_289, xirr(cftest, -0.3), atol=0, rtol=1e-4)
if __name__ == '__main__':
 test()

After vectorisation, including the replacement root-finding approach,

import numpy as np
from scipy import optimize
# https://en.wikipedia.org/wiki/Annual_percentage_rate
# APR is commonly defined with a 365-day year, even though that isn't how time works
FISCAL_YEAR = np.timedelta64(365, 'D')
def xnpv(rate: float, dates: np.ndarray, cashflows: np.ndarray) -> float:
 return cashflows.dot(
 (1 + rate)**((dates.min() - dates)/FISCAL_YEAR)
 )
def npv_error(rate1: float, years: np.ndarray, cashflows: np.ndarray) -> float:
 return cashflows.dot(rate1**years)
def xirr(
 dates: np.ndarray, cashflows: np.ndarray,
 lb: float=-0.9, ub: float=0.9, tol: float = 1e-15,
) -> float:
 years = (dates.min() - dates)/FISCAL_YEAR
 result = optimize.root_scalar(
 f=npv_error, args=(years, cashflows), method='toms748', 
 bracket=(1+lb, 1+ub), xtol=tol,
 )
 if not result.converged:
 raise ValueError(result.flag)
 return result.root - 1
def test_irr(
 dates: np.ndarray, cashflows: np.ndarray,
 npv: float, irr: float, npv_rate: float = 0.05,
) -> None:
 npv_actual = xnpv(rate=npv_rate, dates=dates, cashflows=cashflows)
 assert np.isclose(npv, npv_actual, atol=0, rtol=1e-4)
 irr_actual = xirr(dates=dates, cashflows=cashflows)
 assert np.isclose(irr, irr_actual, atol=0, rtol=1e-4)
 root_error = xnpv(irr_actual, dates, cashflows)
 assert np.isclose(0., root_error, rtol=0, atol=1e-9)
def test() -> None:
 # OP scenario 1
 dates = np.array(('2001-12-05', '2007-12-05', '2017-12-05'), dtype='datetime64[D]')
 cashflows = np.array((-2000, -10, 20))
 test_irr(dates, cashflows, npv=-1998.30, irr=-0.251_269)
 # OP scenario 2
 dates = np.array(('2001-12-05', '2000-12-05', '2017-12-05'), dtype='datetime64[D]')
 cashflows = np.array((-2000, -1000, 2000))
 test_irr(dates, cashflows, npv=-2032.634, irr=-0.024_507)
 # OP scenario 3
 dates = np.array(('2001-12-05', '2007-12-05', '2017-12-05'), dtype='datetime64[D]')
 cashflows = np.array((-2000, -1000, 20))
 test_irr(dates, cashflows, npv=-2736.958, irr=-0.334_289)
if __name__ == '__main__':
 test()