[R] Plot of 3d stock price density

[Maximilian Lklweryc]

Hi,
I have thought a lot about my problem and also posted it on
stackexchange, the post can be found here:
http://stats.stackexchange.com/questions/44276/how-to-plot-3d-gbm
but I have got no useful answer, that's why I trry the last
possibility, which is to ask you professional guys.

I want to recreate the following picture with my own code, the picture
can be found in my other post:
http://stackoverflow.com/questions/13387119/how-can-i-recreate-this-3d-histogram

The picture is about stock price densities, dependent on time and S.
So the underlying model is the famous geometric brownian motion, the
stock prices are therefore log-normal distributed. With increasing
time the mean increases and also the variance, for simulation I always
use the discrete version, you can see it here:
http://quant.stackexchange.com/questions/4589/how-to-simulate-stock-prices-with-a-geometric-brownian-motion

My question is now: How can I get this 3d plot? I tried first of all a
more simpler formula (x^2+y^2), where I adjusted color and things like
that:
farbe<-rgb(85, 141, 85, maxColorValue=255)

x <- seq(-40, 40, length= 10)
y <- c(1:22)
f <- function(x,y) { r <- x^2+y^2; r }
z <- outer(x, y, f)

z[is.na(z)] <- 1
op <- par(bg = "gray80")
persp(x, y, z, theta = -60, phi = 30, expand = 0.5, col = farbe,
     ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Z"
)

now I tried to put in the log-normal density formula with an
increasing mean and variance. That means the mean and variance
increase linear or? I am not really sure about the variance and how it
behaves in the GBM, I implemented it this way:
sigma<- 0.1
mu<-0.1

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

x <- seq(-40, 40, length= 10)
y <- c(1:22)
f <- function(x,y) { r <-
1/(sqrt(2*pi)*sigma*y*x)*exp(-(log(x)-mu*y)^2/2*(sigma*y)^2); r }
z <- outer(x, y, f)

z[is.na(z)] <- 1
op <- par(bg = "gray80")
persp(x, y, z, theta = -60, phi = 30, expand = 0.5, col = farbe,
     ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Z"
)

The sigma*y and mu*y should be the increasing mean and variance, I am
not sure about this. My problem is first, that the plot looks
completely wrong and the second problem is I am not sure about the
implementation of the increasing mean and variance.


It would be great if you can help me. My aim is to plot the
theoretical distribution (log-normal) and later on to plot the
empirical result of Monte Carlo simulation. So my next problem is
going to be how to put a second surface into the plot. But this is
another problem :-)

Thanks a lot for your help