Why my implementation of CRR model does not converge?

Question

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.

Answer 1

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