Monte carlo pricing of European call option

I am trying to price a European call option with an closed form expression for binomial pricing, the following formula:

Which is think i implemented correctly:

T = 2 # Number of periods
S0 = 10 # 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
p = 0.2 # risk of up and down

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)

Now I'm struggling with the following, I want to generate N sample paths of the stock price in the multi-period binomial setting, under the risk-neutrality assumption. Where each path i contains stock price Si_1, Si_T. Where based on this path i i compute the payoff in period T.

And the i get the value of the option by averaging the outcome of all sample paths:

Now I'm struggling to implement this text, I have thep=0.2 so i know the probability of the stock price going up and down, but then im not sure how to calculate the new sample path and then use that to price the option for an arbitrary number of N runs, for example 200.

You can do it this way:

# monte carlo
from random import random
value = 0
numPaths = 200
for j in range(numPaths):
    S = S0
    for i in range(T):
        S *= u if random() < q else d
    value += max(S - K, 0)
value /= numPaths * (1 + r) ** T
print(value)