# Distribution of Geometric Brownian Motion

Please let me know where I have been mistaken!

Let the SDE satisfied by the GBM $S(t)$ be $$\frac{dS(t)}{S(t)} = \mu dt + \sigma dW(t).$$

Then, the underlying BM $X(t)$ will satisfy $$dX(t) = \left( \mu - \frac{1}{2} \sigma^2 \right)dt + \sigma dW(t).$$

To simulate the GBM for times $t_0 < t_1 < \ldots < t_n$, generate $n$ iid $\mathcal{N}(0,1)$ RVs, $Z_i$, $\quad i = 1,2,\ldots, n$ and set

$$S(t_{i+1}) = S(t_i)\exp\left(\left( \mu - \frac{1}{2} \sigma^2 \right)(t_{i+1} -t_i) + \sigma \sqrt{t_{i+1} - t_i} Z_{i+1} \right).$$ Then $$\frac{S(t)}{S(0)} = \exp\left(\underbrace{\left(\mu - \frac{1}{2} \sigma^2\right) t + \sigma \sqrt{t}Z}_{\mathcal{N}\left(\left(\mu - \frac{\sigma^2}{2}\right)t, \sigma^2 t\right)} \right)$$ and so $$\log \frac{S(t)}{S(0)} \sim \mathcal{N}\left(\left(\mu - \frac{\sigma^2}{2}\right)t, \sigma^2 t\right).$$

I think I'm missing something because when I use this mean and variance (in the equation directly above) for testing the normality of the log returns ($\log \frac{S(t)}{S(0)}$), I get ridiculous answers.

To be specific, with $S_0 = 20$, $\mu = 2$, $\sigma^2 = 1$ and partitioning $[0,1]$ into 100 subintervals, generating the GBM at these 100 points gives a range of values from 15.399 to 97.1384 for the $S(t_i)$. Then the log returns, $\log S(t_i)/S_0$, range from -0.26143 to 1.5804. The means I'm using for each $\log S(t_i)/S_0$ are (in increasing order of $i$) $0.015, 0.03, 0.045, \ldots, 1.5$ and the variances are $0.01,0.02, \ldots, 1$ since these are the parameters in the supposed normal distribution as described above. Finally, when these log returns are normalized using this mean and variance their values range from -2.7643 to 0.080404, which is clearly not $\mathcal{N}(0,1)$-distributed.

UPDATE: The test I am using to test for normality (Anderson-Darling) relies on independent samples from a (supposed) normal distribution, and as a couple people have pointed out in the comments, $\log S(t_i)/S_0$ is dependent on $\log S(t_{i-1})/S_0$. Indeed, changing to testing returns of the form $\log S(t_i)/S(t_{i-1})$ results in a "pass" for the Anderon-Darling test (giving $A^2 = 0.244$ for those familiar).

• could you perhaps elaborate on how you are testing the normality :) I can't really see a problem with the derivation so perhaps the test itself is a bit off – Probilitator Mar 9 '14 at 15:54
• @ bcf : Hi in what context are you testing normality ? Is it on time series of quoted index or stock, or is it in a Monte Carlo simulation that you have done yourself ? Best regards – TheBridge Mar 9 '14 at 17:20
• Anderson-Darling test. I'm partitioning $[0,1]$ into 100 intervals and generating the GBM at those points, then performing the Anderson-Darling test for these 100 random variables. – bcf Mar 9 '14 at 17:58
• Just added specifics to my normality test in the problem description. – bcf Mar 9 '14 at 18:09
• I'm not familiar with the specifics of the Anderson-Darling test but it seems that it is based on the empirical cdf of an iid sample. I'm wondering if you took into account that by considering the $\log(S_i/S_0)$ instead of the $\log(S_i/S_{i-1})$, your sample is not iid. Even if you normalise your data so that each follows an $N(0,1)$, they will never be independant since the log-return during the period $[0,i/100]$ clearly depends on the log-return during $[0,(i-1)/100]$. – AFK Mar 10 '14 at 1:59

Overall you are not mistaken, although it is worth revisiting a few steps in your question.

We assume $S$ follows the SDE $$\dfrac{dS}{S} = \mu\:dt+ \sigma\:dW^\mathbb{P}(t)$$ under the physical measure $\mathbb{P}$. If we change to the risk neutral measure $\mathbb{Q}$ (using Girsanov's theorem) then $\mu \to r$ and we have the following SDE $$\dfrac{dS}{S} = r\:dt+ \sigma\:dW^\mathbb{Q}(t).$$ Under the Black-Scholes model we would assume that $r$, $\mu$, and $\sigma$ were constants (with $\sigma > 0$) and we could integrate the SDE in the range $[t,T]$ find that $$S(T) = S(t)\exp\left(\left(\mu - \dfrac{\sigma^2}{2}\right)(T-t) + \sigma \left(W^\mathbb{Q}(T) - W^\mathbb{Q}(t)\right)\right).$$ Now to simulate this process we only need to sample $\left(W^\mathbb{Q}(T) - W^\mathbb{Q}(t)\right)$ which is achieved by computing $\sqrt{T-t\:}Z$ where $Z\sim N(0,1)$. Following this approach you find that $\log \left(\dfrac{S(T)}{S(t)}\right)$ has the desired normal distribution.

The alternative approach you seem to have implemented is to generate the entire path for $S$ which is unnecessary in this case. However, it would be necessary if we needed to know the path at intermediate times (e.g. for an exotic derivative such as an Asian option). At this point there are many ways to simulate the path and the simplest (which you have implemented) is the Euler-Maruyama method (alternatives could include the Milstein method), where $$S(t + \Delta t) = S(t) + S(t)\times\left(\mu \Delta t + \sigma \sqrt{\Delta t\:}Z\right)$$ and would involve implementing this iteratively. This is much more expensive to compute and so should be avoided if possible, and gives $O(\Delta t)$ convergence in $S(T)$.

I give the following 2-line implementation in Python2.7.10

import numpy as np
import scipy.stats as stats
s_0, mu, sigma, T, paths, path_steps = [20.0, 2.0, 1.0, 1.0,  100, 100]
s_T = s_0 * np.exp((mu-0.5*(sigma**2))*T + sigma*np.sqrt(T)*np.random.normal(0.0, 1.0, 1000))
print("\"True\" Mean: \t" + str((mu-0.5*(sigma**2))*T))
print("Computed Mean: \t" + str(np.log(s_T/s_0).mean()))
print("\"True\" Standard Deviation: \t\t" + str(sigma*np.sqrt(T)))
print("Computed Standard Deviation: \t" + str(np.log(s_T/s_0).std()))
print("Shaprio-Wilk Normality test, p-value: \t" + str(stats.shapiro(np.log(s_T/s_0))[0]))
print("Anderson-Darling Test:")
print("\t" + str(stats.anderson(np.log(s_T/s_0), dist="norm")))


which has the following output:

"True" Mean:    1.5
Computed Mean:  1.48206641711
"True" Standard Deviation:      1.0
Computed Standard Deviation:    0.983089657019
Shaprio-Wilk Normality test, p-value:   0.998949408531
Anderson-Darling Test:
AndersonResult(statistic=0.26799531674055288, critical_values=array([ 0.574,  0.653,  0.784,  0.914,  1.088]), significance_level=array([ 15. ,  10. ,   5. ,   2.5,   1. ]))


The high p-value of the Shapiro-Wilk test supports the Null hypothesis that the result is normally distributed. I am not familiar with the Anderson-Darling test, so see the package description and interpret the result for yourself.

I hope this helps.

What you are missing is the sum of errors. Though the sum is i.i.d but if you want to match your numbers, then your need to sum them correctly. sqrt(ti+1 - ti)zi + sqrt(ti - ti-1)zi-1 is not equal to sqrt(ti+1 - ti-1)(zi-zi-1). Still getting familiar with the math symbols, but I think you get the idea. Log(St/S0) is definitely normal but to compare numerically, you need to take into account the correct sum of errors.

• $$\sqrt{t_{i+1}-t_i}z_i+\sqrt{t_{i}-t_{i-1}}z_{i-1}\neq \sqrt{t_{i+1}-t_{i-1}}\left(z_i - z_{i-1}\right)$$ is that what you meant? If so hover over the comment and see mathjax markdown :). – Chinny84 Jan 14 '15 at 21:41
• I think the comments have given the specified solution already... – SmallChess Mar 16 '15 at 4:26