# Convergence of GBM mean after simulation?

As a follow up of my previous question, I am now simulating the GBM step by step for $n$ steps.

I am using the following implementation for the simulation:

$$S_{t+1} = S_t \exp \left[ \left(\mu-\frac{\sigma^2}{2}\right) \Delta t + \sigma \sqrt{\Delta t} Z_t \right], ~ Z_t \sim \mathcal{N}(0,1)$$

Each step represents a unit of time so $\Delta t = 1$.

I use the following MATLAB code:

function paths = gbm_exp(mu,vol,s0,nbr_steps,nbr_paths)

shocks = randn(nbr_steps,nbr_paths/2);
shocks_ant = [shocks, -shocks];

paths = zeros(nbr_steps+1,nbr_paths);

paths(1,:) = s0;

for i=1:nbr_paths
for j=1:nbr_steps
paths( j + 1, i ) = paths( j, i ) * exp( (mu - vol^2/2) + vol * shocks_ant( j, i ) );
end
end


As you can see, I use antithetic path to try to reduce the overall variance.

The thing is, again, I should have $\mathbb{E}(S_t) = S_0 ~ \forall t$ if I set $\mu=0$.

So I do the following:

>> test=gbm_exp(0,.3,100,300,2000);
>> mean(test(end,:))


This means that I simulate 300 steps, with $\sigma=0.3$, $\mu=0$ and $S_0 = 100$.

The mean I get though, is something quite small, around 20 on average. So, not at all the excpected 100.

So, I tried increasing my number of paths to 200k and I get a mean of roughly 40 on average.

So I'm suprised by this behavior, I would expect it to converge much more quickly, especially with antithetic paths.

Did I miss something obvious again?

## 3 Answers

You should look at confidence interval. Normally, your confidence interval size is proportional to the standard deviation, looking something like: with probability $p$ your value will be in the interval: $$[\bar{S} - k*StdDev, \bar{S} + k*StdDev]$$

Then, getting back to your simulation, we can say that your time step is very big (1 year) and you simulate 300 points, meaning that you simulate outcome in 300 years.

$$\mathbb{E}(S_t) = S_0e^{\mu t} \\ \mathbb{Var}(S_t) = S_0^2e^{2\mu t} (e^{\sigma^2t}-1)$$

In your case, e.g. $\{\mu = 0\}$ we obtain: $$\mathbb{E}(S_t) = S_0 \\ \mathbb{Var}(S_t) = S_0^2(e^{\sigma^2t}-1)$$

Plugging in $\{S_0 = 100, \sigma =.3, t = 300\}$ we obtain: $$\mathbb{E}(S_t) = 100 \\ \mathbb{Var}(S_t) = 100^2(e^{27}-1)$$

I tested similar code with Mathematica (without antithetic) for 300 time steps with 1 year step and got {Paths, Mean, StdDev, {Min, Max}} = {10^6, 60.2203, 16363.3, {2.40216*10^-15, 1.2311*10^7}}. As you can see the standard deviation is very important, it means that your confidence interval will be also big.

Just try to test for 1 year and you should get good average: {Paths, Mean, StdDev, {Min, Max}} = {10^6, 100.028, 30.6978, {22.5584, 370.411}}

• Ok, this is a good answer, but it's still incomplete in my opinion. The variance and mean of the process are key. According to the central limit theorem, the mean of all $S_{300}$ is normally distributed with mean $\mathbb{E}[S_{300}]$ and variance $\frac{Var[S_{300}]}{n}$. These values can be computed as you illustrate. – SRKX Sep 5 '14 at 6:07
• I would, however, remove you last paragraph, as I think it's missing the point. Of course it works for low-volatility or short periods of times, I'm in particular trying to test what happens in more extreme cases. – SRKX Sep 5 '14 at 6:08

The GBM is a continuous model, so using large integer time steps naturally leaves large discretization error (which vanishes when you increase the number of steps).

Use small time step 0.001:

paths(j + 1,i) = paths(j,i) * exp((mu - vol^2/2)*0.001 + vol * 0.001^0.5*shocks_ant(j,i));


Then the mean is almost exactly 100 as expected.

appearantly your sampling variance is too large. I reimplemented your example in R. What I first saw is, that the mean got worse if I took more time steps (you take $300$). Your volatility is $0.3$ which is $30\%$ per year and you sample $300$ years. What you should do is the following:

• define a variable  nbr_steps_peryear
• choose the number of  years
• then you take  nbr_steps = years*nbr_steps_peryear
• define a variable  dt = 1/nbr_steps_peryear

The dt will control your variance. I get nice results with the following set-up:

nbr_steps_peryear = 20
years = 5
nbr_steps = years*nbr_steps_peryear
dt = 1/nbr_steps_peryear
nbr_paths = 10^3
S0 = 100
mu = 0.0
vol = 0.3


The dt term balances between the drift term $\mu-\sigma^2/2$ and the volatility term.

The main source of the problem is the large number of years - as already mentioned above. With a reasonable number of years and intermediate time steps to balance th drift it should look fine. nbr_steps_peryear = 20 years = 5 nbr_steps = years*nbr_steps_peryear dt = 1/nbr_steps_peryear nbr_paths = 10^3 S0 = 100 mu = 0.0 vol = 0.3

shocks <- replicate(nbr_paths/2,rnorm(nbr_steps))
shocks_ant = cbind(shocks,-shocks)

paths = matrix(0,nrow=nbr_steps+1,ncol=nbr_paths)
paths[1,]  = S0

for (j in 1:nbr_steps){
paths[j+1,] = paths[j,]*exp((mu-vol^2/2)*dt + vol*sqrt(dt)*shocks_ant[j,])
}

mean(paths[nbr_steps+1,])
summary(paths[nbr_steps+1,])
hist(paths[nbr_steps+1,])

• Well, I agree that high volatility is the source, but I don't think it's a problem. In the end, the model should work alright even with high volatilities and long samples. I was trying to understand whether results looked rational for that kind of $\sigma$... – SRKX Sep 4 '14 at 10:56
• It is not the value $0.3$ ... but the discretization error if you sample full years. You have to discretize the year. You have $300$ one-year steps. Try to cut a year in, sy $20$ steps and simulate first $1,2,3$ years and look what happens. You have $\Delta t$ in you formula but you don't use it in the code. – Richard Sep 4 '14 at 11:39
• I understand what you say, but I think the problem would be exactly the same if I take $n = 300 \cdot \frac{1}{\Delta t}$ steps... (i assume $\Delta t$ is expressed as a fraction of a year here...) – SRKX Sep 5 '14 at 5:59
• Besides, I'm not using the euler approximation here. So the number of steps I take won't change anything; I have the exact solution anyway... – SRKX Sep 5 '14 at 6:10
• If you choose $\Delta t$ small, then you sample more often during one year - right. I will reread some works about MC simulation of stochastic processes/SDEs but I think smaller $\Delta t$ is better. if not - then why not take just a single time step of 300 years? – Richard Sep 5 '14 at 7:33