# 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…

Categories R