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

版权声明:本文为博主原创文章,未经博主允许不得转载。

蒙特卡洛模拟计算看涨期权

蒙特卡洛估值是一个计算量较高的一个算法,对于简单的问题也需要海量的计算,因此要高效的实现蒙特卡洛算法.
下面使用不同的方法实现蒙特卡洛算法
1.Scipy
2.纯Python
3.Numpy(1)
4.Numpy(2)

方法一:Scipy

欧式看涨期权的定价公式Black-Scholes-Merton(1973):
  公式1:      
        C(St,K,t,T,r,σ)=StN(d1)er(Tt)KN(d2)
       
        N(d)=12πde12x2dx
         
         d1=log(StK)+(r+σ22)(Tt)σTt

         d2=log(StK)+(rσ22)(Tt)σTt

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

# 倒入用到的库
from math import log, sqrt, exp
from scipy import stats
# 期权的定价计算,根据公式1.
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)表示的布朗运动.
随机微分方程:
 
        dSt=rStdt+σStdZt

Zt是一个服从布朗运动的随机变量 
下面是蒙特卡洛估值公式2,中的参数与公式1中的定义相同:

      St=StΔtexpr12σ2Δt+σΔtzt
 
变量 zt 是服从正态分布的随机变量,Δt是一个足够小的一个时间间隔.期权的到期年限T满足0<5<=T(Hilpisch的著作中)
 
蒙特卡洛数值化计算方法:

1.将到期日之前的时间间隔[0,T],分割为多个等距的Δt子间隔
2.进行 I 次的模拟,i {1,2,…I}
 每次模拟,循环M个子间隔, t {Δt,2Δt,3Δt…T}
 (1),每个间隔点,取一个随机数 zt(i)
 (2),根据离散化公式3,计算 ST
 (3),计算T时刻,期权的内在价值 hT:
    
             hT(ST(i))=max(ST(i)K,0)

3.根据I次的模拟,计算期权到期价值:

              C0=erT1Ti=0hT(ST(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

前30条模拟路径

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

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

使用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条模拟路径,每条路径50个时间步数
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

前20条模拟路径

# 前20条模拟路径
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

对数欧拉方法计算蒙特卡洛估值(不使用任何循环实现蒙特卡洛估值计算)
公式3的对数版本:
    
      logSt=logStΔt+(r12σ2)Δt+σΔtzt

才用递增方法,计算每个时间步数,指数水平增长的百分比(比例),最后在乘以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

前20条模拟路径

# 前20条模拟路径
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)    # 在两万次模拟中超过一万次到期期权内在价值为0
10748

结果对比

1.Scipy, 估值结果:8.021352,耗时:0.001010s
2.Python,估值结果:8.159995,耗时:2.384639s
3.Numpy,估值结果:7.993282,耗时:0.166926s
4.Numpy,估值结果:8.113643,耗时:0.156061s
Scipy,用时最短是因为,沒有进行20000次的模拟估值.其他三个方法进行了20000次的模拟,基于Numpy的计算方法速度比较快.

发布了209 篇原创文章 · 获赞 217 · 访问量 49万+
展开阅读全文

没有更多推荐了,返回首页

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

分享到微信朋友圈

×

扫一扫,手机浏览