For the last LFPUG meeting, I have developed a simple example showing how Cython could be used to improve the performance of Python. The idea here is not to find the best way to do it but just to try to get a great speed up without digging into the hard stuff.

1. The model

The theoretical binomial model can be used to value options. “Paul Wilmott’s introduces Quantitative Finance” shows a VB example of a call valuation with such model. The following code is inspired by the VB code and gives you a pure Python version of it.

''' Call option valuation using the theoretical binomial model. Inspired from Wilmott, Paul, "Paul Wilmott Introduces Quantitative Finance", 2nd ed., John Wiley & Sons, 2007 ''' from __future__ import division from numpy import array, zeros, exp, sqrt, linspace, zeros_like def payoff(S, K): return (S-K) if S > K else 0 def binomial_price(asset, volatility, interest_rate, strike, expiry, steps): """ """ s = zeros((steps+1)) # asset array v = zeros((steps+1)) # option array time_step = (1.0 * expiry) / steps discount_factor = exp(-interest_rate * time_step) temp1 = exp((interest_rate + volatility * volatility) * time_step) temp2 = 0.5 * (discount_factor + temp1) u = temp2 + sqrt(temp2 * temp2 -1) d = 1/u p = (exp(interest_rate * time_step) - d ) / (u-d) s[0] = asset for i in xrange(1, steps+1): for j in xrange(i, 0, -1): s[j] = u * s[j-1] s[0] = d * s[0] for j in xrange(steps+1): v[j] = payoff(s[j], strike) for n in xrange(steps): for j in xrange(steps): v[j] = ((p * v[j+1]) + ((1-p) * v[j])) * discount_factor return v[0]

You can execute it for some volatility values and get a plot of it with execution time.

if __name__ == '__main__': import time from pylab import plot, show v_rates = linspace(.01, 0.2, 1000) option_prices = zeros_like(v_rates) asset = 100 strike = 100 interest_rate = .05 expiry = 4 steps = 12 t = time.time() for i, volatility in enumerate(v_rates): option_prices[i] = binomial_price(asset, volatility, interest_rate, strike, expiry, steps) print 'Execution time:', time.time() - t plot(v_rates, option_prices) show()

It gives you an execution time of around 0.7 seconds and the following results

2. Cython improvements

Cython can then used on the original file renamed to pyx by :

- typedef’ing as much as possible into the file

- inlining some function

- using function from the math library in place of the Python function

This gives the following Cython file :

import numpy cimport cython cimport numpy as np cdef extern from "math.h": double sqrt(double x) double exp(double x) cdef inline double payoff(double S, double K): return (S-K) if S > K else 0. @cython.boundscheck(False) @cython.wraparound(False) def binomial_price(double asset, double volatility, double interest_rate, double strike, unsigned int expiry, unsigned int steps): """ Compute the option price for the given arguments using a binomial model """ cdef Py_ssize_t bsteps = steps + 1 cdef np.ndarray[double, ndim=1] s = numpy.zeros(dtype="d", shape=(bsteps)) # asset array cdef np.ndarray[double, ndim=1] v = numpy.zeros(dtype="d", shape=(bsteps)) # option array cdef double time_step = 1.0 * expiry / bsteps cdef double discount_factor = exp(-interest_rate * time_step) cdef double temp1 = exp((interest_rate + volatility * volatility) * time_step) cdef double temp2 = 0.5 * (discount_factor + temp1) cdef double u = temp2 + sqrt(temp2 * temp2 -1) cdef double d = 1./u cdef double p = (exp(interest_rate * time_step) - d ) / (u-d) s[0] = asset cdef Py_ssize_t i, j, n for i in range(1, bsteps): for j in range(i, 0, -1): s[j] = u * s[j-1] s[0] = d * s[0] for j in range(bsteps): v[j] = payoff(s[j], strike) for n in range(steps): for j in range(steps): v[j] = ((p * v[j+1]) + ((1-p) * v[j])) * discount_factor return v[0]

All of that seems not that hard if you do not know about interfacing Python with C or even C.

The execution time for the Cython version of the file is 0.009 seconds using the same code as the one producing the plot for the pure Python version. That is 77 times faster. Quiet nice !