9

I'm pretty new to Python, but for a paper in University I need to apply some models, using preferably Python. I spent a couple of days with the code I attached, but I can't really help, what's wrong, it's not creating a random process which looks like standard brownian motions with drift. My parameters like mu and sigma (expected return or drift and volatility) tend to change nothing but the slope of the noise process. That's my problem, it all looks like noise. Hope my problem is specific enough, here is my coode:

import math
from matplotlib.pyplot import *
from numpy import *
from numpy.random import standard_normal

'''
geometric brownian motion with drift!

Spezifikationen:

    mu=drift factor [Annahme von Risikoneutralitaet]
    sigma: volatility in %
    T: time span
    dt: lenght of steps
    S0: Stock Price in t=0
    W: Brownian Motion with Drift N[0,1] 
'''

T=1
mu=0.025
sigma=0.1
S0=20
dt=0.01

Steps=round(T/dt)

t=(arange(0, Steps))
x=arange(0, Steps)
W=(standard_normal(size=Steps)+mu*t)### standard brownian motion###
X=(mu-0.5*sigma**2)*dt+(sigma*sqrt(dt)*W) ###geometric brownian motion####
y=S0*math.e**(X)

plot(t,y)

show()
Ray
  • 2,382
  • 16
  • 21
Sebastian Schrön
  • 103
  • 1
  • 1
  • 3

2 Answers2

22

According to Wikipedia,

enter image description here

So it appears that

X=(mu-0.5*sigma**2)*t+(sigma*W) ###geometric brownian motion####

rather than

X=(mu-0.5*sigma**2)*dt+(sigma*sqrt(dt)*W)

Since T represents the time horizon, I think t should be

t = np.linspace(0, T, N)

Now, according to these Matlab examples (here and here), it appears

W = np.random.standard_normal(size = N) 
W = np.cumsum(W)*np.sqrt(dt) ### standard brownian motion ###

not,

W=(standard_normal(size=Steps)+mu*t)

Please check the math, however, I could be wrong.

So, putting it all together:

import matplotlib.pyplot as plt
import numpy as np

T = 2
mu = 0.1
sigma = 0.01
S0 = 20
dt = 0.01
N = round(T/dt)
t = np.linspace(0, T, N)
W = np.random.standard_normal(size = N) 
W = np.cumsum(W)*np.sqrt(dt) ### standard brownian motion ###
X = (mu-0.5*sigma**2)*t + sigma*W 
S = S0*np.exp(X) ### geometric brownian motion ###
plt.plot(t, S)
plt.show()

yields

enter image description here

unutbu
  • 711,858
  • 148
  • 1,594
  • 1,547
  • 1
    Man, you research much faster than I do :) +1 – RocketDonkey Nov 02 '12 at 21:20
  • Oh well, according to the literature it was my formula, but it looks a lot better this way, thanks a lot! Is there a way to increase the steps to make it look more like a continuous motion? If I just increase the step the exp-funktion will create enormous values for the stock price – Sebastian Schrön Nov 02 '12 at 21:31
  • Additionally, the parameters like sigma don't show off the desired changes to the graph :/ – Sebastian Schrön Nov 02 '12 at 21:43
  • Should `t` range from 0 to `Steps`? Or should `t` range from 0 to `T`? – unutbu Nov 02 '12 at 22:05
  • Hm, yea thanks I guess that's wrong and should be t=linspace(0, T, Steps) – Sebastian Schrön Nov 02 '12 at 22:54
  • 1
    Thanks for the support so far, now it seems that my last problem is the distribution of W, the jumps within the random variable are way to high to sample a stock price movement, different attempts to change the distribution didn't really succeed though – Sebastian Schrön Nov 02 '12 at 23:45
  • 1
    That's it! Thanks a lot!! The changes from step W(t) to W(0) can be displayed by the sum of all random variables over the time horizon, so it's mathematically correct too, you saved my weekend, I'm really grateful! :-) – Sebastian Schrön Nov 03 '12 at 10:00
  • One niggle. It looks like you added noise to your initial condition such that `S(0)` is not `S0`. An important property of the Wiener process is that `W(0) = 0`. – horchler Jul 30 '13 at 15:38
  • In the final code: X = (mu-0.5*sigma**2)*t + sigma*W is wrong. It should be multiplied by dt, not t. So, the correct line is X = (mu-0.5*sigma**2)*dt + sigma*W – Vadikus Nov 13 '13 at 20:41
  • 1
    Can you provide proof that it is wrong? Your claim contradicts Wikipedia and the Matlab implementation linked to in my post. – unutbu Nov 13 '13 at 20:51
  • @unutbu how can multivariate GBM be done, combined with Cholesky decomposition to target the correlation structure between the individual processes, in python? – develarist Sep 29 '20 at 21:31
0

An additional implementation using the parametrization of the gaussian law though the normal fonction (instead of standard_normal), a bit shorter.

import numpy as np

T = 2
mu = 0.1
sigma = 0.01
S0 = 20
dt = 0.01
N = round(T/dt)
# reversely you can specify N and then compute dt, which is more common in financial litterature

X = np.random.normal(mu * dt, sigma* np.sqrt(dt), N)
X = np.cumsum(X)
S = S0 * np.exp(X)
Thabris
  • 24
  • 5
  • It is completely unclear what you are asking. Please improve your question so that there is no investigative effort necessary to find out what you want to know. – Vroomfondel Nov 01 '17 at 12:44
  • not asking anything, just proposing another way of doingit – Thabris Nov 01 '17 at 14:51
  • Sorry, SO presented me this for review and I was under the impression that this was a question. I didn't understand that I was reviewing an answer * blushing face * - sorry again, will look twice next time. – Vroomfondel Nov 01 '17 at 18:07