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…