Convergence of GBM mean after simulation?

Question

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?

Answer 1

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}}

Answer 2

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.

Answer 3

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,])