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 !
