这里有一些代码的重新编写，可以使S的符号更直观，并允许您检查答案的合理性。

起始点：在代码中，第二个deltat应该替换为np.sqrt(deltat)。Sourcehere(是的，我知道这不是最官方的，但下面的结果应该让人放心)。

关于短期利率和西格玛值的非年度化的评论可能不正确。这与你看到的向下漂移无关。你需要保持这些年利率。这些将始终是连续复合(恒定)速率。

首先，这里是一个GBM路径生成函数，它来自Yves Hilpisch-Python for Finance，chapter 11。链接中解释了这些参数，但设置与您的非常相似。

def gen_paths(S0, r, sigma, T, M, I):

dt = float(T) / M

paths = np.zeros((M + 1, I), np.float64)

paths[0] = S0

for t in range(1, M + 1):

rand = np.random.standard_normal(I)

paths[t] = paths[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt +

sigma * np.sqrt(dt) * rand)

return paths

设置初始值(但使用N=252，一年中的交易日数作为时间增量)：S0 = 100.

K = 100.

r = 0.05

sigma = 0.50

T = 1

N = 252

deltat = T / N

i = 1000

discount_factor = np.exp(-r * T)

然后生成路径：np.random.seed(123)

paths = gen_paths(S0, r, sigma, T, N, i)

现在，要检查：paths[-1]将在过期时获取结束的St值：np.average(paths[-1])

Out[44]: 104.47389541107971

现在的回报是最大的(St - K, 0)：CallPayoffAverage = np.average(np.maximum(0, paths[-1] - K))

CallPayoff = discount_factor * CallPayoffAverage

print(CallPayoff)

20.9973601515

如果绘制这些路径(很容易使用pd.DataFrame(paths).plot())，您将看到它们不再是向下的趋势，而是St近似对数正态分布。

最后，这里有一个通过BSM进行的健康检查：class Option(object):

"""Compute European option value, greeks, and implied volatility.

Parameters

==========

S0 : int or float

initial asset value

K : int or float

strike

T : int or float

time to expiration as a fraction of one year

r : int or float

continuously compounded risk free rate, annualized

sigma : int or float

continuously compounded standard deviation of returns

kind : str, {'call', 'put'}, default 'call'

type of option

Resources

=========

http://www.thomasho.com/mainpages/?download=&act=model&file=256

"""

def __init__(self, S0, K, T, r, sigma, kind='call'):

if kind.istitle():

kind = kind.lower()

if kind not in ['call', 'put']:

raise ValueError('Option type must be \'call\' or \'put\'')

self.kind = kind

self.S0 = S0

self.K = K

self.T = T

self.r = r

self.sigma = sigma

self.d1 = ((np.log(self.S0 / self.K)

+ (self.r + 0.5 * self.sigma ** 2) * self.T)

/ (self.sigma * np.sqrt(self.T)))

self.d2 = ((np.log(self.S0 / self.K)

+ (self.r - 0.5 * self.sigma ** 2) * self.T)

/ (self.sigma * np.sqrt(self.T)))

# Several greeks use negated terms dependent on option type

# For example, delta of call is N(d1) and delta put is N(d1) - 1

self.sub = {'call' : [0, 1, -1], 'put' : [-1, -1, 1]}

def value(self):

"""Compute option value."""

return (self.sub[self.kind][1] * self.S0

* norm.cdf(self.sub[self.kind][1] * self.d1, 0.0, 1.0)

+ self.sub[self.kind][2] * self.K * np.exp(-self.r * self.T)

* norm.cdf(self.sub[self.kind][1] * self.d2, 0.0, 1.0))

option.value()

Out[58]: 21.792604212866848

在GBM设置中为i使用更高的值应该会导致更紧密的收敛。