Which formula for GBM is correct? I am trying to write a simple GBM simulator. Unfortunately, the task has turned rather difficult.
The first approach I looked into was the most obvious. I could use the analytic solution for the GBM found here:
$S_t = S_0exp((\mu - \frac{\sigma^2}{2})t + \sigma W_t)$
This is relatively straight-forward in implementation except the $W_t$, the $t$-th observation from a brownian process. Following the direction of this tutorial on brownian motion left me lost as to it's implementation. It seemed to be:


*

*Simulate a brownian process from $W_0$ to $W_t$

*For each $x$ from $0$ to $t$ use $W_x$ to solve for $S_x$


This seemed excessively complicated compared to what I remembered from doing simulations a few years back. Stumbling on this article I came across an interesting formula I've seen before:
$\frac{\Delta S}{S} = \mu \Delta t + \sigma \epsilon \sqrt{\Delta t}$
Where $\epsilon$ is a shock drawn from the standard normal distribution $N(0, 1)$. $\Delta  t$ is obviously the change in time from $S_t$ to $S_{t + \Delta t}$. Perhaps for daily data $\Delta t = 1$. This is a little closer to what I remember - though there's no brownian component here - $W_t$. Where did it go?
I am unsure which of these I should use that would be more correct. Currently I'm relearning stochastic calculus after many years and am inclined to use the SDE (if I could figure out how to simulate brownian motion easily!). However, the simpler second solution intrigues me as well if I could figure out how it is derived.
 A: The geometric Brownian motion is defined by the SDE:
$$dS_t = \mu S_tdt + \sigma S_tdW_t$$
The first formula you show is the exact, analytic solution to this equation. The second formula is an approximation to the SDE known as the Euler-Maruyama method in general. Basically, replace infinitesimal increments by small deltas, and evaluate the drift and diffusion functions at the start of the interval (i.e. at $t$), like this:
$$S_{t+\Delta t} - S_t= \mu S_t\Delta t + \sigma S_t \left(W_{t+\Delta t}-W_t\right)$$
By definition of the Brownian motion, the increments are distributed as:
$$W_{t+\Delta t}-W_t \sim \mathcal{N}(0,\Delta t)$$
So if we let $\varepsilon_t \sim^{\text{iid}} \mathcal{N}(0,1)$, this is the same as what you have:
$$S_{t+\Delta t} - S_t= \mu S_t\Delta t + \sigma S_t \left(\sqrt{\Delta t}\varepsilon_t\right)$$
That is, the $\sqrt{\Delta t}\varepsilon_t$ is "where the $W_t$ went".
In the case of this specific equation, the Euler-Maruyama method is not necessary because, as you mentioned, we have an exact formula. The only issue is how to simulate $W_t$. This is not particularly difficult since the increments $W_{t+\Delta t} - W_t$ are independent and normally distributed, so you can just draw independent normal variates with variance $\Delta t$, and then compute the cumulative sum to get $W_t$ for every $t$ that you need. Once you have a path for $W_t$, plug it in to the analytic solution to get a path for $S_t$.
A: The best course of action is simulate Brownian motion via
$
W_0 = 0,
$
$
W_{t+\Delta t} = W_t + \sqrt{\Delta t}\ \rm{N(0,1)} 
$
and then plug it into the very first formula that you mentioned:
$
S_t = S_0 \exp((μ − σ^2 / 2)t + \sigma W_t).
$
With $\Delta t$ small enough you will get quite realistic trajectories, like these.
