Simulation of Geometric Brownian Motion in R

Question

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

Answer 1

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]