# incremental simulation of GBM

(I asked this question in stackoverflow.com, but I am now thinking my mistake may be mathematical rather than programming).

I am simulating geometric brownian motion, using closed-form solution for the value after an arbitrary (not necessarily infinitesimal) step dt. My drift = 0, standard deviation = 1, initial value = 1.

Of course, the expected value at any future point of such GBM is 1.

But, after 100 steps, I find my GBM process to be in the range between 10^-30 and 10^-15. Typically, it's around 10^-23, and falling consistently. It seems like each step my process tends to fall by a constant factor of around 1.7. If I double my standard deviation, my process falls twice faster (~3.4 times each step).

I'm sure something is wrong with my setup. Here's my Python code, which is nearly readable as pseudocode:

import math
import random
p = 1
dt = 1
mu = 0
sigma = 1
for k in range(100): # repeat loop 100 times
dW = random.normalvariate(0, math.sqrt(dt)) # draw a standard normal random value
p = p * exp((mu - sigma * sigma / 2) * dt + sigma * dW)
print(p)


Thanks...

Edit: I fixed normalvariate to use sqrt(dt) rather than dt * dt. Thanks El Moro.

## 2 Answers

Ummm I think your dW should have a standard deviation of sqrt(dt) ( variance of dt not dt*dt).

• Yes... thank you. But since dt = 1, my result is unchanged, and hence the problem remains... :( – max Apr 25 '11 at 5:47
• I think you function takes the variance as argument it should be dt instead of dt*dt or sqrt(dt) :). Did you try to see if your normalvariate func works normally? if it gives you independent normal variates ? (because if they happen to be related the algo fails) – El Moro Apr 25 '11 at 6:12
• It does take standard deviation, and in any case dt = 1. I also replaced the pseudo random number generator with a true random, with the same result. I just figured out what's wrong though (see my answer below). +1 for fixing my problem with variance. – max Apr 26 '11 at 0:42

I found the answer. There is no problem with the code. It's just that the resulting lognormal distribution has enormous scale parameter = 1 * sqrt(100) = 10. With scale of 10, the skewness is insane.

Thus, even though the mean of the distribution is 1.0, it will take me billions of iterations (if not billions of billions) to see a single number greater than 1.0.