Simulating Geometric Brownian Motion in Python and Checking Its Moments

For constant drift

\mu

and volatility

\sigma

, sample independent standard normals

Z_i

and set

S_{i+1}=S_i\exp\!\left((\mu-\tfrac12\sigma^2)h+\sigma\sqrt h\,Z_i\right).

This is the exact transition at grid times. Compare the simulated terminal mean with

S_0e^{\mu T}

and the simulated variance with

S_0^2e^{2\mu T}(e^{\sigma^2T}-1)

Derive the transition and the theoretical moments

For

dS_t=\mu S_t\,dt+\sigma S_t\,dW_t

, apply Itô's lemma to the log to get

d\log S_t=(\mu-\tfrac12\sigma^2)\,dt+\sigma\,dW_t.

Independent Brownian increments supply the normal steps in the simulation rule. Summing those increments yields a terminal log price with mean $\log S_0+(\mu-\tfrac12\sigma^2)T$ and variance $\sigma^2T$ . Exponentiating and applying the normal moment-generating function gives

\mathbb E[S_T^k]=S_0^k \exp\!\left(k\mu T+\tfrac12k(k-1)\sigma^2T\right).

Use $k=1,2$ and subtract the squared first moment to obtain the variance. The terminal distribution is independent of how many grid times you choose. A finer grid matters for path summaries such as monitored barrier crossings.

Set an error bar before looking at the simulation

With

S_0=100

\mu=0.08

\sigma=0.25

, and

T=1

, the theoretical mean is about

108.33

and the standard deviation is about

27.51

. A Monte Carlo mean from 40,000 independent terminal prices has a standard error of roughly

27.51/\sqrt{40{,}000}\approx0.14

. A result a few cents away from

108.33

is consistent with the model. The terminal standard deviation measures price dispersion and the much smaller standard error measures the uncertainty of the estimated mean.

You can isolate two classes of error with separate tests. First simulate terminal prices directly with one normal draw per path and compare their moments with theory. Then simulate the path with many increments and compare the terminal moments again. Both are exact on the grid for constant coefficients. If the path version changes the terminal distribution appreciably as you add steps, check the

\sqrt h

factor and the sum over the time axis. For a barrier option, finer monitoring changes which paths hit the barrier even with exact price transitions, so that is a different source of error.

Vectorised simulation in Python

import numpy as np

rng = np.random.default_rng(42)
s0, mu, sigma, horizon = 100.0, 0.08, 0.25, 1.0
paths, steps = 40_000, 52
dt = horizon / steps

log_steps = (
    (mu - 0.5 * sigma**2) * dt
    + sigma * np.sqrt(dt) * rng.standard_normal((paths, steps))
)
prices = s0 * np.exp(np.cumsum(log_steps, axis=1))
terminal = prices[:, -1]

expected_mean = s0 * np.exp(mu * horizon)
expected_var = expected_mean**2 * np.expm1(sigma**2 * horizon)
simulated_sd = terminal.std(ddof=1)

print(f"Mean: {terminal.mean():.3f} vs {expected_mean:.3f}")
print(f"Variance: {simulated_sd**2:.3f} vs {expected_var:.3f}")
print(f"Mean standard error: {simulated_sd / np.sqrt(paths):.3f}")

For terminal values alone, draw one normal per path with

h=T

. Keep the full grid when a payoff or statistic needs intermediate observations.

Simulating Geometric Brownian Motion in Python and Checking Its Moments

Derive the transition and the theoretical moments

Set an error bar before looking at the simulation

Vectorised simulation in Python

Change the drift and volatility of geometric Brownian motion

References and next steps