Difference between closed form binomial option value and monte carlo simulation

Question

I am trying to calculate the price of a European call option using both the the closed form expression and a monte carlo simulation. But the value's I get from both these methods are not the same:

Closed form expression:

$$q = \frac{(1+r)-d}{u-d}$$

$$C \frac{1}{(1+r)^T} * \left [\sum \limits_{i=0}^T \binom{T}{i}*q^i*(1-q)^{T-i}*max(u^i*d^{T-i}*S_0-K, 0) \right ]$$

Python implementation of closed from expression:

import math
T = 10 # Number of periods
S0 = 8 # Starting price of stock
K = 9 # Strike price of option
r = 0.2 # Risk free interest rate
u = 1.5 # Up factor
d = 0.5 # Down factor


C = 0 #Value of call

risk_free = 1 / (1 + r)**T

q = ((1 + r) - d) / (u - d)


for i in range(T+1):
    prob = math.comb(T, i)*(q**T)*(1-q)**(T-i)
    ST = max(((u**i)*(d**(T-i))*S0)-K, 0) 
    C += ST*prob
    
    
print(risk_free*C)

Output: 4.945275514422904

Python implementation of monte carlo simulation:

import random
T = 10 # Number of periods
S0 = 8 # Starting price of stock
K = 9 # Strike price of option
r = 0.2 # Risk free interest rate
u = 1.5 # Up factor
d = 0.5 # Down factor


n = 20000 # Number of runs
for j in range(n):
    S = S0
    for i in range(T):
        S *= u if random() < q else d
    value += max(S - K, 0)
value /= n * (1 + r) ** T
print("For {} runs the value is {}".format(n, value))

Output: 6.876698097695621

I don't understand what causes this difference, because the code does produce the same values when I set T=2 and S0=10, but that input does have a different p value of 0.2 while the current input has a p value of 0.25, but i don't understand what the p value means at it is not used in the formula..

Answer 1

Okay i found the problem, my implementation of binomial pricing was wrong.

This python implementation:

T = 10 # Number of periods
S0 = 8 # Starting price of stock
K = 9 # Strike price of option
r = 0.2 # Risk free interest rate
u = 1.5 # Up factor
d = 0.5 # Down factor

C = 0

q = ((1+r) - d) / (u - d)
risk_free = 1 / ((1 + r)**T)

for i in range(0, T+1):
    prob = math.comb(T, i) * (q**i) * (1-q)**(T-i)
    ST = (u**i) * (d**(T-i)) * S0
    max_value = max(ST - K, 0)
    C += max_value * prob

print(C * risk_free)

Outputs: 6.836045774062984

Which is a lot closer to the MC output