Why Python?

Python has a reputation primarily as a scripting language, functioning as "glue" between other codebases. The QuantStart team, however, believe that with tools such as Numpy and SciPy it is perfectly capable of being utilised to price options. However, it is yet another language to learn - so why should you invest the time? We have listed the primary benefits below:

• For the same task, Python usually requires far less coding than, say, Java or C++.
• Python contains some great libraries out of the box - the so called "batteries included" philosophy.
• Python is straightforward to grasp; some people believe that it is almost like reading pseudo-code.
• By integrating tools such as NumPy, Numba and Cython it is becoming faster all the time.

The articles to follow will concentrate on simplicity, rather than optimisation. We will follow Daniel Duffy's philosophy of "First we get it working, then we get it right, then we optimise it". To this end, we will initially be focusing on clear and explicit syntax. It should be obvious how we have gone from an algorithm to its implementation, with the minimum of head-scratching. After we have code working, we can make sure it is producing the correct values. Finally, we can optimise it and put it into production.

Black-Scholes Pricing Formula

It has already been outlined that the reader is to be familiar with the Black-Scholes formula for the pricing of European Vanilla Calls and Puts. For completeness, the price of a European Vanilla Call, \(C(S,t)\) is given below, where \(S\) is the underlying asset, \(K\) is the strike price, \(r\) is the risk-free rate, \( T \) is the time to maturity and \(\sigma\) is the (constant) volatility of the underlying \(S\) (\(N\) is described below):

\[ C(S,t) = SN(d_1) - Ke^{-rT} N(d_2) \]

With \(d_1\) and \(d_2\) defined as follows:

\[ d_1 = \frac{log(S/K) + (r+\frac{\sigma^2}{2})T}{\sigma \sqrt{T} } \]

\[ d_2 = d_1 - \sigma \sqrt{T} \]

Thanks to Put-Call Parity, we are also able to price a European Vanilla Put \(P(S,t)\) with the following formula:

\[ P(S,t) = Ke^{-rT} - S + C(S,t) = Ke^{-rT} - S + (SN(d_1) - Ke^{-rT} N(d_2)) \]

The remaining function we have yet to describe is \(N\). This is the cumulative distribution function of the standard normal distribution. The formula for \( N \) is given by:

\[ N(x) = \frac{1}{\sqrt{2 \pi}} \int^x_{-\infty} e^{-t^{2} /2} dt \]

In addition, we would like to have closed form solutions for the "Greeks", which are the option price sensitivities to the underlying variables and parameters. For this reason we also need the formula for the probability density function of the standard normal distribution which is given below:

\[ f(x) = \frac{1}{\sqrt{2 \pi}} e^{-x^{2}/2}\]

Our task is now to utilise Python to implement these functions and provide us with values for the closed-form solution to the price of a European Vanilla Call or Put with their associated sensitivities.

Python Implementation of Statistical Functions

Open a new Python file in your favourite IDE and name it statistics.py. We require some mathematical functions so we will use the math library and import the natural logarithm, the exponential function and the mathematical constant \(\pi\):

from math import exp, log, pi

The first function to define is the standard normal distribution probability density function:

def norm_pdf(x):
    """
    Standard normal probability density function
    """
    return (1.0/((2*pi)**0.5))*exp(-0.5*x*x)

This function is fairly self-explanatory. The next function to code up is the CDF of the normal distribution. However, it contains an integral with \(-\infty\) as one of its limits - so how can we appoximate this on the computer? The method utilised (26.2.17) comes from the famous Abramovitz & Stegun textbook on mathematical functions. The Wikipedia article on the Normal Distribution sheds more light on this and other methods. Here is the Python listing for the algorithm:

def norm_cdf(x):
    """
    An approximation to the cumulative distribution
    function for the standard normal distribution:
    N(x) = \frac{1}{sqrt(2*\pi)} \int^x_{-\infty} e^{-\frac{1}{2}s^2} ds
    """
    k = 1.0/(1.0+0.2316419*x)
    k_sum = k * (0.319381530 + k * (-0.356563782 + \
        k * (1.781477937 + k * (-1.821255978 + 1.330274429 * k))))

    if x >= 0.0:
        return (1.0 - (1.0 / ((2 * pi)**0.5)) * exp(-0.5 * x * x) * k_sum)
    else:
        return 1.0 - norm_cdf(-x)

Notice that this function calls itself if $x$ is negative.

We now have the two statistical functions necessary for calculating the closed-form options prices for European vanilla calls and puts.

Python Implementation of Closed-Form European Vanilla Call-Put Prices

We need to create a second file, which we will call closed_form.py. It should reside in the same file directory as the statistics.py file in order for the following import syntax to work correctly:

from math import exp, log, sqrt
from statistics import norm_pdf, norm_cfd

Referring back to the Black-Scholes formula for the price of a call, \(C(S,t)\), we can see that it requires two separate functions named \(d_j\), \( j \in \{1,2\} \). The \(d_j\) take \(S\), \(T\), \(K\), \(r\) and \(\sigma\) (v in the code) as arguments. It is possible to write the \(d_j\) as a single function which takes an additional argument \(j\). The code for this function is provided below:

def d_j(j, S, K, r, v, T):
    """
    d_j = \frac{log(\frac{S}{K})+(r+(-1)^{j-1} \frac{1}{2}v^2)T}{v sqrt(T)}
    """
    return (log(S/K) + (r + ((-1)**(j-1))*0.5*v*v)*T)/(v*(T**0.5))

All that remains now is to code up the functions for \(C(S,t)\) and \(P(S,t)\). The Python code for \(C(S,t)\) is provided below:

def vanilla_call_price(S, K, r, v, T):
    """
    Price of a European call option struck at K, with
    spot S, constant rate r, constant vol v (over the
    life of the option) and time to maturity T
    """
    return S * norm_cdf(d_j(1, S, K, r, v, T)) - \
        K*exp(-r*T) * norm_cdf(d_j(2, S, K, r, v, T))

Finally, the code for \(P(S,t)\):

def vanilla_put_price(S, K, r, v, T):
    """
    Price of a European put option struck at K, with
    spot S, constant rate r, constant vol v (over the
    life of the option) and time to maturity T
    """
    return -S * norm_cdf(-d_j(1, S, K, r, v, T)) + \
        K*exp(-r*T) * norm_cdf(-d_j(2, S, K, r, v, T))

This concludes the coding of formulae for the statistical distribution functions in statistics.py and the vanilla call option prices in closed_form.py. At this stage it is prudent to check that the formulae produce the correct results and satisfy the known bounds on the prices (such as Put-Call Parity). That will be the subject of a future article.