# Probability Density Function for Prices from a GBM Process in R

We’ve been working on visualisation tools to make option pricing models intuitive to the non-mathematician.

Fundamental to such an exercise is a way to model the random nature of asset price processes. The Geometric Brownian Motion (GBM) model is a ubiquitous way to do this.

We can represent the price of an asset at time t as the state x(t) of a GBM process.

x(t) satisfies an Ito differential equation dx(t) = \mu x(t) dt + \sigma x(t) dw(t) where w(t) follows a wiener process with drift \mu and volatility \sigma .

The probability distribution of future prices follows a log-normal distribution [(\mu – \frac{\sigma^2}{2}) + \log{x_0}, \sigma\sqrt{t}]

OK, nerd, but how do I get the probability distribution of future prices from the starting price of the asset, and assumptions about the return distribution?

I couldn’t work that out quickly, so I asked Wolfram Mathematica that question.

PDF[GeometricBrownianMotionProcess[\[Mu], \[Sigma], Subscript[x, 0]][t], x]\frac{exp{\frac{(-t(\mu – \frac{\sigma^2}{2}) + \log{x} – log{x_0})^2}{2t\sigma^2}}}{\sqrt{2\pi}\sqrt{t}x\sigma}, x > 0

Woah! Amazing…

Now I can code that in R (which I’m using for my modelling):

# Calculate probability density of GBM price process for price S at time t # S - price to get probability for # mu - mean of return process (for 1 step) # sd - sd of return process (for 1 step) # t - number of steps gbmpdf <- function(x, mu, sig, x0, t) { if (x == 0) return(0) exp(-(((-t*(mu - (0.5*sig^2))) + log(x) - log(x0))^2) / (2*t*sig^2)) / (sqrt(2*pi) * sqrt(t) * x * sig) }

Let’s use this function to plot the probability distribution of 30 step ahead prices for the case where:

- starting price = $100
- expected returns = 0%
- annualised volatility = 30%

# Example, stock at $100 with no drift, annualised vol 30% in 30 days time x0 <- 100 mu <- 0 sig <- 0.3 / sqrt(252) t <- 30 x <- seq(50, 150, by = 1) y <- sapply(x,gbmpdf, mu=mu, sig=sig, x0=x0, t=t) data <- data.frame(x=x, y=y) data %>% ggplot(aes(x=x, y=y)) + geom_line() + ggtitle(paste('PDF of Prices from GBM Process: x0 =', x0, ', mu =', mu, ', sig =', round(sig,2), ', t =', t))

*Note the slightly positively skewed distribution, and the peak of the distribution occurring slightly lower than the starting price…*