GBM 3d plot with R

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?

-
I'm starting to doubt the graph given on Wikipedia. The peakedness of the density seems to correspond to a volatility of 0.005. –  Bob Jansen Dec 2 '12 at 21:14
the density function can go to 1.5?? –  SRKX Dec 3 '12 at 21:36
@SRKX the function can certainly go to 1.5 or arbitrarily high with a narrow distribution -- the integral over the $S$ dimension simply has to be 1.0. Think of Dirac delta functions. –  Brian B Dec 4 '12 at 4:06
Or, less exotic, uniform over a very small interval. –  Bob Jansen Dec 5 '12 at 6:01

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"
)

-
big thanks, but your the code you posted (the long, which follows "generated by") is for what? It just calculates the density of a lognormal in a different way right? So if I inserv mu and sigma equal to 0.1 it is the same right? Thanks a lot again! –  user1690846 Dec 5 '12 at 11:03
Yes, I couldn't find any problems with your implementation, and apparently it is correct. Since dlnorm is just the log normal density in R, with the right $\mu$ and $\sigma$ they are equivalent. –  Bob Jansen Dec 5 '12 at 11:18
It looks like the wikipedia graph was made with a much smaller value for the volatility $\sigma$. Try 0.01.