Probability Density Function for Prices from a GBM Process in R

Posted on May 07, 2020 by Robot James
No Comments

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…

Leave a Comment