python：几何布朗运动模拟

Question

我在做一个简单的GBM（几何布朗运动）模拟，但好像不太对劲。我的代码总是输出小于1e-20的值，而我希望得到的是围绕1.0随机分布的值：

import math
import random

p = 1
dt = 1
mu = 0
sigma = 1
for k in range(100):
    p *= math.exp((mu - sigma * sigma / 2) * dt +
         sigma * random.normalvariate(0, dt * dt))
print(p)

我使用的环境是：

ActivePython 3.1.2.3（来自ActiveState Software Inc.），基于 Python 3.1.2（r312:79147，2010年3月22日，12:30:45）[MSC v.1500 64位（AMD64）]，在 win32上运行。

我的操作系统是Windows 7专业版，使用的是i7-930 CPU（64位）。

如果需要，我很乐意在我的机器上运行其他测试来找出问题所在。

Answer 1

看起来没问题：

import math
import random

def compute():
    p = 1
    dt = 1
    mu = 0
    sigma = 1
    for k in range(100):
        p *= math.exp((mu - sigma * sigma / 2) * dt +
                      sigma * random.normalvariate(0, dt * dt))
    return p

print [compute() for x in range(20)]

结果是：

[118.85952235362008, 7.3312246638859059e-14, 29.509674994408684, 1.8720575806397408, 1.5882398997219834e-05, 2967.524471541024, 0.0031130343677571093, 19942.669293314699, 0.00011878818261757498, 5382.80288111769, 0.22867624175360118, 0.028535167622775418, 12.6324011631628, 20.604456159054738, 0.0034504567371058613, 6.5804828930878056e-06, 6398.0493448486704, 0.0014978345496292245, 721546.38343724841, 1285.546939393781]

这是在使用 Python 2.6.1 版本

Answer 2

在Python 3.1中，使用浮点数除法（//）而不是整数除法（/）是可以正常工作的：

import math
import random

p = 1
dt = 1
mu = 0
sigma = 1
for k in range(100):
    p *= math.exp((mu - sigma * sigma // 2) * dt +
         sigma * random.normalvariate(0, dt * dt))
print(p)

在一次示例运行中，我得到了以下数字：

0.0989269233704
2536660.91466
2146.09989782
0.502233504924
0.43052439984
14.1156450335

Answer 3

我找到了答案。代码没有问题。问题在于，得到的对数正态分布的规模参数非常大，等于 1 * sqrt(100) = 10。规模参数为10时，分布的偏斜度非常大。

所以，尽管这个分布的平均值是1.0，但我可能需要进行数十亿次的迭代（如果不是数百亿次的话），才能看到一个大于1.0的数字。