python模拟B-S期权定价模型

随机变量模拟

在当前股票指数水平S0给定的情况下，未来某个日期T的股票指数水平ST可以根据Black-Scholes-Merton公式模拟：
$S_{T}=S_{0}*exp((r-0.5 \sigma^{^{2}})T+\sigma \sqrt{T}z)$

r ：无风险利率
σ ：S的恒定波动率（= 收益率的标准差）
z：标准正态分布随机变量

#通过 standard_normal模拟几何布朗运动，随机变量呈对数正态分布
S0 = 100
r = 0.05
sigma = 0.25
T = 2.0
I = 10000
ST1 = S0 * np.exp((r - 0.5 * sigma ** 2) * T +sigma * math.sqrt(T) * npr.standard_normal(I))  #如果使用 lognormal 函数直接得出随机变量值，必须向函数提供均值和标准差，通过 lognormal模拟几何布朗运动：
ST2 = S0 * npr.lognormal((r - 0.5 * sigma ** 2) * T,sigma * math.sqrt(T), size=I)  #使用 scipy.stats子库和 print_statistics 比较模拟结果的分布特性
import scipy.stats as scs
def print_statistics(a1, a2):"""Print selected statistics"""sta1 = scs.describe(a1)sta2 = scs.describe(a2)print('%14s %14s %14s' % ('statistic', 'data set 1', 'data set 2'))print(45 * '-')print('%14s %14.3f %14.3f' % ('size', sta1[0], sta2[0]))print('%14s %14.3f %14.3f' % ('min', sta1[1][0], sta2[1][0]))print('%14s %14.3f %14.3f' % ('max', sta1[1][1], sta2[1][1]))print('%14s %14.3f %14.3f' % ('mean', sta1[2], sta2[2]))print('%14s %14.3f %14.3f' % ('std', np.sqrt(sta1[3]), np.sqrt(sta2[3])))print('%14s %14.3f %14.3f' % ('skew', sta1[4], sta2[4]))print('%14s %14.3f %14.3f' % ('kurtosis', sta1[5], sta2[5]))
print_statistics(ST1, ST2)
#差异主要是由于模拟中的所谓采样误差

  statistic     data set 1     data set 2
---------------------------------------------size      10000.000      10000.000min         23.167         27.035max        377.909        433.432mean        110.586        110.652std         40.483         40.817skew          1.090          1.215kurtosis          1.888          2.663

动态形式--几何布朗运动

Black-Scholes-Merton模型的动态形式由下面的随机微分方程（SDE）描述。Zt是标准布朗运动，SDE 被称作几何布朗运动。St的值呈对数正态分布，边际收益 dSt/St 呈正态分布。

公式 Black-Scholes-Merton设置中的随机微分方程
$dS_t=rS_td_t+\sigma S_tdZ_t$

SDE可以由一个欧拉格式精确地离散化，Δt是固定的离散化间隔， zt是标准正态分布随机变量，在Black-Scholes-Merton中动态模拟指数水平：

$S_t=S_{t-\Delta t}exp\left( \left( r-\frac{1}{2}\sigma^2 \right)\Delta t + \sigma \sqrt{\Delta t}Z_t \right)$

I = 10000
M = 50
dt = T / M
S = np.zeros((M + 1, I))
S[0] = S0
for t in range(1, M + 1):S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt +sigma * math.sqrt(dt) * npr.standard_normal(I))  #展示前10条模拟路径：
plt.figure(figsize=(10, 6))
plt.plot(S[:, :10], lw=1.5)
plt.xlabel('time')
plt.ylabel('index level')

均值回归过程--平方根扩散

均值回归过程用于建立短期利率或者波动性过程模型。平方根扩散对应的随机微分方程：

$dx_t=k(\theta-x_t)dt+\sigma \sqrt{x_t}dZ_t$

xt ：日期 t 的过程水平
k ：均值回归因子
θ ：长期过程均值
σ ：恒定波动率
Z ：标准布朗运动
xt的值可以使用正态分布进行离散化和近似计算。平方根扩散的欧拉离散化：

$\tilde{x}_t=\tilde{x}_s+k(\theta- \tilde{x}_s^+ )\Delta t+\sigma \sqrt{ \tilde{x}_s^+ } \sqrt{ \Delta t z_t }$ ， $x_t=\tilde{x}_t^+$ 、

到期日的模拟平方根扩散：

#用表示短期利率模型的值参数化后续模拟所用的模型
x0 = 0.05
kappa = 3.0
theta = 0.02
sigma = 0.1
I = 10000
M = 50
dt = T / Mdef srd_euler():xh = np.zeros((M + 1, I))x1 = np.zeros_like(xh)xh[0] = x0x1[0] = x0for t in range(1, M + 1):xh[t] = (xh[t - 1]+ kappa * (theta - np.maximum(xh[t - 1], 0)) * dt+ sigma * np.sqrt(np.maximum(xh[t - 1], 0)) * np.sqrt(dt)* random.standard_normal(I))x1 = np.maximum(xh, 0)return x1x1 = srd_euler()
#前十条路径，说明得出的平均偏离值为负值(x0>θ=0.02)，收敛于 θ
plt.plot(x1[:, :10], lw=1.5)
plt.xlabel('time')
plt.ylabel('index level')
plt.grid(True)

波动率随机的情况

B-S模型假设波动率恒定。而波动率是随机的。 Heston随机波动率模型的随机徽分方程：

$dS_t=rS_tdt+\sqrt{v_t}S_tdZ_t^1$

$dv_t=k_v(\theta_v-v_t)dt+\sigma_v\sqrt{v_t}dZ_t^2$

$dZ_t^1dZ_t^2=\rho$

ρ 代表两个标准布朗运动 Zt1 ，Zt2之间的瞬时相关性。用于解释杠杆效应：波动性在衰退市场中上升，而在牛市下降。

S0 = 100
r = 0.05
v0 = 0.1
kappa = 3.0
theta = 0.25
sigma = 0.1
rho = 0.6
T = 1.0
#确定两个随机过程相关矩阵的柯列斯基分解：
corr_mat = np.zeros((2, 2))
corr_mat[0, :] = [1.0, rho]
corr_mat[1, :] = [rho, 1.0]
cho_mat = np.linalg.cholesky(corr_mat)
cho_mat
# array([[ 1. ,  0. ],
#        [ 0.6,  0.8]])
# 为指数过程和波动性过程生成整组随机数
M = 50
I = 10000
ran_num = random.standard_normal((2, M + 1, I))
# 以平方根扩散过程类型建模波动性过程
dt = T / M
v = np.zeros_like(ran_num[0])
vh = np.zeros_like(v)
v[0] = v0
vh[0] = v0
for t in range(1, M + 1):ran = np.dot(cho_mat, ran_num[:, t, :])vh[t] = (vh[t - 1] + kappa * (theta - np.maximum(vh[t - 1], 0)) * dt+ sigma * np.sqrt(np.maximum(vh[t - 1], 0)) * np.sqrt(dt)* ran[1])v = np.maximum(vh, 0)
#使用几何布朗运动的欧拉格式建模指数水平过程
S = np.zeros_like(ran_num[0])
S[0] = S0
for t in range(1, M + 1):ran = np.dot(cho_mat, ran_num[:, t, :])S[t] = S[t - 1] * np.exp((r - 0.5 * v[t]) * dt +np.sqrt(v[t]) * ran[0] * np.sqrt(dt))#指数水平过程和波动性过程的模拟结果
#到期日的模拟随机波动率模型
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 5))
ax1.hist(S[-1], bins=50)
ax1.set_xlabel('index level')
ax1.set_ylabel('frequency')
ax1.grid(True)
ax2.hist(v[-1], bins=50)
ax2.set_xlabel('volatility')
ax2.grid(True)

#模拟随机波动率模型路径
fig, (ax1, ax2) = plt.subplots(2,1,sharex=True, figsize=(7, 6))
ax1.plot(S[:,:10],lw=1.5)
ax1.set_ylabel('index level')
ax1.grid(True)
ax2.plot(v[:,:10],lw=1.5)
ax2.set_xlabel('time')
ax2.set_ylabel('volatility')
ax2.grid(True)
#每个过程的前10条模拟路径表明波动性过程的平均漂移值为正数，收敛于 θv=0.25