Why does changing the step size in my Binomial Tree changes the final stock prices so much?

Question

I am trying to price a convertible bond by using a binomial tree. For this, I wrote a binomial tree for the stock price. I noticed that changing the step size (timesteps), changes the final value of my stock prices significantly. Intuitively this is wrong as the up and down movements should be scaled accordingly. I can imagine this will also change the price of the convertible bond.

I thought setting $u = \text{e}^{\sigma\sqrt{dt}}$ would do the trick. Could anyone tell me what I am doing wrong? Thanks!

import numpy as np
import math as math

S0 = 100 #Initial stock price
T = 5 #Maturity
timesteps = 16 #Amount of steps in the three
dt = T/timesteps #Step size
sigma = 0.28 #Vol
r = 0.01 #Interest rate
N = 300 #Notional amount (for convertible bond)
kappa = N/S0 #Conversion rate (for convertible bond)
c = 0.05 #Coupomn rate (for convertible bond)

u = np.exp(sigma*math.sqrt(dt))
d = 1/u
p = (np.exp(r*dt)-d)/(u-d)


S = np.zeros((timesteps,timesteps))
for i in range(timesteps):
    for j in range(timesteps):
        S[j,i] = S0*(u**(i-j))*(d**j)

S = np.triu(S)

Answer 1

You only got one minor bug, but let me explain why the range increases.

Let us denote $n:=timesteps$, then

1. You are looping one iteration too little when filling your $S$ matrix array, causing you to have S(T-dt) and not S(T) as terminal values. This is because you are not accounting for your starting position, i.e. you need $1+n$ iterations in each dimension.

...
S = np.zeros((1+timesteps, 1+timesteps)) # Include starting position too (S0).
for i in range(S.shape[1]):
    for j in range(S.shape[0]):
        S[j,i] = S0*(u**(i-j))*(d**j)

S = np.triu(S)

2. The range of possible values of $S(T)$ does correctly increase with $n$. The range is given by $[S_0 d^n, S_0 u^n] = [S_0 e^{-\sigma\sqrt{nT}}, S_0 e^{\sigma\sqrt{nT}}]$. Note that while the range increases, the probability for ending at these extreme values decreases. In fact, according to the Central Limit Theorem, in the limit as $n\to \infty$ the binomial distribution will become the continuous log-Normal distribution, which indeed has infinite positive support $(0,+\infty)$.
3. What the CRR binomial model ensures is that the mean and variance of the discrete binomial model matches those of the continuous model. The mean of the stock will match exactly for any $n$, whereas the variance will asymptotically approach the continuous case. The reason for the variance not being exactly matched for any $n$ is that when finding the parameter $u$ by matching the variance, approximations were made by only keeping first-order terms in a Taylor series expansion.

# ====================================================
# === Compare moments to the continuous exact ones ===
# ====================================================

from scipy import stats as stats

# Binomial model for number of down moves
pd = 1-p
dist = stats.binom(n=timesteps, p=pd)

# Stock terminal mean and variance E[S(T)/S0] and V[S(T)/S0]
# 1st moment for S(T)
bin_m1 = dist.expect(lambda k: S[k.astype(int),-1]/S0)
# 2nd moment for S(T)
bin_m2 = dist.expect(lambda k: (S[k.astype(int),-1]/S0)**2)
# Var[S_T] = E[S_T^2] - E[S_T]^2
bin_var = bin_m2 - bin_m1**2
print(f'Binomial Model: E[S/S0]={bin_m1:.4f}, V[S/S0]={bin_var:.4f}')

# Continuous S(T) moments, https://en.wikipedia.org/wiki/Geometric_Brownian_motion#Properties
gbm_m1 = np.exp(r*T)
gbm_var = np.exp(2*r*T)*(np.exp(sigma**2*T)-1)
gbm_m2 = gbm_var + gbm_m1**2
print(f'Continuous GBM: E[S/S0]={gbm_m1:.4f}, V[S/S0]={gbm_var:.4f}')

Binomial Model: E[S/S0]=1.0513, V[S/S0]=0.5218

Continuous GBM: E[S/S0]=1.0513, V[S/S0]=0.5304