Improving Python speed using Cython : Binomial option valuation example

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
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 !

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: