# Simulation of Geometric Brownian Motion in R

Using R, I would like to simulate a sample path of a geometric Brownian motion using

$$\begin{equation*} S(t) = S(0) \exp\left(\left(\mu - \frac{\sigma^{2}}{2}\right)t + \sigma B_{t}\right), \end{equation*}$$

where $$(B_t)$$ is the Wiener process, i.e. $$B_t\sim N(0,t)$$ for all $$t$$.

I would like to compare this path with the one that I get using the Euler- Maruyama scheme:

$$\begin{equation*} S(i+1) = S(i) + mu*S(i)*delta_t + sigma*S(i)*B_{t} \end{equation*}$$

I would like to reproduce the graph at page 534 in the paper Higham (2001)"An algorithmic introduction to numerical simulation of SDE":

I got a wrong result using the code:

rm(list=ls())
#Simulating  Geometric Brownian motion (GMB)
tau <- 1 #time to expiry
N <- 1000 #number of sub intervals
dt <- tau/N #length of each time sub interval
time <- seq(from=0, to=tau, by=dt) #time moments in which we simulate the process
length(time) #it should be N+1

mu <- 0.05 #GBM parameter 1
sigma <- 0.9 #GBM parameter 2
X0 <- 10 #initial condition

#simulate 1 Geometric Brownian motion path
Z <- rnorm(N, mean = 0, sd = 1) #standard normal sample of N elements
dW <- Z*sqrt(dt) #Brownian motion increments
W <- c(0, cumsum(dW)) #Brownian motion at each time instant N+1 elements

#Analytic solution
X_analytic <- numeric(N+1) #vector of zeros, N+1 elements
X_analytic[1] <- X0 #first element of X_analytic is X0. with the for loop we find the other N elements

for(i in 2:length(X_analytic)){
X_analytic[i] <- X_analytic[1]*exp(mu - 0.5*sigma^2*i*dt + sigma*W[i-1])
}

#plot X against time
plot(time, X_analytic, type = "l", main = "GBM path with analytical solution",
xlab = expression("t"[i]), ylab = expression("W"[t[i]]))

#Euler-Maruyama scheme
X_EM <- numeric(N+1) #vector of zeros, N+1 elements
X_EM[1] <- X0 #first element of X_EM is X0. with the for loop we find the other N elements

for(i in 2:length(X_EM)){
X_EM[i] <- X_EM[i-1] + mu*X_EM[i-1]*dt + sigma*dW[i-1]
}

#plot X against time
plot(time, X_EM, type = "l", main = "GBM path with Euler-Maruyama scheme",
xlab = expression("t"[i]), ylab = expression("W"[t[i]]))

#plot W against time
matplot(time, cbind(X_analytic, X_EM), type = "l", main = "GBM",
xlab = expression("t"[i]), ylab = expression("X"[t[i]]))


This is the result:

I don’t know which one is wrong and why

The issue is that you do not plot one sample path but for each time point $$t$$, you simply plot one possible realisation of the random variable $$S_t(\omega)$$. Thus, you don't get a connected path.

(Just as a minor, you would need brackets in the exponential in your for loop, i.e.

X_analytic[i] <- X_analytic[1]*exp((mu - 0.5*sigma^2)*time[i] + sigma*Z[i-1]*sqrt(time[i]))


Nonetheless, in order to simulate a sample path of a geometric Brownian motion, note that \begin{align*} S_{t_i}=S_{t_{i-1}}\cdot\exp\left(\left(\mu-\frac{1}{2}\sigma^2\right)(t_i-t_{i-1})+\sigma B_{t_i-t_{i-1}}\right) \end{align*} In your case, you chose a fixed step size $$\Delta t=t_i-t_{i-1}$$ for all $$i$$ such that \begin{align*} S_{t_i}=S_{t_{i-1}}\cdot\exp\left(\left(\mu-\frac{1}{2}\sigma^2\right)\Delta t+\sigma \sqrt{\Delta t}Z\right), \end{align*} where $$Z\sim N(0,1)$$.

Thus, change your for loop line to

X_analytic[i] <- X_analytic[i-1]*exp((mu - 0.5*sigma^2)*dt + sigma*Z[i-1]*sqrt(dt))


Furthermore, you may want to change the line time <- seq(from=0, to=1, by=dt) to time <- seq(from=0, to=tau, by=dt) such that you can actually make use of tau. Finally, the line sigma <- 0.9 is a bit ambitious, a volatility of 90% is rather high. If you're modelling stock prices, a value of 0.1 to 0.4 is more appropriate.

## Edit in response to your updated code

You do not need the first for loop to compute the analytical'' solution. Simply use

X_analytic = X0 * exp((mu-0.5*sigma^2)*time+sigma*W)


Secondly, in your Euler approximation, you missed $$S$$ in the last term, so the for loop step should read as

X_EM[i] <- X_EM[i-1] + mu*X_EM[i-1]*dt + sigma*X_EM[i-1]*dW[i-1]

• Why would you use the euler discretization when you know the analytical distribution? – Andrew Dec 1 '19 at 12:46
• @Andrew as I said in the answer, the approach above which is indeed a version of the Euler Maruyama algorithm, ensures that you can plot the sample path afterwards and it indeed looks like a geometric Brownian motion. The initial proposal leads to completely disconnected realisations of a geometric Brownian motion. – Kevin Dec 1 '19 at 13:06
• how would the code change for simulating multivariate correlated Brownian motion time series using Cholesky method, where some of the assets can be set correlated to one another – develarist Dec 1 '19 at 13:30
• @develarist Sorry but I fail to see how this is related the OPs question in any way? There is only one single Brownian motion driving the process? One applies the Cholesky decomposition to the covariance matrix to generate sample paths of several correlated processes but this has nothing to do with the original question? – Kevin Dec 1 '19 at 13:37
• @KeSchn sure thanks, probably I'll get back to you for pricing barrier options with Monte Carlo during the week! I'm teaching a course at University about computational finance for the first time so I'm writing all the codes :) – luca dibo Dec 1 '19 at 17:09