GBM 3d plot with R

Question

I want to plot the density of the GBM in a 3d plot. So I have on one axis the stock price, on the other the time and on the z axis the density. At the end I want to produce this graph.

The formula I tried to implement can be found on Wikipedia.

Here is my approach:

mu <- 0.1
sigma <- 0.1
S0 <- 100

color <- rgb(85, 141, 85, maxColorValue=255)

x <- seq(100, 112, length=40)
y <- seq(0.25, 1.1, length=25)

f <- function(s, t) {
  dlnorm(s, meanlog=log(S0) + ((mu - 1/2 * sigma^2) * t), 
         sdlog=sigma * sqrt(t))
}

z <- outer(x, y, f)

persp(x, y, z, theta=160, phi=25, expand=0.75, col=color,
      ticktype="detailed", xlab="s", ylab="time", zlab="density"
)

But it looks clearly wrong. So where is my mistake?

Answer 1

I think there is no mistake on your part, if you set sigma <- 0.0045 and

x <- seq(100, 112, length=100) // Lower values produce jagged edges
y <- seq(0.25, 1.1, length=60)

you'll get this:

With these parameters the density has about the same peak and the maximum of the density function also has a similar direction. Alas, a number of things are wrong with this plot: the sigma parameter has been changed and the maximum of the density function seems to decrease more slowly. However, the code produced is correct, since we can assume that plnorm is implemented correctly and the sdlog parameter is obviously correct. The mean parameter is also correct, the proof of that is left as an exercise ;)

I can imagine you're not satisfied with the above argument but the plot from Wikipedia must be wrong. The volatility of a lognormal is given by $\sqrt{(e^{\sigma^2}-1) e^{2 \mu + \sigma^2}}$. For $t=1$ this evaluates to $11.08$, this is clearly much wider than the plot on Wikipedia, maybe the author forgot to include the stock price in his calculation of $\mu$. Compare with this generated by

mu <- 0.1
sigma <- 0.1

S0 <- 100

color <- rgb(85, 141, 85, maxColorValue=255)

x <- seq(80, 130, length=100)
y <- seq(0.25, 1.1, length=60)

f <- function(s, t) {
  dlnorm(s, meanlog=log(S0) + ((mu - 1/2 * sigma^2) * t), 
         sdlog=sigma * sqrt(t))
}

z <- outer(x, y, f)

persp(x, y, z, theta=180, phi=25, expand=0.75, col=color,
      ticktype="detailed", xlab="s", ylab="time", zlab="density"
)

Answer 2

It looks like the wikipedia graph was made with a much smaller value for the volatility $\sigma$. Try 0.01.