Knowledge Base/Finance/Option Pricing in Python and Simple English

From Thalesians
Jump to: navigation, search

Under construction.

Introduction

Option pricing does not need to be complicated. At least the basics of it are easily comprehensible. The goal of this tutorial is not to prove the mathematics behind option pricing. Rather it is to show how to implement it in a simple and natural programming language — Python.

The interested reader can then refer to the relevant papers to fill in the details.

Models and methods

When talking about option pricing one should distinguish between models and methods. A model is a simplified statement of our view of the dynamics and randomness that we believe determine or drive the price of a particular asset. Having postulated the model, we then use one of the several methods from our arsenal to determine the price of an asset given this particular model.

The commonly used models include Black-Scholes, local volatility, stochastic volatility, local-stochastic volatility mixtures. The more refined the model, the better it represents the world, and therefore it should give us a more realistic price.

However, the more complex the model, the less mathematically tractable it is, and we are likely to be more limited in our choice of available methods.

The methods include analytic and numerical procedures, such as trees, finite difference methods and Monte Carlo simulations. The analytic method consists in using mathematical analysis to derive the pricing formulae from the mathematical statement of our model. The resulting formulae can be used immediately and are computationally cheap. Just plug in the numbers, and you get the price back! Such formulae are known as closed form solutions. For more complex models the analytic solutions are often very difficult to find.

The numerical procedures tend to be more computationally expensive because we have to iterate. We use them when we can't find analytic solutions. Their advantage is that they make many difficult problems tractable. We may not be able to find the closed form solutions but at least we can obtain the prices to a reasonable degree of accuracy. 

import math
import numpy
import matplotlib.pyplot
import scipy.stats
 
# See [Hull], Chapter 13. The pricing formulae are given in Section 13.8
 
def d1(S0, K, r, sigma, T):
    return (math.log(S0 / K) + (r + sigma**2 / 2) * T) / (sigma * math.sqrt(T))
 
def d2(S0, K, r, sigma, T):
    return (math.log(S0 / K) + (r - sigma**2 / 2) * T) / (sigma * math.sqrt(T))
 
def BlackScholesEuropeanCallPrice(S0, K, r, sigma, T):
    return S0 * scipy.stats.norm.cdf(d1(S0, K, r, sigma, T)) - K * math.exp(-r * T) * scipy.stats.norm.cdf(d2(S0, K, r, sigma, T))
 
def BlackScholesEuropeanPutPrice(S0, K, r, sigma, T):
    return K * math.exp(-r * T) * scipy.stats.norm.cdf(-d2(S0, K, r, sigma, T)) - S0 * scipy.stats.norm.cdf(-d1(S0, K, r, sigma, T))
 
S0 = 100.0
K = 110.0
r = 0.03
sigma = 0.1
T = 0.5
 
print "S0\tstock price at time 0:", S0
print "K\tstrike price:", K
print "r\tcontinuously compounded risk-free rate:", r
print "sigma\tvolatility of the stock price per year:", sigma
print "T\ttime to maturity in trading years:", T
# This is usually Number of trading days until option maturity / 252, see
# [Hull], Section 13.4
 
c_BS = BlackScholesEuropeanCallPrice(S0, K, r, sigma, T)
p_BS = BlackScholesEuropeanPutPrice(S0, K, r, sigma, T)
 
print
print "c_BS\tBlack-Scholes European call price:", c_BS
print "p_BS\tBlack-Scholes European put price:", p_BS
 
"""
S0s = numpy.arange(90.0, 120.0, 0.5)
cs_1D = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=1.0 / 252.0) for S0 in S0s]
cs_3M = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=0.25) for S0 in S0s]
cs_6M = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=0.5) for S0 in S0s]
cs_1Y = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=1.0) for S0 in S0s]
 
matplotlib.pyplot.plot(S0s, cs_1D, 'r-', lw=2)
matplotlib.pyplot.plot(S0s, cs_3M, 'y-', lw=2)
matplotlib.pyplot.plot(S0s, cs_6M, 'b-', lw=2)
matplotlib.pyplot.plot(S0s, cs_1Y, 'g-', lw=2)
matplotlib.pyplot.xlabel("$S_0$")
matplotlib.pyplot.ylabel("$c$")
matplotlib.pyplot.title("European call price versus $S_0$, $T = \\frac{1}{252}, \\frac{1}{4}, \\frac{1}{2}, 1$")
matplotlib.pyplot.grid(True)
matplotlib.pyplot.show()
"""
 
def BinomialTreeEuropeanCallPrice(S0, K, r, sigma, T, N=30):
    deltaT = T / N
 
    u = math.exp(sigma * math.sqrt(deltaT))
    d = 1.0 / u
 
    # Initialise our f_{i,j} tree with zeros
    fs = [[0.0 for j in xrange(i + 1)] for i in xrange(N + 1)]
 
    a = math.exp(r * deltaT)
 
    p = (a - d) / (u - d)
    oneMinusP = 1.0 - p
 
    # Compute the leaves, f_{N, j}
    for j in xrange(i+1):
        fs[N][j] = max(S0 * u**j * d**(N - j) - K, 0.0)
 
    for i in xrange(N-1, -1, -1):
        for j in xrange(i + 1):
            fs[i][j] = math.exp(-r * deltaT) * (p * fs[i + 1][j + 1] + oneMinusP * fs[i + 1][j])
 
    # print fs
 
    return fs[0][0]
 
c_BT = BinomialTreeEuropeanCallPrice(S0, K, r, sigma, T)
print "c_BT\tBinomial tree European call price:", c_BT
 
def MonteCarloEuropeanCallPrice(S0, K, r, sigma, T, N=100, pathCount=5000):
    deltaT = T / N
 
    fs = [0.0 for i in xrange(pathCount)]
 
    for i in xrange(pathCount):
        S = S0
        epsilons = scipy.stats.norm.rvs(size=N)
        for j in xrange(N):
            S += r * S * deltaT + sigma * S * epsilons[j] * math.sqrt(deltaT)
        fs[i] = max(S - K, 0.0)
 
    return scipy.mean(fs)
 
c_MC = MonteCarloEuropeanCallPrice(S0, K, r, sigma, T)
print "c_MC\tMonte Carlo European call price:", c_MC
Retrieved from "http://www.thalesians.com/finance/index.php/Knowledge_Base/Finance/Option_Pricing_in_Python_and_Simple_English"
Personal tools