Binomial Asset Pricing Model – Choosing Parameters

Binomial Tree Model / By

Implementation of a simple slow and fast binomial pricing model in python. We will treat binomial tree as a network with nodes (i,j) with i representing the time steps and j representing the number of ordered price outcome.
We will be implementing the following for the binomial_tree_slow algorithm, but it will work with the fast algorithm we made previously:

  • Cox, Ross and Rubinstein (CRR)
  • Jarrow and Rudd (JR)
  • Equal probabilities (EQP)
  • Trigeorgis (TRG)

Binomial Tree Representation

Stock tree can be represented using nodes (i,j) and intial stock price \(S_0\)
\(S_{i,j} = S_0u^{j}d^{i-j}\)
\(C_{i,j}\) represents contract price at each node (i,j). Where \(C_{N,j}\) represents final payoff function that we can define.
For this tutorial will will price a European Call, so \(C_{N,j} = max(S_{N,j}-K,0)\)

import numpy as np

# Initialise parameters
S0 = 100      # initial stock price
K = 110       # strike price
T = 0.5       # time to maturity in years
r = 0.06      # annual risk-free rate
N = 100       # number of time steps
sigma = 0.3   # Annualised stock price volatility
opttype = 'C' # Option Type 'C' or 'P'

Cox, Ross and Rubinstein (CRR) Method

Here we choose equal jump sizes

def CRR_method(K,T,S0,r,N,sigma,opttype='C'):
    #precomute constants
    dt = T/N
    u = np.exp(sigma*np.sqrt(dt))
    d = 1/u
    q = (np.exp(r*dt) - d) / (u-d)
    disc = np.exp(-r*dt)
    
    # initialise asset prices at maturity - Time step N
    S = np.zeros(N+1)
    S[0] = S0*d**N
    for j in range(1,N+1):
        S[j] = S[j-1]*u/d
    
    # initialise option values at maturity
    C = np.zeros(N+1)
    for j in range(0,N+1):
        if opttype == 'C':
            C[j] = max(0, S[j]-K)
        else:
            C[j] = max(0, K - S[j])
        
    # step backwards through tree
    for i in np.arange(N,0,-1):
        for j in range(0,i):
            C[j] = disc * ( q*C[j+1] + (1-q)*C[j] )
    
    return C[0]

CRR_method(K,T,S0,r,N,sigma,opttype='C')

Jarrow and Rudd (JR) Method

Here we choose equal risk-neutral probabilities

def JR_method(K,T,S0,r,N,sigma,opttype='C'):
    #precomute constants
    dt = T/N
    nu = r - 0.5*sigma**2
    u = np.exp(nu*dt + sigma*np.sqrt(dt))
    d = np.exp(nu*dt - sigma*np.sqrt(dt))
    q = 0.5
    disc = np.exp(-r*dt)
    
    # initialise asset prices at maturity - Time step N
    S = np.zeros(N+1)
    S[0] = S0*d**N
    for j in range(1,N+1):
        S[j] = S[j-1]*u/d
    
    # initialise option values at maturity
    C = np.zeros(N+1)
    for j in range(0,N+1):
        if opttype == 'C':
            C[j] = max(0, S[j]-K)
        else:
            C[j] = max(0, K - S[j])
        
    # step backwards through tree
    for i in np.arange(N,0,-1):
        for j in range(0,i):
            C[j] = disc * ( q*C[j+1] + (1-q)*C[j] )
    
    return C[0]

JR_method(K,T,S0,r,N,sigma,opttype='C')

Equal Probabilities (EQP) Method

Here we choose equal risk-neutral probabilities, under logarithmic asset pricing tree

def EQP_method(K,T,S0,r,N,sigma,opttype='C'):
    #precomute constants
    dt = T/N
    nu = r - 0.5*sigma**2
    dxu = 0.5*nu*dt + 0.5*np.sqrt(4*sigma**2 * dt - 3*nu**2 * dt**2)
    dxd = 1.5*nu*dt - 0.5*np.sqrt(4*sigma**2 * dt - 3*nu**2 * dt**2)
    pu = 0.5
    pd = 1-pu
    disc = np.exp(-r*dt)
    
    # initialise asset prices at maturity - Time step N
    S = np.zeros(N+1)
    S[0] = S0*np.exp(N*dxd)
    for j in range(1,N+1):
        S[j] = S[j-1]*np.exp(dxu - dxd)
    
    # initialise option values at maturity
    C = np.zeros(N+1)
    for j in range(0,N+1):
        if opttype == 'C':
            C[j] = max(0, S[j]-K)
        else:
            C[j] = max(0, K - S[j])
        
    # step backwards through tree
    for i in np.arange(N,0,-1):
        for j in range(0,i):
            C[j] = disc * ( pu*C[j+1] + pd*C[j] )
    
    return C[0]

EQP_method(K,T,S0,r,N,sigma,opttype='C')

Trigeorgis (TRG) Method

Here we choose equal jump sizes, under logarithmic asset pricing tree

def TRG_method(K,T,S0,r,N,sigma,opttype='C'):
    #precomute constants
    dt = T/N
    nu = r - 0.5*sigma**2
    dxu = np.sqrt(sigma**2 * dt + nu**2 * dt**2)
    dxd = -dxu
    pu = 0.5 + 0.5*nu*dt/dxu
    pd = 1-pu
    disc = np.exp(-r*dt)
    
    # initialise asset prices at maturity - Time step N
    S = np.zeros(N+1)
    S[0] = S0*np.exp(N*dxd)
    for j in range(1,N+1):
        S[j] = S[j-1]*np.exp(dxu - dxd)
    
    # initialise option values at maturity
    C = np.zeros(N+1)
    for j in range(0,N+1):
        if opttype == 'C':
            C[j] = max(0, S[j]-K)
        else:
            C[j] = max(0, K - S[j])
        
    # step backwards through tree
    for i in np.arange(N,0,-1):
        for j in range(0,i):
            C[j] = disc * ( pu*C[j+1] + pd*C[j] )
    
    return C[0]

TRG_method(K,T,S0,r,N,sigma,opttype='C')

Comparision of Methods

Now we will compare convergence as function of time steps

from py_vollib.black_scholes import black_scholes as bs
import matplotlib.pyplot as plt

## call option with different steps 
CRR, JR, EQP, TRG = [],[],[],[]

periods = range(10,500,10)
for N in periods:
    CRR.append(CRR_method(K,T,S0,r,N,sigma,opttype='C'))
    JR.append(JR_method(K,T,S0,r,N,sigma,opttype='C'))
    EQP.append(EQP_method(K,T,S0,r,N,sigma,opttype='C'))
    TRG.append(TRG_method(K,T,S0,r,N,sigma,opttype='C'))
    
BS = [bs('c', S0, K, T, r, sigma) for i in periods]

plt.plot(periods, CRR, label='Cox_Ross_Rubinstein')
plt.plot(periods, JR, label='Jarrow_Rudd')
plt.plot(periods, EQP, label='EQP')
plt.plot(periods, TRG, 'r--',label='Trigeorgis')
plt.plot(periods, BS, 'k',label='Black-Scholes')
plt.legend(loc='upper right')
plt.show()