3
$\begingroup$

Recall that CRR (Cox-Ross-Rubinstein) model for option pricing is the usual binomial tree model with $u$ (up-factor) and $p$ (one of the risk-neutral probabilities) defined as follows: $$u = e^{\sigma\sqrt{\Delta t}},$$ $$p = \frac{e^{r\Delta t} - e^{-\sigma\sqrt{\Delta t}}}{e^{\sigma\sqrt{\Delta t}} - e^{-\sigma\sqrt{\Delta t}}},$$ where $\sigma$ is volatility, $r$ is interest rate, $\Delta t = \frac{T}{M}$ (time step), then $d = \frac{1}{u}$ (down-factor) and $q = 1- p$, everything else as usual.

Below is my implementation of CRR model using Python 3:

# Implementation of Cox-Ross-Rubenstein option's pricing model.

import math

T = 0.25 # time horizon
M = 2 # quantity of steps
t = T/M # step
sigma = 0.1391*math.sqrt(0.25) # volatility
r = 0.0214*0.25 # interest rate
u = math.exp(sigma*math.sqrt(t)) # up-factor
d = 1.0/u # down-factor
S0 = 2890.30 # initial underlying stock price
K = 2850 # strike

# compute risk-neutral probabilities
p = (math.exp(r*t)-math.exp(-sigma*math.sqrt(t)))/(math.exp(sigma*math.sqrt(t))-math.exp(-sigma*math.sqrt(t))) # up
q = 1 - p # down

# profit from call option
def call(stock_price, K):
    price = max(stock_price - K, 0)
    return price

# profit from put option
def put(stock_price, K):
    price = max(K - stock_price, 0)
    return price

# price for European style
def european():
    price = 1.0/(1+r)*(p*option_prices[i+1][j+1]+q*option_prices[i+1][j])
    return price

# price for American style, specify call or put in argument
def american(style):
    price = max(style, european())
    return price

stock_final_prices = []
option_final_prices = []

# create dictionary, containing lists of options prices at every time step
option_prices = {}
for i in range(0,M+1):
    option_prices[i] = [None] * (i + 1)

# calculate possible final stock prices
for i in range(0,M+1):
    stock_final_prices.append(S0*math.pow(u,i)*math.pow(d,M-i))

# calculate possible option final prices -- choose call or put function
for i in range(0,M+1):
    option_final_prices.append(put(stock_final_prices[i], K))

option_prices[M] = option_final_prices

# going backwards -- uncomment european or american function, choose call or put for american style
for i in range(M-1,-1,-1):
    for j in range(0,i+1):
        option_prices[i][j] = european()
        #option_prices[i][j] = american(call(S0*math.pow(u,j)*math.pow(d,i-j), K))

print('The price is ${0} for {1} steps.'.format(option_prices[0][0], M))

You can try to play with number of time steps variable $M$ and see that when $M$ grows option's price goes to zero, that does not make any sense. However, if you assign numerical values of $p$ and $u$ manually with $M = T$, it will become the usual binomial tree model (Black-Scholes-Merton) which works perfectly.

So, why my implementation of CRR does not converge to some meaningful non-zero price? Where I made a mistake? I'm really stuck and couldn't find it. Any help with code review will be very appreciated.

$\endgroup$
2
+50
$\begingroup$
  1. from the look of it your discounting is incorrect because as you increase M you should discount with 1/(1+r0*t) (assuming r0=0.0214 is the annual interest rate where as you seem to discount by 1/(1+r0*T))
  2. also seems there is a typo in your "r" variable, probably you want it to be r=r0=0.0214 instead of r=0.0214*0.25.

now a working rewrite of your code for the tree is

# Implementation of Cox-Ross-Rubenstein option's pricing model.

from scipy.stats import norm
import numpy as np

#OP inputs as i understand them
T = 0.25 # time horizon
M = 2 # quantity of steps
sigma = 0.1391*np.sqrt(0.25) # volatility
r0 = 0.0214
S0 = 2890.30 # initial underlying stock price
K = 2850 # strike

#size M+1 grid of stock prices simulated at time T
def stock_prices(S0,T,sigma,M):    
    res = np.zeros(M+1)    
    t = T*1.0/M # step
    u = np.exp(sigma*np.sqrt(t)) # up-factor
    d = 1.0/u
    dn = d/u
    res[0] = S0*np.power(u,M) 
    for i in range(1,M+1):
        res[i] = res[i-1] * dn
    return res

# terminal payoff from call option
def payoff(stock_price, K,kind='call'):
    epsilon = 1.0 if kind == 'call' else -1.0
    price = np.maximum(epsilon*(stock_price - K), 0)
    return price

# price for European style option using CRR
def european_crr(S0,K,T,r,sigma,M,kind='call'):
    #terminal payoff
    option_price = payoff(stock_prices(S0,T,sigma,M),K,kind)
    t = T*1.0/M # time_step
    df = np.exp(-r*t) #discount factor
    u = np.exp(sigma * np.sqrt(t))
    d = 1/u
    p = (np.exp(r*t)-d)/(u-d) # risk neutral probability for up-move        
    q=1-p    
    for time_idx in range(M): #move backward in time
        for j in range(M-time_idx):
            option_price[j] = df*(p*option_price[j]+q*option_price[j+1])
    return option_price[0]

#analytical check Black-Scholes formula (no dividend nor repo)
def european_bs(S0,K,T,r,sigma,kind='call'):
    df = np.exp(-r*T)
    F = np.exp(r*T)*S0                   # forward price with no dividend or repo
    m = np.log(F/K)/(sigma * np.sqrt(T)) #moneyness
    epsilon = 1.0 if kind == 'call' else -1.0
    d1 = m + 0.5*sigma*np.sqrt(T)
    d2 = d1 - sigma * np.sqrt(T)
    Nd1 = norm.cdf(epsilon*d1)
    Nd2 = norm.cdf(epsilon*d2)
    return epsilon*df*(F*Nd1-K*Nd2)
  • you can check that for your inputs, both analytical (european_bs) and tree (european_crr) methods converge to 17.97 for the put and 73.48 for the call using your provided inputs
  • note also that for the tree it is not necessary to hold a 2 dimensional array. It's sufficient to use a single array of size M+1 to hold the temporary results
  • when i perform both changes to your code explained in 1+2 then your results match with my results for the tree and analytical version
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.