Python代码：几何布朗运动 - 有什么问题？

Question

我对Python还比较陌生，但是为了一篇大学论文，我需要应用一些模型，并且最好使用Python。我花了几天时间研究了我附加的代码，但是我无法得出正确的结果，这段代码并没有生成一个看起来像带漂移的标准布朗运动的随机过程。我的参数如mu和sigma（期望回报或漂移和波动率）在噪声过程中只会改变斜率，而不会对其它方面产生影响。这就是我的问题，所有东西看起来都像是噪声。希望我的问题已经足够具体明确，以下是我的代码：

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

- Sebastian Schrön

---

尽量让代码易读。 - Mikhail

请看这里的Python几何布朗运动模拟：http://stackoverflow.com/questions/5774433/python-geometric-brownian-motion-simulation - gliptak

Answer 1

根据维基百科的介绍，

所以看起来

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

而不是

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

---

由于T代表时间范围，我认为t应该是

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

---

根据这些Matlab示例（here 和 here），似乎

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

不是，

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

请检查数学计算，但是我可能错了。

---

因此，将所有内容放在一起：

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

产出

- unutbu

---

t 应该从 0 到 Steps 还是从 0 到 T？ - unutbu

感谢迄今为止的支持，现在似乎我的最后一个问题是W的分布，随机变量内部的跳跃太高了，无法对股票价格进行采样，尝试更改分布的不同尝试并没有真正成功。 - Sebastian Schrön

就是这样！非常感谢！！ 从步骤W(t)到W(0)的变化可以通过时间跨度内所有随机变量的总和来显示，因此在数学上也是正确的，你救了我的周末，我非常感激！ :-) - Sebastian Schrön

在最终的代码中：X = (mu-0.5sigma**2)t + sigmaW 是错误的。它应该乘以dt，而不是t。因此，正确的行是 X = (mu-0.5sigma**2)dt + sigmaW。 - Vadikus

你能提供证据证明它是错误的吗？你的说法与维基百科和我帖子中链接到的Matlab实现相矛盾。 - unutbu

Answer 2

使用高斯定律的参数化进行额外实现，但使用正常函数（而不是standard_normal），代码稍微更短。

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