版权声明：本文为博主原创文章，未经博主允许不得转载。
蒙特卡洛模拟计算看涨期权
蒙特卡洛估值是一个计算量较高的一个算法，对于简单的问题也需要海量的计算，因此要高效的实现蒙特卡洛算法．
下面使用不同的方法实现蒙特卡洛算法
１．Scipy
２．纯Python
３．Numpy(1)
４．Numpy(2)
方法一：Scipy
欧式看涨期权的定价公式Black-Scholes-Merton(1973):
公式１：
C
(
S
t
,
K
,
t
,
T
,
r
,
σ
)
=
S
t
∗
N
(
d
1
)
−
e
−
r
(
T
−
t
)
∗
K
∗
N
(
d
2
)
C(S_t,K,t,T,r,\sigma) = S_t * N(d_1) - e^{-r(T-t)}*K*N(d_2)
C(St,K,t,T,r,σ)=St∗N(d1)−e−r(T−t)∗K∗N(d2)
N
(
d
)
=
1
2
π
∫
−
∞
d
e
−
1
2
x
2
d
x
N(d) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{d}e^{-\frac{1}{2}x^2}dx
N(d)=2π1∫−∞de−21x2dx
d
1
=
l
o
g
(
S
t
K
)
+
(
r
+
σ
2
2
)
(
T
−
t
)
σ
T
−
t
d_1 = \frac{log(\frac{S_t}{K})+(r+\frac{\sigma^2}{2})(T-t)}{\sigma\sqrt{T-t}}
d1=σT−tlog(KSt)+(r+2σ2)(T−t)
d 2 = l o g ( S t K ) + ( r － σ 2 2 ) ( T − t ) σ T − t d_2 = \frac{log(\frac{S_t}{K})+(r－\frac{\sigma^2}{2})(T-t)}{\sigma\sqrt{T-t}} d2=σT−tlog(KSt)+(r－2σ2)(T−t)
S
t
S_t
St, t时刻股票或者指数的水平
K
K
K, 行权价格
T
T
T, 期权的到期年限(距离到期日时间间隔)
r
r
r, 无风险利率
σ
\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
下面的计算仍基于BSM(balck-scholes-merton),模型中的高风险标识(股票指数)在风险中立的情况下遵循以随机微分方程(SDE)表示的布朗运动．
随机微分方程：
d
S
t
=
r
S
t
d
t
+
σ
S
t
d
Z
t
dS_t = rS_tdt+\sigma S_tdZ_t
dSt=rStdt+σStdZt
Z
t
Z_t
Zt是一个服从布朗运动的随机变量
下面是蒙特卡洛估值公式２,中的参数与公式１中的定义相同：
$ 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$
变量
z
t
z_t
zt 是服从正态分布的随机变量，
Δ
t
\Delta t
Δt是一个足够小的一个时间间隔．期权的到期年限T满足0<5<=T(Hilpisch的著作中)
蒙特卡洛数值化计算方法：
１．将到期日之前的时间间隔[0,T],分割为多个等距的
Δ
t
\Delta t
Δt子间隔
２．进行 Ｉ 次的模拟，$ i \in$ {１，２，…Ｉ}
每次模拟，循环M个子间隔, $t \in $ {
Δ
t
,
2
Δ
t
,
3
Δ
t
\Delta t,2\Delta t,3 \Delta t
Δt,2Δt,3Δt…T}
(1),每个间隔点，取一个随机数
z
t
(
i
)
z_t(i)
zt(i)
(2),根据离散化公式３，计算 $ S_T$
(3),计算T时刻，期权的内在价值 $ h_T
:
∗
∗
: **
: ∗∗ h_T(S_T(i)) = max(S_T(i) - K, 0)$**
３．根据Ｉ次的模拟，计算期权到期价值：
$ 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])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-fQVCcGNh-1680055513923)(output_14_0.png)]
方法三：Numpy
使用numpy的一些数组计算，减少for循环使用
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])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ziysETVD-1680055513924)(output_20_0.png)]
到期指数模拟水平
# 到期时所有模拟指数水平的频率直方图
%matplotlib inline
plt.hist(S[-1], bins=50)
plt.grid(True)
plt.xlabel('index level')
plt.ylabel('frequency')
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1U3vhyjx-1680055513924)(output_21_1.png)]
到期期权内在价值
# 模拟期权到期日的内在价值
%matplotlib inline
plt.hist(np.maximum(S[-1]-K, 0), bins=50)
plt.grid(True)
plt.xlabel('option inner value')
plt.ylabel('frequency')
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-L5MZi1jy-1680055513925)(output_22_1.png)]
方法四：Numpy
对数欧拉方法计算蒙特卡洛估值(不使用任何循环实现蒙特卡洛估值计算)
公式３的对数版本:
$ log S_t = log S_{t - \Delta t} + (r - \frac{1}{2} \sigma^2)\Delta t + \sigma \sqrt{\Delta t} z_t $
才用递增方法，计算每个时间步数，指数水平增长的百分比(比例），最后在乘以 S 0 S_0 S0
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])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-HhGiR7jA-1680055513925)(output_28_1.png)]
模拟到期指数水平
# 到期时所有模拟指数水平的频率直方图
import matplotlib.pyplot as plt
%matplotlib inline
plt.hist(S[-1], bins=50)
plt.grid(True)
plt.xlabel('index level')
plt.ylabel('frequency')
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-SYINIDdQ-1680055513926)(output_29_1.png)]
模拟到期期权内在价值
# 模拟期权到期日的内在价值
%matplotlib inline
plt.hist(np.maximum(S[-1]-K, 0), bins=50)
plt.grid(True)
plt.xlabel('option inner value')
plt.ylabel('frequency')
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WbP8w4FV-1680055513927)(output_30_1.png)]
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的计算方法速度比较快．
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
