Distribution of Geometric Brownian Motion

Question

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.

Thanks in advance!

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).

Answer 1

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.

Answer 2

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.