A Try of Binomial Tree for Option Valuation¶

In [246]:
import numpy as np
import pandas as pd
In [247]:
#timing
from time import time
from functools import wraps

def timing(f):
    @wraps(f)
    def wrap(*args, **kw):
        ts=time()
        result=f(*args, **kw)
        te=time()
        print('func:%r args:[%r, %r] took: %2.4f sec'% \
             (f.__name__, args, kw, te-ts))
        return result
    return wrap

Initalise Parameters¶

In [248]:
S0 = 100
K = 100
T = 1
r = 0.06
N = 100
u = 1.1
d = 1/u
opttype = 'C' # option type is call or put. C or P

1. A Slow Method¶

In [249]:
@timing
def binomial_slow(S0, K, T, r, N, u, d, opttype):
    # 0. Pre-compute Constants
    dt = T/N
    prob_u = (np.exp(r*dt)-d)/(u-d)
    df = np.exp(-r*dt)
    
    # 1. initialise asset price at maturity - Time step = N
    S = np.zeros(N+1) # create a N+1 empty array 
    S[0] = S0 * d**N # the very bottom the right column
    for j in range(1,N+1):
        S[j]=S[j-1]*u/d
        
    # 2. initialise the option value at maturity
    C=np.zeros(N+1)
    for j in range(0, N+1):
        C[j]=max(0, S[j]-K )
    
    # 3. Step backward through trees
    for i in np.arange(N,0,-1): #work backward (Inital, End, Step)
        for j in range(0, i):
            C[j] = df * (C[j]*(1-prob_u) + C[j+1]*prob_u)
        #Collapse to the bottom until arriving C[0]
    return C[0]


binomial_slow(S0, K, T, r, N, u, d, opttype)

func:'binomial_slow' args:[(100, 100, 1, 0.06, 100, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0050 sec
Out[249]:
38.45357509478931

2. A Fast Way¶

In [250]:
@timing
def binomial_fast(S0, K, T, r, N, u, d, opttype):
    # 0. Pre-compute Constants
    dt = T/N
    prob_u = (np.exp(r*dt)-d)/(u-d)
    df = np.exp(-r*dt)
    
    # The time-saving part is we apply np.array instead of looping.
    
    # 1. We directly apply, ST = S0 * u^i * d^(T-i)
    S = S0 * d**(np.arange(N,-1,-1)) * u**np.arange(0,N+1,1)
    
    # 2. We also use np.array to find the max of option in the tree
    # c = max (S-K, 0)   in an array
    C=np.maximum(S - K, np.zeros(N+1))
    
    # 3. Step backward through trees
    for i in np.arange(N,0,-1):
        C = df * (prob_u*C[1:i+1] + (1-prob_u)*C[0:i])
    return C


binomial_fast(S0, K, T, r, N, u, d, opttype)

func:'binomial_fast' args:[(100, 100, 1, 0.06, 100, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0010 sec
Out[250]:
array([38.45357509])

3. Comparison of Time¶

In [116]:
for N in [3,5,10,100,5000]:
    binomial_slow(S0, K, T, r, N, u, d, opttype)
    binomial_fast(S0, K, T, r, N, u, d, opttype)

func:'binomial_slow' args:[(100, 100, 1, 0.06, 3, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0001 sec
func:'binomial_fast' args:[(100, 100, 1, 0.06, 3, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0001 sec
func:'binomial_slow' args:[(100, 100, 1, 0.06, 5, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0001 sec
func:'binomial_fast' args:[(100, 100, 1, 0.06, 5, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0003 sec
func:'binomial_slow' args:[(100, 100, 1, 0.06, 10, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0001 sec
func:'binomial_fast' args:[(100, 100, 1, 0.06, 10, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0005 sec
func:'binomial_slow' args:[(100, 100, 1, 0.06, 100, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0047 sec
func:'binomial_fast' args:[(100, 100, 1, 0.06, 100, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0010 sec
func:'binomial_slow' args:[(100, 100, 1, 0.06, 5000, 1.1, 0.9090909090909091, 'C'), {}] took: 7.5409 sec
func:'binomial_fast' args:[(100, 100, 1, 0.06, 5000, 1.1, 0.9090909090909091, 'C'), {}] took: 0.0453 sec

4. Volatility included¶

$$ p u+(1-p)v = e^{\mu\ \delta t}$$$$ p=\frac{e^{\mu\ \delta t}-v}{u-v} $$
$$ u\approx 1+\sigma \sqrt{\delta t} +\frac{1}{2}\sigma^2 \delta t$$$$ v\approx 1-\sigma \sqrt{\delta t} +\frac{1}{2}\sigma^2 \delta t$$

and, $$ p\approx \frac{1}{2}+\frac{(\mu-\frac{1}{2}\sigma^2)\sqrt{\delta t}}{2\sigma} $$

In [214]:
S0=100
N=4
T=1/3
r=0.1
sigma=0.2
In [350]:
def binomial_fast_volatility(S0, K, T, r, N, sigma, opttype):
    # 0. Pre-compute Constants
    # P.S. Different Methods are applied to get those constants, but similar logic
    dt = T/N
    df = np.exp(-r*dt)
    temp1 = np.exp((r + sigma**2)*dt)
    temp2 = 0.5 * (df + temp1)
    u = temp2 + np.sqrt(temp2**2 - 1)
    d = 1/u
    prob_u = ( np.exp(r*dt) - d )/( u - d)
    
    # The time-saving part is we apply np.array instead of looping.
    
    # 1. We directly apply, ST = S0 * u^i * d^(T-i)
    S = S0 * d**(np.arange(N,-1,-1)) * u**np.arange(0,N+1,1)
    
    # 2. We also use np.array to find the max of option in the tree
    # c = max (S-K, 0)   in an array
    C=np.maximum(S - K, np.zeros(N+1))
    
    # 3. Step backward through trees
    for i in np.arange(N,0,-1):
        C = df * (prob_u*C[1:i+1] + (1-prob_u)*C[0:i])
    return C

binomial_fast_volatility(S0, K, T, r, N, sigma, opttype)
Out[350]:
array([6.13447247])
In [254]:
U=50
a=np.zeros(U)
for i in range(1,U+1):
    a[i-1]=binomial_fast_volatility(S0, K, T, r, i, sigma, opttype)
pd.Series(a).plot()
Out[254]:
<AxesSubplot:>

5. American Style¶

In [263]:
S0=100
K=100
N=4
T=1/3
r=0.1
sigma=0.2
In [359]:
def binomial_fast_volatility_American(S0, K, T, r, N, sigma, opttype):
    # 0. Pre-compute Constants
    # P.S. Different Methods are applied to get those constants, but similar logic
    dt = T/N
    df = np.exp(-r*dt)
    temp1 = np.exp((r + sigma**2)*dt)
    temp2 = 0.5 * (df + temp1)
    u = temp2 + np.sqrt(temp2**2 - 1)
    d = 1/u
    prob_u = ( np.exp(r*dt) - d )/( u - d)
    
    # The time-saving part is we apply np.array instead of looping.
    
    # 1. We directly apply, ST = S0 * u^i * d^(T-i)
    S = S0 * d**(np.arange(N,-1,-1)) * u**np.arange(0,N+1,1)
    
    # 2. We also use np.array to find the max of option in the tree
    # c = max (S-K, 0)   in an array
    C=np.maximum(S - K, np.zeros(N+1)) #len(C)=N+1
    
    # 3. Step backward through trees
    for i in np.arange(N,0,-1):
        C = df * (prob_u*C[1:i+1] + (1-prob_u)*C[0:i])
        
        # for american
        #S = S0 * d**(np.arange(i-1,-1,-1)) * u**np.arange(0,i,1) 
        # the above row has same meaning as the row below, but much more time consuming
        payoff = S[:i]/d - K
        C=np.maximum(C, payoff) # the second term is the payoff
    return C



binomial_fast_volatility_American(S0, K, T, r, N, sigma, opttype)
Out[359]:
array([13.26244768])
In [360]:
S0=100
K=100
N=120
T=1
r=0.1
sigma=0.2

print(binomial_fast_volatility(S0, K, T, r, N, sigma, opttype),binomial_fast_volatility_American(S0, K, T, r, N, sigma, opttype))

[13.26244768] [13.26244768]

by Fanyu Zhao