Fastest possible Cython for Black-Scholes algorithm

Question

I started with a pure python implementation, and have been trying to get the performance as close to native C as possible using numpy, numexpr, and cython. Here is the the numpy version that I compile with cython:

import numpy as np
cimport cython
cimport numpy as np
from cython.parallel cimport prange, parallel
from scipy.special import erf as sp_erf
from libc.math cimport log, exp, sqrt, erf
from libc.stdlib cimport malloc, free

DTYPE = np.float64
ctypedef np.float64_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
def black_scholes_cython(int nopt,
                         np.ndarray[DTYPE_t, ndim=1] price,
                         np.ndarray[DTYPE_t, ndim=1] strike,
                         np.ndarray[DTYPE_t, ndim=1] t,
                         double rate,
                         double vol):


    cdef np.ndarray d1 = np.zeros(nopt, dtype=DTYPE)
    cdef np.ndarray d2 = np.zeros(nopt, dtype=DTYPE)
    cdef np.ndarray call = np.zeros(nopt, dtype=DTYPE)
    cdef np.ndarray put = np.zeros(nopt, dtype=DTYPE)

    d1 = (np.log(price / strike) + (rate + vol * vol / 2.) * t) / (vol * np.sqrt(t))
    d2 = (np.log(price / strike) + (rate - vol * vol / 2.) * t) / (vol * np.sqrt(t))

    cdef np.ndarray n_d1 = 0.5 + 0.5 * sp_erf(d1 / np.sqrt(2))
    cdef np.ndarray n_d2 = 0.5 + 0.5 * sp_erf(d2 / np.sqrt(2))
    cdef np.ndarray neg_d1 = np.negative(n_d1)
    cdef np.ndarray neg_d2 = np.negative(n_d1)
    cdef np.ndarray e_rt = np.exp(-rate * t)

    call = price * n_d1 - e_rt * strike * n_d2
    put = e_rt * strike * neg_d2 - price * neg_d1

    return call, put

This code takes about 56 seconds to compute 4,194,304 options. C code with MKL takes about 7 seconds. I wanted to try to take advantage of parallelism, but I read here that you have to remove all python objects from a block to run the code in parallel. So I tried rewriting the function like this:

@cython.boundscheck(False)
@cython.wraparound(False)
def black_scholes_cython_parallel(int nopt,
                                  double[:] price,
                                  double[:] strike,
                                  double[:] t,
                                  double rate,
                                  double vol,
                                  bint ret_call=1):

    cdef double[:] call = np.zeros(nopt, dtype=DTYPE)
    cdef double[:] put = np.zeros(nopt, dtype=DTYPE)

    cdef double *d1 = <double *>malloc(nopt * sizeof(double))
    cdef double *d2 = <double *>malloc(nopt * sizeof(double))

    cdef int i

    cdef DTYPE_t n_d1
    cdef DTYPE_t n_d2
    cdef DTYPE_t neg_d1
    cdef DTYPE_t neg_d2
    cdef DTYPE_t s
    cdef DTYPE_t p
    cdef DTYPE_t e_rt

    with nogil, parallel():
        for i in prange(nopt, schedule='guided'):
            d1[i] = (log(price[i] / strike[i]) + (rate + vol * vol / 2.) * t[i]) / (vol * sqrt(t[i]))
            d2[i] = (log(price[i] / strike[i]) + (rate - vol * vol / 2.) * t[i]) / (vol * sqrt(t[i]))
            n_d1 = 0.5 + 0.5 * erf(d1[i] / sqrt(2))
            n_d2 = 0.5 + 0.5 * erf(d2[i] / sqrt(2))
            e_rt = exp(-rate * t[i])

            neg_d1 = -n_d1
            neg_d2 = -n_d1
            s = strike[i]
            p = price[i]

            call[i] = p * n_d1 - e_rt * s * n_d2
            put[i] = e_rt * s * neg_d2 - p * neg_d1

    free(d1)
    free(d2)

    return call if ret_call else put

This compiles and runs, but I don't see any real speedup over the numpy version. How can I take advantage of parallelism to make this code run as fast as possible?

Edit Here is the testing code:

import numba as nb
import numexpr as ne
import numpy as np
import time

from math import log, sqrt, exp
from random import seed, uniform
from scipy.stats import norm
from scipy.special import erf

SEED = 7777777
S0L = 10.0
S0H = 50.0
XL = 10.0
XH = 50.0
TL = 1.0
TH = 2.0
RISK_FREE = 0.1
VOLATILITY = 0.2

def gen_data(nopt):
    seed(SEED)
    price = []
    strike = []
    time = []

    for i in range(0, nopt):
        price.append(uniform(S0L, S0H))
        strike.append(uniform(XL, XH))
        time.append(uniform(TL, TH))
    return price, strike, time

if __name__ == '__main__':
    size = 4 * 1024 * 1024
    n = 100
    nthreads = 2
    iterations = range(n)

    def run(alg, nopt=size, nparr=False):
        for i in range(0, 4):
            price, strike, t = gen_data(nopt)
            if nparr:
                price = np.array(price, dtype=np.float64)
                strike = np.array(strike, dtype=np.float64)
                t = np.array(t, dtype=np.float64)
            t0 = time.clock()
            for j in iterations:
                alg(nopt, price, strike, t, RISK_FREE, VOLATILITY)
            t1 = time.clock()
            exec_time = t1 - t0
            print("Size: {}\nTime: {}".format(nopt, exec_time))
            nopt = nopt * 2

    import os
    if os.path.exists('./bs.so'):
        import bs
        print("Cython")
        run(bs.black_scholes_cython, nparr=True)
        print("Parallel Cython")
        run(bs.black_scholes_cython_parallel, nparr=True)
    else:
        print("Must run 'python setup.py build_ext --inplace' first.")

setup.py is just setup(ext_modules=cythonize("bs.pyx"),).

Edit: I added the -openmp compiler flag and now it goes much faster. However, something is wrong. When I run the code, the timing reports that it took 111 seconds, but if I use wall time while it's running, it completes in about 12 seconds. Any ideas about why there's a discrepancy?

Answer 1

Well it's no surprise that the first function doesn't perform differently than the Python version, as it will still call into Python and to NumPy - nothing to be gained by doing that via Cython.

The setup.py is probably something closer to this:

from distutils.core import setup

from Cython.Build import cythonize

setup(ext_modules=cythonize("bs.pyx"),)

In test.py the generated .so file isn't necessarily named that way, I'd rather try the import and catch the exception, that's a more fool proof way:

try:
    import bs
except ImportError:
    print("Must run 'python setup.py build_ext --inplace' first.")
else:
    print("Cython")
    run(bs.black_scholes_cython, nparr=True)
    print("Parallel Cython")
    run(bs.black_scholes_cython_parallel, nparr=True)

• Possible bug: n_d2/neg_d2 should probably also use d2 instead of d1?
• Memory should be free'd in any case, so the main loop should be protected with try/finally.
• I think it's more customary to use time.time instead of time.clock.

Now, I got some better timings by manually decomposing the terms in the parallel version. Next remove the d1/d2 arrays, as the values stored in there aren't actually used again; that also gets rid of the memory allocation.

Finally, if only one of those values is going to be returned, only compute one of them; assuming that the if isn't much slower than the double calculation it might be worth it to put a conditional there.

Assuming I didn't mess up the rewrite, I have now the following code, which is a bit faster than the initial version:

@cython.boundscheck(False)
@cython.wraparound(False)
def black_scholes_cython_parallel2(int nopt,
                                   double[:] price,
                                   double[:] strike,
                                   double[:] t,
                                   double rate,
                                   double vol,
                                   bint ret_call=1):

    cdef double[:] ret = np.zeros(nopt, dtype=DTYPE)

    cdef int i

    cdef DTYPE_t d1, d2, n_d1, n_d2, s, p, e_rt, v, x, y

    with nogil, parallel():
        for i in prange(nopt, schedule='guided'):
            s = strike[i]
            p = price[i]
            v = vol * sqrt(t[i])
            x = log(p / s) + rate * t[i]
            y = vol * vol / 2. * t[i]

            d1 = (x + y) / v
            d2 = (x - y) / v

            n_d1 = 0.5 + 0.5 * erf(d1 / sqrt(2))
            n_d2 = 0.5 + 0.5 * erf(d2 / sqrt(2))

            e_rt = exp(-rate * t[i])

            if ret_call:
                ret[i] = p * n_d1 - e_rt * s * n_d2
            else:
                ret[i] = e_rt * s * -n_d2 - p * -n_d1

    return ret

For building the extension, check that optimisations are on, i.e. -O2 or -O3 is passed as arguments during compilation; there might be other interesting flags as well.

Looking at the code samples from the Intel page there's probably only so much you can do without getting into vectorised operations.

Answer 2

Not sure why your Cython code is so slow. In native Python using numpy and scipy.stats import norm you can easily price 10M options in 3s, just vectorize it. It is 3 lines of code! d1=... d2=... result = ... just pass the call/put as a 1 or -1 and it will be this compact.