import numpy as np
from numpy.random import default_rng
import matplotlib.pyplot as plt
Nsim = 5
t0 = 0
t1 = 1
dt = 1/250
sqrtdt=np.sqrt(dt)
Nt = int((t1-t0)/dt)
mu=.5
sigma=0.2
S0 = 1
t = np.linspace(t0,t1,Nt)
S = np.zeros([Nsim,Nt])
S[:,0] = S0

rng = default_rng()
for j in range(0, Nt-1):
    S[:,j+1] = S[:,j]*(1+sigma*sqrtdt*rng.normal(0,1, Nsim)+mu*dt)
for j in range(0,Nsim):
    plt.plot(t,S[j,:])

Nsim=1000
S = np.zeros([Nsim,Nt])
S[:,0] = S0
for j in range(0, Nt-1):
    S[:,j+1] = S[:,j]*(1+sigma*sqrtdt*rng.normal(0,1, Nsim)+mu*dt)

mean=0
for j in range(0,Nsim):
    mean+=S[j,-1]
print('mean =',mean/Nsim)
print(S0*np.exp(mu*t1))

plt.show()