Simulating geometric Brownian motion
3 August 2026
Geometric Brownian motion (GBM) is a frequently used modelling tool in finance and economics (especially Ergodicity Economics), but it can be tricky to understand and implement. This post compares several approaches to simulating GBM in discrete time.
GBM has the following stochastic differential equation (SDE):
This SDE describes an infinitesimal change in , during an infinitesimal slice of time, . We cannot simulate infinitesimals and so we need to discretise the equation. This means finding a way to calculate finite changes of in finite slices of time: .
There are a few different approaches we can take and comparing them helps us to understand GBM.
Each approach will use an identical source of noise, . We will take samples, and use exactly 1 sample per unit-time of the simulation. Here is a plot of these samples:
We will keep the source of noise identical across each simulation approach. This will help us compare approaches without the idiosyncratic differences arising from randomness alone.
If follows GBM, then follows Brownian motion.
We start by simulating regular Brownian motion of :
- is the drift and is the variance of the BM process of
- is a sample selected from our source of noise, which is scaled by .
We can then exponentiate to simulated GBM of :
We can plot both the trajectory of and on a symlog scale:
Just as Brownian motion is repeated addition of a normally distributed variable, GBM is repeated multiplication of a log-normally distributed variable:
This form of the equation is functionally identical to exponentiated Brownian motion, but it emphasises being the product of repeated multiplication.
We can plot the trajectory calculated by Approach 2 and compare it with Approach 1, offset slightly, and see that they are identical:
Euler-Maruyama (EM) is a workhorse approach for simulating all sorts of random processes. EM is useful when there is no closed-form expression for the process we want to simulate. However EM can only approximate GBM.
Like Approach 2, EM calculates sequential values of as the product of repeated multiplication. EM differs in that the multiplier is a normally distributed random variable that only approximates the true log-normally distributed one.
The approximation is achieved by:
- Adding to the multiplier, shifting its distribution up. This approximates the process of exponentiation around as .
- Adding the Itô correction, , to the drift term .
We can simplify this because . Substituting:
Here is our approximate Approach 3 compared to the exact Approaches 1 & 2.
This approach breaks down if is compared to or , as some of the samples from this pseudo-lognormal distribution may be . When this happens, the trajectory is multiplied by a negative number and flips into negative territory, something that is not possible in GBM.
We can plot the distribution of the lognormal multipliers (from Approach 1) and Euler-Maruyama multipliers and see how well their distributions compare.
The problems of Euler-Maruyama can be mitigated by making smaller. This concentrates the distribution of the multipliers around the mean and lessens the chance of drawing negative values from the distribution. However, making smaller requires more samples of for the same period of time.
This approach is not open to us as we have committed to use exactly 1 sample of each step of and to use the same values of for each time-step across the different simulations.
Instead we will sub-divide into steps, during which we keep constant.
We can then use the following approach to calculate sequential values of :
| Note |
|---|
| See that the Itô correction in this modified EM expression is being scaled by . |
As , then this multiplier converges on the log-normal distribution of Approaches 1 & 2.
Here is Approach 4 compared to the Approaches 1 & 2.
const mu = gamma + 0.5 * sigma ** 2const dt = 1const samples = 500let xi = (function () {
let output = [];
for (let i = 0; i < samples; i++) {
output.push(realiseNormal());
}
return output;
})();function realiseNormal() {
const u1 = Math.random();
const u2 = Math.random();
const bump = bumper;
// Box-Muller transform
const z0 = Math.sqrt(-2.0 * Math.log(u1)) * Math.cos(2.0 * Math.PI * u2);
// Scale and shift
return z0;
}const colour = d3.scaleOrdinal(d3.schemeObservable10)setLargeSigmaButton = () => htl.html`<button
onclick=${() => {
viewof sigma.value = 0.75;
viewof sigma.dispatchEvent(new Event("input", { bubbles: true }));
}}>large</button>`let approach1 = (function () {
let output = [{t: 0, x: 1}]
for (let i = 1; i <= samples; i++){
// Note: we need to log the previous value to move it back into additive space
let previousLnX = Math.log(output[i-1].x)
let nextLnX = previousLnX + gamma * dt + sigma * Math.sqrt(dt) * xi[i-1]
let nextX = Math.exp( nextLnX )
output.push({t: i * dt, x: nextX, lnX: nextLnX})
}
return output
})()let approach2 = (function () {
let output = [{t: 0, x: 1}]
for (let i = 1; i <= samples; i++){
let previousX = output[i-1].x
let nextX = previousX * Math.exp(gamma * dt + sigma * Math.sqrt(dt) * xi[i-1] )
output.push({t: i * dt, x: nextX})
}
return output
})()let approach3 = (function () {
let output = [{t: 0, x: 1}]
for (let i = 1; i <= samples; i++){
let previousX = output[i-1].x
let nextX = previousX * (1 + mu * dt + sigma * Math.sqrt(dt) * xi[i-1])
output.push({t: i * dt, x: nextX})
}
return output
})()let approach2Factors = ( function () {
let output = []
for(let i = 0; i<= samples; i++){
output.push({t: i, factor: Math.exp(gamma * dt + sigma * Math.sqrt(dt) * xi[i])})
}
return output
})()let approach3Factors = ( function () {
let output = []
for(let i = 0; i<= samples; i++){
output.push({t: i, factor: 1 + mu * dt + sigma * Math.sqrt(dt) * xi[i]})
}
return output
})()let approach4Factors = ( function () {
let output = []
for(let i = 0; i<= samples; i++){
output.push({t: i, factor: Math.pow(1+(gamma + ((sigma ** 2)/(2 * m)) * dt + sigma * Math.sqrt(dt) * xi[i]) / m, m)})
}
return output
})()