蒙特卡洛模拟

蒙特卡洛模拟是一类通过随机采样来进行数值计算的方法，其在物理、金融等领域具有非常广泛的应用。

其中的典型问题是计算某个复杂的随机变量的期望：

很多重要的计算问题可以化归为这个计算期望的问题。

举个例子。假如我们希望计算一个 -维函数 的积分:

被积函数可以理解为 ，其中 为 个独立的均匀分布随机变量 在其取值空间 的概率密度函数。

那么我们就有：

对于期权定价而言，期权价格即为折现回报在风险中性测度下的期望。这种情况下，在 这个问题中， 可以理解为折现回报。

解决这个问题的一种自然的思路是：利用样本平均来近似数学期望。即采样随机变量 的 个独立同分布样本 ，然后计算样本均值：

根据大数定律，如果 的方差有界，那么 依概率收敛到 ，即对于任意的 :

那么收敛的速度是多少呢？这个问题可以通过中心极限定理解答。根据中心极限定理，当 趋于无穷大时，

趋于标准正态分布，其中 为 的标准差。

从这个结果可以看出来:

其中 为一个标准正态分布的随机变量。

所以可以说蒙特卡洛模拟的收敛阶数为 。

通过这个结果，我们也可以构建 的置信区间。

令 为标准正态分布的累积分布函数， 为其逆函数。那么 以概率 落在区间

那么，近似地， 以概率 落在区间

一般我们取 置信区间，即 ，这时

那么 置信区间为：

由于 一般比较复杂， 也常需要用样本估计：

接下来以Black-Scholes模型下的欧式看涨期权定价为例对以上过程进行编程。

import numpy as np
from scipy.stats import norm

# 欧式看涨期权参数
r, sig, S0, K, T = 0.02, 0.3, 100, 100, 1.0

# 准确值
d1 = (np.log(S0/K) + (r+sig**2/2)*T)/sig/np.sqrt(T)
d2 = d1 - sig*np.sqrt(T)
bs = S0*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2) 

# 样本的次数
ns = 100*np.arange(100, 5001, 100)

mc = np.zeros(len(ns)) # 蒙特卡洛模拟的估计值
cf_lower = np.zeros(len(ns)) # 95%置信区间下边界
cf_upper = np.zeros(len(ns)) # 95%置信区间上边界

for i in range(len(ns)):
    n = ns[i]
    
    # 到期标的资产价格的n个样本
    ST = S0*np.exp((r-sig**2/2)*T + sig*np.sqrt(T)*norm.rvs(size=n))
    
    # 到期回报的样本
    Y = np.exp(-r*T)*np.maximum(ST-K, 0)
    
    # 估计值
    mc[i] = np.mean(Y)
    
    # 标准差估计值
    sd = np.sqrt(np.sum((Y - mc[i])**2)/(n-1))
    
    # 置信区间的上下界估计
    cf_lower[i] = mc[i] - 1.96*sd/np.sqrt(n)
    cf_upper[i] = mc[i] + 1.96*sd/np.sqrt(n)

# 作图
from matplotlib import pyplot as plt

plt.plot(ns, mc, label="Monte Caro Estimator")
plt.plot(ns, bs*np.ones_like(ns), label="Exact Value")
# 阴影部分为置信区间
plt.fill_between(ns, cf_lower, cf_upper, color='b', alpha=0.1) 
plt.legend()
plt.xlabel('$n$')
plt.ylabel('Price')
plt.savefig('Monte_Carlo_BS.jpg')
plt.show()