蒙特卡洛估值几种不同的计算方式(Python)

１．Scipy
２．纯Python
３．Numpy(1)
４．Numpy(2)

方法一：Scipy

公式１：
$C\left({S}_{t},K,t,T,r,\sigma \right)={S}_{t}\ast N\left({d}_{1}\right)-{e}^{-r\left(T-t\right)}\ast K\ast N\left({d}_{2}\right)$$C(S_t,K,t,T,r,\sigma) = S_t * N(d_1) - e^{-r(T-t)}*K*N(d_2)$

$N\left(d\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^{d}{e}^{-\frac{1}{2}{x}^{2}}dx$$N(d) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{d}e^{-\frac{1}{2}x^2}dx$

${d}_{1}=\frac{log\left(\frac{{S}_{t}}{K}\right)+\left(r+\frac{{\sigma }^{2}}{2}\right)\left(T-t\right)}{\sigma \sqrt{T-t}}$$d_1 = \frac{log(\frac{S_t}{K})+(r+\frac{\sigma^2}{2})(T-t)}{\sigma\sqrt{T-t}}$

${d}_{2}=\frac{log\left(\frac{{S}_{t}}{K}\right)+\left(r－\frac{{\sigma }^{2}}{2}\right)\left(T-t\right)}{\sigma \sqrt{T-t}}$$d_2 = \frac{log(\frac{S_t}{K})+(r－\frac{\sigma^2}{2})(T-t)}{\sigma\sqrt{T-t}}$

${S}_{t}$$S_t$, t时刻股票或者指数的水平
$K$$K$, 行权价格
$T$$T$, 期权的到期年限(距离到期日时间间隔)
$r$$r$, 无风险利率
$\sigma$$\sigma$,波动率(收益标准差)

# 倒入用到的库
from math import log, sqrt, exp
from scipy import stats
# 期权的定价计算,根据公式１．
def bsm_call_value(S_0, K, T, r, sigma):
S_0 = float(S_0)
d_1 = (log(S_0 / K) + (r + 0.5 *sigma **2) *T)/(sigma * sqrt(T))
d_2 = (log(S_0 / K) + (r - 0.5 *sigma **2) *T)/(sigma * sqrt(T))
C_0 = (S_0 * stats.norm.cdf(d_1, 0.0, 1.0) - K * exp(-r * T) * stats.norm.cdf(d_2, 0.0, 1.0))
return C_0

# 计算的一些初始值
S_0 = 100.0    # 股票或指数初始的价格;
K = 105        #  行权价格
T = 1.0        #  期权的到期年限(距离到期日时间间隔)
r = 0.05       #   无风险利率
sigma = 0.2    # 波动率(收益标准差)
# 到期期权价值
%time print bsm_call_value(S_0, K, T, r, sigma)

8.02135223514
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 1.01 ms


8.0213522作为蒙特卡洛估值的基准值,计算耗时1.01ms,等于0.00101s

方法二：Python

$d{S}_{t}=r{S}_{t}dt+\sigma {S}_{t}d{Z}_{t}$$dS_t = rS_tdt+\sigma S_tdZ_t$

${Z}_{t}$$Z_t$是一个服从布朗运动的随机变量

${S}_{t}={S}_{t-\mathrm{\Delta }t}exp⟮⟮r-\frac{1}{2}{\sigma }^{2}⟯\mathrm{\Delta }t+\sigma \sqrt{\mathrm{\Delta }t}{z}_{t}⟯$$S_t = S_{t-\Delta t} exp\lgroup\lgroup r- \frac{1}{2}\sigma^2\rgroup \Delta t + \sigma\sqrt{\Delta t}z_t \rgroup$

１．将到期日之前的时间间隔[0,T],分割为多个等距的$\mathrm{\Delta }t$$\Delta t$子间隔
２．进行 Ｉ 次的模拟，$i\in$$i \in$ {１，２，…Ｉ}
每次模拟，循环M个子间隔, $t\in$$t \in$ {$\mathrm{\Delta }t,2\mathrm{\Delta }t,3\mathrm{\Delta }t$$\Delta t,2\Delta t,3 \Delta t$…T}
(1),每个间隔点，取一个随机数 ${z}_{t}\left(i\right)$$z_t(i)$
(2),根据离散化公式３，计算 ${S}_{T}$$S_T$
(3),计算T时刻，期权的内在价值 ${h}_{T}$$h_T$:

${h}_{T}\left({S}_{T}\left(i\right)\right)=max\left({S}_{T}\left(i\right)-K,0\right)$$h_T(S_T(i)) = max(S_T(i) - K, 0)$

３．根据Ｉ次的模拟，计算期权到期价值：

${C}_{0}={e}^{-rT}\frac{1}{T}\sum _{i=0}^{Ｉ}{h}_{T}\left({S}_{T}\left(i\right)\right)$$C_0 = e^{-rT}\frac{1}{T}\sum_{i=0}^{Ｉ}h_T(S_T(i))$

Python实现

# 纯Python实现
from time import time
from math import exp, sqrt, log
from random import gauss, seed
seed(2000)
# 计算的一些初始值
S_0 = 100.0    # 股票或指数初始的价格;
K = 105        #  行权价格
T = 1.0        #  期权的到期年限(距离到期日时间间隔)
r = 0.05       #   无风险利率
sigma = 0.2    # 波动率(收益标准差)
M = 50         # number of time steps
dt = T/M       # time enterval
I = 20000       # number of simulation
start = time()
S = []     #
for i in range(I):
path = []    # 时间间隔上的模拟路径
for t in range(M+1):
if t==0:
path.append(S_0)
else:
z = gauss(0.0, 1.0)
S_t = path[t-1] * exp((r-0.5*sigma**2) * dt + sigma * sqrt(dt) * z)
path.append(S_t)
S.append(path)
# 计算期权现值
C_0 = exp(-r * T) *sum([max(path[-1] -K, 0) for path in S])/I
total_time = time() - start
print 'European Option value %.6f'% C_0
print 'total time is %.6f seconds'% total_time

European Option value 8.159995
total time is 2.384639 seconds


# 选取部分模拟路径可视化

import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(10,7))
plt.grid(True)
plt.xlabel('Time step')
plt.ylabel('index level')
for i in range(30):
plt.plot(S[i])

方法三：Numpy

import numpy as np
from time import time
# 计算的一些初始值
S_0 = 100.0    # 股票或指数初始的价格;
K = 105        #  行权价格
T = 1.0        #  期权的到期年限(距离到期日时间间隔)
r = 0.05       #   无风险利率
sigma = 0.2    # 波动率(收益标准差)
M = 50         # number of time steps
dt = T/M       # time enterval
I = 20000       # number of simulation
# 20000条模拟路径，每条路径５０个时间步数
S = np.zeros((M+1, I))
S[0] = S_0
np.random.seed(2000)
start = time()
for t in range(1, M+1):
z = np.random.standard_normal(I)
S[t] = S[t-1] * np.exp((r- 0.5 * sigma **2)* dt + sigma * np.sqrt(dt)*z)
C_0 = np.exp(-r * T)* np.sum(np.maximum(S[-1] - K, 0))/I
end = time()
# 估值结果
print 'total time is %.6f seconds'%(end-start)
print 'European Option Value %.6f'%C_0

total time is 0.166926 seconds
European Option Value 7.993282


# 前２０条模拟路径
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(10,7))
plt.grid(True)
plt.xlabel('Time step')
plt.ylabel('index level')
for i in range(20):
plt.plot(S.T[i])

# 到期时所有模拟指数水平的频率直方图
%matplotlib inline
plt.hist(S[-1], bins=50)
plt.grid(True)
plt.xlabel('index level')
plt.ylabel('frequency')

# 模拟期权到期日的内在价值
%matplotlib inline
plt.hist(np.maximum(S[-1]-K, 0), bins=50)
plt.grid(True)
plt.xlabel('option inner value')
plt.ylabel('frequency')

方法四：Numpy

$log{S}_{t}=log{S}_{t-\mathrm{\Delta }t}+\left(r-\frac{1}{2}{\sigma }^{2}\right)\mathrm{\Delta }t+\sigma \sqrt{\mathrm{\Delta }t}{z}_{t}$$log S_t = log S_{t - \Delta t} + (r - \frac{1}{2} \sigma^2)\Delta t + \sigma \sqrt{\Delta t} z_t$

import numpy as np
from time import time
# 计算的一些初始值
S_0 = 100.0    # 股票或指数初始的价格;
K = 105        #  行权价格
T = 1.0        #  期权的到期年限(距离到期日时间间隔)
r = 0.05       #   无风险利率
sigma = 0.2    # 波动率(收益标准差)
M = 50         # number of time steps
dt = T/M       # time enterval
I = 20000       # number of simulation
np.random.seed(2000)
start = time()
# 生成一个随机变量的数组，M+1行，I列
# 同时计算出没一条路径，每一个时间点的指数水平的增量
# np.cumsum(axis=0)，在列的方向上进行累加得到每一个时间步数上的指数水平
S = S_0 * np.exp(np.cumsum((r - 0.5*sigma **2) *dt +sigma *np.sqrt(dt) *np.random.standard_normal((M+1, I)),axis=0))
S [0] = S_0
C_0 = np.exp(-r * T) * np.sum(np.maximum(S[-1] - K, 0))/I
end = time()
print 'toatl time is %.6f seconds'%(end-start)
print 'Europaen Option Value %.6f'%C_0

toatl time is 0.156061 seconds
Europaen Option Value 8.113643


# 前２０条模拟路径
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(10,7))
plt.grid(True)
plt.xlabel('Time step')
plt.ylabel('index level')
plt.plot(S[:,:20])


# 到期时所有模拟指数水平的频率直方图
import matplotlib.pyplot as plt
%matplotlib inline
plt.hist(S[-1], bins=50)
plt.grid(True)
plt.xlabel('index level')
plt.ylabel('frequency')

# 模拟期权到期日的内在价值
%matplotlib inline
plt.hist(np.maximum(S[-1]-K, 0), bins=50)
plt.grid(True)
plt.xlabel('option inner value')
plt.ylabel('frequency')

sum(S[-1] < K)    # 在两万次模拟中超过一万次到期期权内在价值为０
10748


结果对比

１．Scipy，　估值结果：8.021352，耗时：0.001010s
２．Python，估值结果：8.159995，耗时：2.384639s
３．Numpy，估值结果：7.993282，耗时：0.166926s
４．Numpy，估值结果：8.113643，耗时：0.156061s
Scipy，用时最短是因为，沒有进行20000次的模拟估值．其他三个方法进行了20000次的模拟，基于Numpy的计算方法速度比较快．

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客