Geometric Brownian Motion Simulation

Why the mean and the median drift apart

With the default drift μ = 8% and volatility σ = 30% a year, the expected value after five years is S₀·e^0.4 ≈ 1.49·S₀. Yet the median path grows at only μ − σ²/2 = 3.5% a year, to about 1.19·S₀, and roughly 40% of paths end below where they started. Both statements are true at once, because the log-normal distribution is skewed: a minority of paths grow enormously and pull the mean up, while the typical path does much less well.

The gap is volatility drag. Gains and losses compound multiplicatively, so a 50% gain followed by a 50% loss leaves 75% of the starting value, not 100%. Over time the typical growth rate is the average log return, μ − σ²/2, not the average percentage return μ. Set σ = 0.5 with μ = 0.08: the mean still grows at 8% a year, but μ − σ²/2 = −4.5%, so the median path shrinks and most paths lose money over long horizons. The fan chart uses a log scale, on which the mean and median are straight lines with slopes μ and μ − σ²/2.

The model and its exact solution

Geometric Brownian motion is defined by the stochastic differential equation

where W is a standard Brownian motion. Itô's lemma applied to ln S removes the multiplicative noise and produces the −σ²/2 correction, which gives the exact solution

So ln(S_t/S₀) is normal with mean (μ − σ²/2)t and variance σ²t. From that: E[S_T] = S₀e^μT, median = S₀e^{(μ−σ²/2)T}, and P(S_T < S₀) = Φ(−(μ − σ²/2)√T/σ).

The simulation steps the exact solution over Δ = T/steps:

This is not an approximation: the values on the time grid have exactly the right joint distribution for any number of steps, which only controls how finely the paths are drawn. The Euler–Maruyama scheme, S(t + Δ) = S(t)(1 + μΔ + σ√Δ·Z), has a bias that shrinks with Δ and can even produce negative prices when σ√Δ is large.

Reading the lab

The histogram of final values is the log-normal density: bounded below by zero, with a long right tail. The green line (sample median) sits left of the orange line (sample mean). The sample mean is also the noisiest statistic here, because it is dominated by the few largest paths — its standard error is in the readouts, and with σ near 1 and long horizons it can miss the theoretical mean by a wide margin even with 10,000 paths. The median and the share of paths below S₀ are much more stable.

This same process, with μ replaced by the risk-free rate, is what Monte Carlo option pricing samples at maturity.

Example result

With S₀ = 100, μ = 0.08, σ = 0.3, a 5-year horizon, 250 steps, 1,000 paths and seed 42, the simulated mean final value is 151.33 (theory 149.18, standard error 3.57) and the median is 124.09 (theory 119.12). 39.0% of paths end below S₀, against a theoretical 39.7%.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Exact discretisation of dS = μS dt + σS dW: over a step Δ,
// S(t + Δ) = S(t)·exp((μ − σ²/2)Δ + σ√Δ·Z), Z ~ N(0, 1). No time-step error,
// because log S is a Brownian motion with drift and its increments are exactly normal.
export function simulatePath(p, rng, out) {
  const dt = p.years / p.steps;
  const drift = (p.mu - (p.sigma * p.sigma) / 2) * dt;
  const vol = p.sigma * Math.sqrt(dt);
  let logS = 0; // log(S/S₀)
  if (out) out[0] = p.s0;
  for (let k = 1; k <= p.steps; k++) {
    logS += drift + vol * rng.normal();
    if (out) out[k] = p.s0 * Math.exp(logS);
  }
  return p.s0 * Math.exp(logS);
}

def simulate_path(p, rng, out=None):
    """Exact step: S(t + dt) = S(t) * exp((mu - sigma^2/2) dt + sigma sqrt(dt) Z)."""
    dt = p["years"] / p["steps"]
    drift = (p["mu"] - p["sigma"] * p["sigma"] / 2) * dt
    vol = p["sigma"] * math.sqrt(dt)
    log_s = 0.0  # log(S / S0)
    if out is not None:
        out[0] = p["s0"]
    for k in range(1, p["steps"] + 1):
        log_s += drift + vol * rng.normal()
        if out is not None:
            out[k] = p["s0"] * math.exp(log_s)
    return p["s0"] * math.exp(log_s)

Run the full model from a terminal:

node simulation.js --s0 100 --mu 0.08 --seed 42
python simulation.py --s0 100 --mu 0.08 --seed 42

Download the model

Assumptions & limitations

• Constant drift and volatility. μ and σ never change. Real volatility clusters: calm and turbulent periods alternate, which models such as GARCH or stochastic volatility capture.
• Log-normal returns, no jumps. Log returns are exactly normal and paths are continuous. Real daily returns have fat tails — crashes of many standard deviations happen far more often than a normal distribution allows — and prices can gap.
• Independent increments. Past returns carry no information about future ones; there is no mean reversion or momentum.
• Frictionless and educational. No dividends, costs or taxes. The model illustrates a stochastic process; it is not a forecast of any asset.
• Pseudo-random normals. Increments use Box–Muller on a seeded generator; see pseudo-random numbers and seeds.

Related simulations

• OS-013Brownian motionParticles jiggle with Gaussian steps; their mean squared displacement grows as 4Dt, as Einstein predicted.
• OS-004Monte Carlo option pricingPrice a European option as a discounted average payoff and check it against Black–Scholes.
• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.
• OS-011Central limit theoremAverage n draws from almost any distribution and the means turn normal — at a speed you can measure.

Background: Pseudo-Random Numbers and Seeds in Simulation

References

1. Black, F. & Scholes, M. (1973). The pricing of options and corporate liabilities. Journal of Political Economy, 81(3), 637–654. doi.org/10.1086/260062
2. Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer.