Python中的几何布朗运动模拟

Posted 2023-02-16

【中文标题】Python中的几何布朗运动模拟

【英文标题】：Geometric Brownian Motion simulation in Python

【发布时间】：2017-12-14 17:37:23

【问题描述】：

我正在尝试在 Python 中模拟几何布朗运动，以通过蒙特卡罗模拟为欧洲看涨期权定价。我对 Python 比较陌生，我收到的答案我认为是错误的，因为它远未收敛到 BS 价格，而且由于某种原因，迭代似乎呈负面趋势。任何帮助将不胜感激。

import numpy as np
from matplotlib import pyplot as plt


S0 = 100 #initial stock price
K = 100 #strike price
r = 0.05 #risk-free interest rate
sigma = 0.50 #volatility in market
T = 1 #time in years
N = 100 #number of steps within each simulation
deltat = T/N #time step
i = 1000 #number of simulations
discount_factor = np.exp(-r*T) #discount factor

S = np.zeros([i,N])
t = range(0,N,1)



for y in range(0,i-1):
    S[y,0]=S0
    for x in range(0,N-1):
        S[y,x+1] = S[y,x]*(np.exp((r-(sigma**2)/2)*deltat + sigma*deltat*np.random.normal(0,1)))
    plt.plot(t,S[y])

plt.title('Simulations %d Steps %d Sigma %.2f r %.2f S0 %.2f' % (i, N, sigma, r, S0))
plt.xlabel('Steps')
plt.ylabel('Stock Price')
plt.show()

C = np.zeros((i-1,1), dtype=np.float16)
for y in range(0,i-1):
    C[y]=np.maximum(S[y,N-1]-K,0)

CallPayoffAverage = np.average(C)
CallPayoff = discount_factor*CallPayoffAverage
print(CallPayoff)

蒙特卡洛模拟示例（股票价格模拟）

我目前使用的是 Python 3.6.1。

提前感谢您的帮助。

【参考方案1】：

这里对代码进行了一些重写，这可能会使S 的符号更直观，并允许您检查您的答案是否合理。

初始点：

在您的代码中，第二个deltat 应替换为np.sqrt(deltat)。来源here（是的，我知道这不是最官方的，但下面的结果应该让人放心）。 关于取消短期利率和 sigma 值的评论可能不正确。这与您看到的向下漂移无关。你需要保持这些年率。这些将始终是连续复合（恒定）利率。

首先，这是一个来自 Yves Hilpisch - Python for Finance 的 GBM 路径生成函数，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，1 年的交易日数，作为时间增量）：

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()，您会发现它们不再呈下降趋势，而是Sts 近似为对数正态分布。

最后，这是通过 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 使用更高的值应该会导致更紧密的收敛。

【参考方案2】：

您似乎使用了错误的公式。

有来自Wikipedia的dS_t = S_t (r dt + sigma dW_t) 和dW_t ~ Normal(0, dt) 来自Wikipedia 所以S_(t+1) = S_t + S_t (r dt + sigma Normal(0, dt))

所以我相信这条线应该是这样的：

S[y,x+1] = S[y,x]*(1 + r*deltat + sigma*np.random.normal(0,deltat))