# Why my implementation of CRR model does not converge?

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, 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.

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 = 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

#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