金融工程里常用：指用随机过程（如布朗运动、随机波动率、跳跃过程）来描述资产价格、利率、波动率等不确定变量的建模方法。def H93_char_func(u, T, r, kappa_v, theta_v, sigma_v, rho, v0): """Heston模型特征函数（Lewis 2001） 参数：与上述5个待校准参数一致，新增r（无风险利率）、T（到期时间） 返回：复数类型的特征函数值 """ c1 = kappa_v * theta_v c2 = -np.sqrt( (rho * sigma_v * u * 1j - kappa_v) ** 2 - sigma_v**2 * (-u * 1j - u**2) ) c3 = (kappa_v - rho * sigma_v * u * 1j + c2) / ( kappa_v - rho * sigma_v * u * 1j - c2 ) H1 = r * u * 1j * T + (c1 / sigma_v**2) * ( (kappa_v - rho * sigma_v * u * 1j + c2) * T - 2 * np.log((1 - c3 * np.exp(c2 * T)) / (1 - c3)) ) H2 = ( (kappa_v - rho * sigma_v * u * 1j + c2) / sigma_v**2 * ((1 - np.exp(c2 * T)) / (1 - c3 * np.exp(c2 * T))) ) char_func_value = np.exp(H1 + H2 * v0) return char_func_value
2.3 数值积分（scipy.integrate.quad 实操）
2.3.1 积分核函数实现
def H93_int_func(u, S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0): """Lewis(2001)定价公式中的积分核函数 作用：构造quad函数需要积分的目标函数 """ # 调用特征函数，注意传入u-1j*0.5（公式要求） char_func_value = H93_char_func(u - 1j * 0.5, T, r, kappa_v, theta_v, sigma_v, rho, v0) # 计算积分核的实部（积分仅需实部） int_func_value = (1 / (u**2 + 0.25) * (np.exp(1j * u * np.log(S0 / K)) * char_func_value).real) return int_func_value
2.3.2 quad函数关键用法（结合文档）
核心功能：计算定积分（支持无穷区间），对应定价公式中的积分计算，关键参数说明：
func：需积分的函数（此处为H93_int_func）
a, b：积分上下限（此处a=0，b=np.inf，即0到无穷大）
limit：自适应算法的子区间上限（示例中设为250，避免积分精度不足）
返回值：[0]为积分结果，[1]为绝对误差估计（实操中仅需取[0]）
2.3.3 看涨期权价格计算（整合特征函数+积分）
def H93_call_value(S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0): """Heston模型看涨期权定价（Lewis 2001） 参数： S0：初始标的价格；K：行权价；T：到期时间（年）；r：无风险利率 其余5个为Heston模型参数 返回：看涨期权价格（float） """ # 调用quad计算积分，lambda函数传递额外参数 int_value = quad( lambda u: H93_int_func(u, S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0), 0, np.inf, limit=250, # 提高积分精度，避免收敛误差 )[0] # 代入定价公式，max(0, ...)避免负价格（边界情况） call_value = max(0, S0 - np.exp(-r * T) * np.sqrt(S0 * K) / np.pi * int_value) return call_value
import pandas as pdimport numpy as np# 加载数据（HDF5格式）h5 = pd.HDFStore("option_data_wqu.h5", "r")data = h5["data"] # 包含欧式看涨/看跌期权数据（3个到期日）h5.close()S0 = 3225.93 # 标的初始价格（EuroStoxx 50指数）# 筛选ATM期权tol = 0.02 # 2%容忍度options = data[(np.abs(data["Strike"] - S0) / S0) < tol].copy()options["Date"] = pd.to_datetime(options["Date"])options["Maturity"] = pd.to_datetime(options["Maturity"])# 计算到期时间T（年）和无风险利率rfor row, option in options.iterrows(): options.loc[row, "T"] = (option["Maturity"] - option["Date"]).days / 365.0 options.loc[row, "r"] = 0.02 # 设定常数无风险利率
步骤2：定义误差函数（目标函数）
核心：计算模型价格与市场价格的MSE，加入参数约束（避免无意义参数值）
i = 0min_MSE = 500 # 初始最小MSE（用于跟踪优化过程）def H93_error_function(p0): """校准的误差函数（MSE） 参数p0：待校准的5个参数，顺序为[kappa_v, theta_v, sigma_v, rho, v0] 返回：模型与市场价格的均方误差（MSE） """ global i, min_MSE kappa_v, theta_v, sigma_v, rho, v0 = p0 # 参数约束（避免无意义值，不符合约束则返回大MSE） if kappa_v < 0.0 or theta_v < 0.005 or sigma_v < 0.0 or rho < -1.0 or rho > 1.0: return 500.0 # Feller条件（保证方差为正）：2*kappa_v*theta_v >= sigma_v^2 if 2 * kappa_v * theta_v < sigma_v**2: return 500.0 # 计算每个期权的模型价格与市场价格的平方误差 se = [] # 平方误差列表 for row, option in options.iterrows(): model_value = H93_call_value( S0, option["Strike"], option["T"], option["r"], kappa_v, theta_v, sigma_v, rho, v0 ) se.append((model_value - option["Call"]) ** 2) # 市场价格为option["Call"] MSE = sum(se) / len(se) # 均方误差 min_MSE = min(min_MSE, MSE) # 更新最小MSE # 每25次迭代打印一次结果（跟踪优化进度） if i % 25 == 0: print("%4d |" % i, np.array(p0), "| %7.3f | %7.3f" % (MSE, min_MSE)) i += 1 return MSE
步骤3：两步优化校准（scipy.brute + scipy.optimize.fmin）
核心逻辑：先粗搜（brute）确定参数合理区间，再精调（fmin）获得最优参数，提升收敛速度和精度（结合文档中两个函数的用法）。
3.3.1 scipy.optimize.brute 用法（粗搜）
功能：暴力搜索，在指定参数范围内构建网格，计算每个网格点的误差函数值，找到全局最优网格点（避免局部最优）。
关键参数（结合文档）：
func：误差函数（H93_error_function）
ranges：每个参数的搜索范围（元组形式，(low, high, step)）
finish：粗搜后是否精调（此处设为None，仅做粗搜）
3.3.2 scipy.optimize.fmin 用法（精调）
功能：使用下山单纯形法（Nelder-Mead），以粗搜结果为初始值，局部精细优化，提升参数精度。
关键参数（结合文档）：
3.3.3 校准函数实现
from scipy.optimize import brute, fmindef H93_calibration_full(): """Heston模型完整校准函数（粗搜+精调） 返回：校准后的最优参数（array类型，顺序为[kappa_v, theta_v, sigma_v, rho, v0]） """ # 第一步：暴力粗搜（确定参数合理区间） p0 = brute( H93_error_function, ( (2.5, 10.6, 5.0), # kappa_v 搜索范围 (0.01, 0.041, 0.01), # theta_v 搜索范围 (0.05, 0.251, 0.1), # sigma_v 搜索范围 (-0.75, 0.01, 0.25), # rho 搜索范围 (0.01, 0.031, 0.01), # v0 搜索范围 ), finish=None, # 暂不精调，仅粗搜 ) # 第二步：局部精调（以粗搜结果为初始值） opt = fmin( H93_error_function, p0, xtol=0.000001, ftol=0.000001, # 收敛精度 maxiter=750, maxfun=900 # 最大迭代/函数调用次数 ) return opt
步骤4：运行校准，获取结果
# 调用校准函数，得到最优参数calibrated_params = H93_calibration_full()print("校准后的Heston模型参数：")print(f"kappa_v = {calibrated_params[0]:.3f}")print(f"theta_v = {calibrated_params[1]:.3f}")print(f"sigma_v = {calibrated_params[2]:.3f}")print(f"rho = {calibrated_params[3]:.3f}")print(f"v0 = {calibrated_params[4]:.3f}")
四、补充任务：实现看跌期权定价（平价公式）
4.1 核心公式（无股息看涨看跌平价）
def H93_put_value_parity(S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0): """ 用看涨看跌平价公式计算Heston模型看跌期权价格 参数：与H93_call_value一致 返回：看跌期权价格（float，负价格钳位为0） """ # 调用已实现的看涨期权定价函数 call_price = H93_call_value(S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0) # 用平价公式计算看跌价格 put_price = call_price - S0 + K * np.exp(-r * T) # 钳位负价格（避免舍入误差导致的微小负值） return max(0.0, put_price)
五、关键注意事项（避坑指南）
参数约束：必须满足Feller条件（$$2\kappa_v\theta_v \geq \sigma_v^2$$），否则方差可能为负，模型无意义。
积分精度：quad函数的limit参数需适当调大（如250），避免无穷区间积分的收敛误差。
数据筛选：优先选择ATM期权，远离虚值/实值期权（流动性差，价格偏离真实预期）。
优化效率：brute的网格步长不宜过小（避免计算量过大），fmin的收敛容忍度不宜过高（避免优化超时）。
实务建议：校准应结合看涨+看跌期权数据，避免看跌-看涨平价不完全成立导致的偏差。
六、核心总结
Heston校准流程：定价函数（Lewis+quad）→ 误差函数（MSE）→ 两步优化（brute粗搜+fmin精调）→ 最优参数。
scipy函数核心作用：quad负责数值积分，brute负责全局粗搜，fmin负责局部精调，三者协同完成校准。
核心目标：让模型价格贴合市场，为奇异期权定价和对冲提供可靠参数。
LESSON 2 - 模型校准：默顿模型（1976年）二、代码实现：从特征函数到期权定价
核心步骤：先定义特征函数 → 定义积分函数 → 求解看涨期权价格，全程基于Python实现，依赖numpy、scipy（积分求解）。
2.1 特征函数实现
import numpy as npfrom scipy.integrate import quaddef M76_char_func(u, T, r, sigma, lamb, mu, delta): """ 默顿（1976）模型的特征函数 参数： u: 特征函数的输入参数 T: 到期时间 r: 无风险利率 sigma: 扩散项波动率 lamb: 跳跃强度 mu: 跳跃幅度期望 delta: 跳跃幅度标准差 返回： char_func_value: 特征函数值（复数） """ omega = r - 0.5 * sigma**2 - lamb * (np.exp(mu + 0.5 * delta**2) - 1) char_func_value = np.exp( ( 1j * u * omega # 1j表示复数单位i - 0.5 * u**2 * sigma**2 + lamb * (np.exp(1j * u * mu - u**2 * delta**2 * 0.5) - 1) ) * T ) return char_func_value
2.2 积分函数实现（Lewis公式核心）
def M76_integration_function(u, S0, K, T, r, sigma, lamb, mu, delta): """ Lewis（2001）定价公式中的积分项 功能：计算积分被积函数的值（取实部，确保积分可解） """ # 特征函数调整：u替换为u - 0.5*1j（对应公式中的z-i/2） char_func = M76_char_func(u - 0.5 * 1j, T, r, sigma, lamb, mu, delta) # 被积函数：实部处理（避免复数积分） value = 1 / (u**2 + 0.25) * (np.exp(1j * u * np.log(S0 / K)) * char_func).real return value
2.3 看涨期权价格计算
def M76_call_value(S0, K, T, r, sigma, lamb, mu, delta): """ 默顿（1976）模型下看涨期权价格（基于Lewis（2001）方法） 步骤：1. 数值积分求解Lewis公式中的积分项；2. 代入公式计算期权价格 """ # 数值积分：积分区间[0,50]，limit提高积分精度 int_value = quad( lambda u: M76_integration_function(u, S0, K, T, r, sigma, lamb, mu, delta), 0, 50, limit=250, )[0] # quad返回(积分值, 误差)，取第一个元素 # 代入Lewis公式计算期权价格，max(0,·)确保期权价格非负 call_value = max(0, S0 - np.exp(-r * T) * np.sqrt(S0 * K) / np.pi * int_value) return call_value
import pandas as pd# 加载市场数据（HDF5格式）h5 = pd.HDFStore("option_data_wqu.h5", "r")data = h5["data"] # 包含欧式看涨/看跌期权数据h5.close()S0 = 3225.93 # 2014年9月30日EuroStoxx 50指数水平# 筛选近ATM期权tol = 0.02options = data[(np.abs(data["Strike"] - S0) / S0) < tol].copy()options["Date"] = pd.to_datetime(options["Date"])options["Maturity"] = pd.to_datetime(options["Maturity"])# 补充：计算到期剩余时间T（年）、设定无风险利率r（常数）options["T"] = (options["Maturity"] - options["Date"]).dt.days / 365options["r"] = 0.05 # 假设无风险利率为5%
def M76_error_function(p0): """ 默顿模型校准的误差函数（RMSE） 参数p0：待校准参数列表 [sigma, lamb, mu, delta] 返回：RMSE值（越小越好） """ global i, min_RMSE # 全局变量，用于跟踪迭代过程中的最小RMSE sigma, lamb, mu, delta = p0 # 约束条件：参数非负（波动率、跳跃强度、跳跃标准差不能为负） if sigma < 0.0 or delta < 0.0 or lamb < 0.0: return 500.0 # 违反约束，返回大误差 # 计算每个期权的模型价格与市场价格的平方误差 se = [] for row, option in options.iterrows(): model_value = M76_call_value( S0, option["Strike"], option["T"], option["r"], sigma, lamb, mu, delta ) se.append((model_value - option["Call"]) ** 2) # 平方误差 RMSE = np.sqrt(sum(se) / len(se)) # 均方根误差 min_RMSE = min(min_RMSE, RMSE) # 更新最小RMSE # 每50次迭代打印一次结果，便于跟踪 if i % 50 == 0: print("%4d |" % i, np.array(p0), "| %7.3f | %7.3f" % (RMSE, min_RMSE)) i += 1 return RMSE
3.3 步骤3：校准优化（两阶段优化）
采用“暴力搜索+局部优化”结合的方式，确保优化收敛到全局最优（避免局部最优解）：
暴力搜索（brute）：在合理参数区间内扫描，找到较优的初始参数（p0）；
局部优化（fmin）：以暴力搜索结果为初始值，进行精细优化，提高校准精度。
def M76_calibration_full(): """ 默顿（1976）模型完整校准函数 返回：校准后的最优参数 [sigma, lamb, mu, delta] """ # 阶段1：暴力搜索（参数区间需合理设定，基于金融常识） p0 = brute( M76_error_function, ( (0.075, 0.201, 0.025), # sigma区间：0.075~0.201，步长0.025 (0.10, 0.401, 0.1), # lambda区间：0.10~0.401，步长0.1 (-0.5, 0.01, 0.1), # mu区间：-0.5~0.01，步长0.1 (0.10, 0.301, 0.1), # delta区间：0.10~0.301，步长0.1 ), finish=None, # 不进行后续局部优化，仅返回暴力搜索结果 ) # 阶段2：局部凸优化（精细调整参数） opt = fmin( M76_error_function, p0, xtol=0.0001, ftol=0.0001, maxiter=550, maxfun=1050 ) return opt# 初始化全局变量，开始校准i = 0min_RMSE = float("inf")opt = M76_calibration_full()
import matplotlib.pyplot as pltdef generate_plot(opt, options): sigma, lamb, mu, delta = opt # 计算所有期权的模型价格 options["Model"] = 0.0 for row, option in options.iterrows(): options.loc[row, "Model"] = M76_call_value( S0, option["Strike"], option["T"], option["r"], sigma, lamb, mu, delta ) # 按到期日分组绘图 mats = sorted(set(options["Maturity"])) options = options.set_index("Strike") for i, mat in enumerate(mats): options[options["Maturity"] == mat][["Call", "Model"]].plot( style=["b-", "ro"], title="到期日：%s" % str(mat)[:10] ) plt.ylabel("期权价格") plt.show()# 调用绘图函数generate_plot(opt, options)
理想拟合效果：红点（模型价格）与蓝线（市场价格）高度重合，说明模型能较好地解释市场期权价格。
LESSON 3 - 结合默顿模型（1976 年）与赫斯顿模型（1993 年）默顿模型（1976）与Lewis（2001）定价方法二、代码实现：从特征函数到期权定价
核心步骤：先定义特征函数 → 定义积分函数 → 求解看涨期权价格，全程基于Python实现，依赖numpy、scipy（积分求解）。
2.1 特征函数实现
import numpy as npfrom scipy.integrate import quaddef M76_char_func(u, T, r, sigma, lamb, mu, delta): """ 默顿（1976）模型的特征函数 参数： u: 特征函数的输入参数 T: 到期时间 r: 无风险利率 sigma: 扩散项波动率 lamb: 跳跃强度 mu: 跳跃幅度期望 delta: 跳跃幅度标准差 返回： char_func_value: 特征函数值（复数） """ omega = r - 0.5 * sigma**2 - lamb * (np.exp(mu + 0.5 * delta**2) - 1) char_func_value = np.exp( ( 1j * u * omega # 1j表示复数单位i - 0.5 * u**2 * sigma**2 + lamb * (np.exp(1j * u * mu - u**2 * delta**2 * 0.5) - 1) ) * T ) return char_func_value
def M76_integration_function(u, S0, K, T, r, sigma, lamb, mu, delta): """ Lewis（2001）定价公式中的积分项 功能：计算积分被积函数的值（取实部，确保积分可解） """ # 特征函数调整：u替换为u - 0.5*1j（对应公式中的z-i/2） char_func = M76_char_func(u - 0.5 * 1j, T, r, sigma, lamb, mu, delta) # 被积函数：实部处理（避免复数积分） value = 1 / (u**2 + 0.25) * (np.exp(1j * u * np.log(S0 / K)) * char_func).real return value
def M76_call_value(S0, K, T, r, sigma, lamb, mu, delta): """ 默顿（1976）模型下看涨期权价格（基于Lewis（2001）方法） 步骤：1. 数值积分求解Lewis公式中的积分项；2. 代入公式计算期权价格 """ # 数值积分：积分区间[0,50]，limit提高积分精度 int_value = quad( lambda u: M76_integration_function(u, S0, K, T, r, sigma, lamb, mu, delta), 0, 50, limit=250, )[0] # quad返回(积分值, 误差)，取第一个元素 # 代入Lewis公式计算期权价格，max(0,·)确保期权价格非负 call_value = max(0, S0 - np.exp(-r * T) * np.sqrt(S0 * K) / np.pi * int_value) return call_value
import pandas as pd# 加载市场数据（HDF5格式）h5 = pd.HDFStore("option_data_wqu.h5", "r")data = h5["data"] # 包含欧式看涨/看跌期权数据h5.close()S0 = 3225.93 # 2014年9月30日EuroStoxx 50指数水平# 筛选近ATM期权tol = 0.02options = data[(np.abs(data["Strike"] - S0) / S0) < tol].copy()options["Date"] = pd.to_datetime(options["Date"])options["Maturity"] = pd.to_datetime(options["Maturity"])# 补充：计算到期剩余时间T（年）、设定无风险利率r（常数）options["T"] = (options["Maturity"] - options["Date"]).dt.days / 365options["r"] = 0.05 # 假设无风险利率为5%
def M76_error_function(p0): """ 默顿模型校准的误差函数（RMSE） 参数p0：待校准参数列表 [sigma, lamb, mu, delta] 返回：RMSE值（越小越好） """ global i, min_RMSE # 全局变量，用于跟踪迭代过程中的最小RMSE sigma, lamb, mu, delta = p0 # 约束条件：参数非负（波动率、跳跃强度、跳跃标准差不能为负） if sigma < 0.0 or delta < 0.0 or lamb < 0.0: return 500.0 # 违反约束，返回大误差 # 计算每个期权的模型价格与市场价格的平方误差 se = [] for row, option in options.iterrows(): model_value = M76_call_value( S0, option["Strike"], option["T"], option["r"], sigma, lamb, mu, delta ) se.append((model_value - option["Call"]) ** 2) # 平方误差 RMSE = np.sqrt(sum(se) / len(se)) # 均方根误差 min_RMSE = min(min_RMSE, RMSE) # 更新最小RMSE # 每50次迭代打印一次结果，便于跟踪 if i % 50 == 0: print("%4d |" % i, np.array(p0), "| %7.3f | %7.3f" % (RMSE, min_RMSE)) i += 1 return RMSE
3.3 步骤3：校准优化（两阶段优化）
采用“暴力搜索+局部优化”结合的方式，确保优化收敛到全局最优（避免局部最优解）：
暴力搜索（brute）：在合理参数区间内扫描，找到较优的初始参数（p0）；
局部优化（fmin）：以暴力搜索结果为初始值，进行精细优化，提高校准精度。
def M76_calibration_full(): """ 默顿（1976）模型完整校准函数 返回：校准后的最优参数 [sigma, lamb, mu, delta] """ # 阶段1：暴力搜索（参数区间需合理设定，基于金融常识） p0 = brute( M76_error_function, ( (0.075, 0.201, 0.025), # sigma区间：0.075~0.201，步长0.025 (0.10, 0.401, 0.1), # lambda区间：0.10~0.401，步长0.1 (-0.5, 0.01, 0.1), # mu区间：-0.5~0.01，步长0.1 (0.10, 0.301, 0.1), # delta区间：0.10~0.301，步长0.1 ), finish=None, # 不进行后续局部优化，仅返回暴力搜索结果 ) # 阶段2：局部凸优化（精细调整参数） opt = fmin( M76_error_function, p0, xtol=0.0001, ftol=0.0001, maxiter=550, maxfun=1050 ) return opt# 初始化全局变量，开始校准i = 0min_RMSE = float("inf")opt = M76_calibration_full()
import matplotlib.pyplot as pltdef generate_plot(opt, options): sigma, lamb, mu, delta = opt # 计算所有期权的模型价格 options["Model"] = 0.0 for row, option in options.iterrows(): options.loc[row, "Model"] = M76_call_value( S0, option["Strike"], option["T"], option["r"], sigma, lamb, mu, delta ) # 按到期日分组绘图 mats = sorted(set(options["Maturity"])) options = options.set_index("Strike") for i, mat in enumerate(mats): options[options["Maturity"] == mat][["Call", "Model"]].plot( style=["b-", "ro"], title="到期日：%s" % str(mat)[:10] ) plt.ylabel("期权价格") plt.show()# 调用绘图函数generate_plot(opt, options)
理想拟合效果：红点（模型价格）与蓝线（市场价格）高度重合，说明模型能较好地解释市场期权价格。
4.1 核心要点
默顿（1976）模型：跳跃-扩散结合，弥补了几何布朗运动无法解释资产价格突发波动的缺陷；
Lewis（2001）方法：通过特征函数+数值积分定价，简化了跳跃模型的定价复杂度；
校准逻辑：两阶段优化（暴力搜索+局部优化），约束参数非负，最小化RMSE，确保拟合效果；
关键参数：$$\sigma$$（扩散波动）、$$\lambda$$（跳跃频率）、$$\mu$$（跳跃幅度）、$$\delta$$（跳跃波动），均需通过校准获得。
4.2 拓展练习
当前校准仅使用看涨期权数据，可修改代码，将看跌期权数据纳入校准（利用看涨-看跌平价公式验证，或直接加入误差计算），提升校准的准确性。
4.3 后续延伸：贝茨模型（Bates, 1996）
默顿模型仅考虑“常数波动率+泊松跳跃”，实际市场中波动率是随机的（如Heston模型），贝茨模型（1996）将赫斯顿模型（1993，随机波动率）与默顿模型（1976，跳跃-扩散）的优势结合，构建“随机波动率+跳跃-扩散”复合模型（Bates, 1996），进一步提升定价和校准的准确性，是金融衍生品定价中更贴近现实的经典模型。
五、贝茨模型（Bates, 1996）详细解析
5.1 模型核心逻辑
贝茨模型的核心思想简洁：将赫斯顿模型（1993）的随机波动率特征与默顿模型（1976）的跳跃-扩散特征融合，同时捕捉资产价格的两大关键特征——波动率的动态变化（时变波动）和价格的突发跳跃，解决了单一模型无法兼顾两者的缺陷。
模型假设在风险中性测度下，资产价格和波动率的动态过程如下：
5.2 关键公式与参数说明
5.3 代码实现：从特征函数到期权定价
贝茨模型的定价的核心是实现其特征函数，再结合此前学习的Lewis（2001）和Carr-Madan（1999，FFT方法）进行期权定价，全程基于Python，依赖numpy、scipy（积分）、matplotlib（可选，可视化）。
（1）赫斯顿模型特征函数实现
def H93_char_func(u, T, r, kappa_v, theta_v, sigma_v, rho, v0): """ 赫斯顿（1993）模型的特征函数 参数： u: 特征函数输入参数 T: 到期时间 r: 无风险利率 kappa_v: 方差调整速度 theta_v: 方差长期均值 sigma_v: 方差的波动率 rho: 资产价格与方差的相关性 v0: 初始方差 返回： char_func_value: 特征函数值（复数） """ c1 = kappa_v * theta_v c2 = -np.sqrt( (rho * sigma_v * u * 1j - kappa_v) ** 2 - sigma_v**2 * (-u * 1j - u**2) ) c3 = (kappa_v - rho * sigma_v * u * 1j + c2) / ( kappa_v - rho * sigma_v * u * 1j - c2 ) H1 = r * u * 1j * T + (c1 / sigma_v**2) * ( (kappa_v - rho * sigma_v * u * 1j + c2) * T - 2 * np.log((1 - c3 * np.exp(c2 * T)) / (1 - c3)) ) H2 = ( (kappa_v - rho * sigma_v * u * 1j + c2) / sigma_v**2 * ((1 - np.exp(c2 * T)) / (1 - c3 * np.exp(c2 * T))) ) char_func_value = np.exp(H1 + H2 * v0) return char_func_value
（2）修正后默顿模型特征函数（仅跳跃项）
def M76J_char_func(u, T, lamb, mu, delta): """ 修正后默顿（1976）模型特征函数：仅保留跳跃项 参数： u: 特征函数输入参数 T: 到期时间 lamb: 跳跃强度 mu: 跳跃幅度期望 delta: 跳跃幅度标准差 返回： char_func_value: 跳跃项特征函数值（复数） """ omega = -lamb * (np.exp(mu + 0.5 * delta**2) - 1) char_func_value = np.exp( (1j * u * omega + lamb * (np.exp(1j * u * mu - u**2 * delta**2 * 0.5) - 1)) * T ) return char_func_value
（3）贝茨模型特征函数（融合两者）
def B96_char_func(u, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta): """ 贝茨（1996）模型特征函数：赫斯顿特征函数 × 修正默顿跳跃项特征函数 参数： 包含赫斯顿和修正默顿模型的所有参数 返回： char_func_value: 贝茨模型特征函数值（复数） """ H93 = H93_char_func(u, T, r, kappa_v, theta_v, sigma_v, rho, v0) M76J = M76J_char_func(u, T, lamb, mu, delta) return H93 * M76J
def B96_int_func(u, S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta): """ Lewis（2001）方法的积分项：适配贝茨模型特征函数 """ # 特征函数调整：u替换为u - 0.5*1j（与Lewis公式一致） char_func_value = B96_char_func( u - 1j * 0.5, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta ) # 计算被积函数（取实部，确保积分可解） int_func_value = ( 1 / (u**2 + 0.25) * (np.exp(1j * u * np.log(S0 / K)) * char_func_value).real ) return int_func_valuedef B96_call_value_Lewis(S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta): """ 贝茨（1996）模型下看涨期权价格（基于Lewis（2001）方法） """ # 数值积分求解Lewis公式中的积分项 int_value = quad( lambda u: B96_int_func( u, S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta ), 0, np.inf, limit=250, )[0] # 代入公式计算期权价格，确保非负 call_value = max(0, S0 - np.exp(-r * T) * np.sqrt(S0 * K) / np.pi * int_value) return call_value
（5）基于Carr-Madan（1999）的FFT定价方法
作为Lewis方法的替代，FFT（快速傅里叶变换）可提高定价效率，核心是将期权价格积分转化为傅里叶变换求解，适配贝茨模型特征函数：
def B96_call_FFT(S0, K, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta): """ 贝茨（1996）模型下看涨期权价格（基于Carr-Madan（1999）FFT方法） """ k = np.log(K / S0) g = 1 # 精度调整因子 N = g * 4096 eps = (g * 150) ** -1 eta = 2 * np.pi / (N * eps) b = 0.5 * N * eps - k u = np.arange(1, N + 1, 1) vo = eta * (u - 1) # 根据期权虚实值调整参数，确保积分可积性 if S0 >= 0.95 * K: # 实值期权（ITM） alpha = 1.5 v = vo - (alpha + 1) * 1j modcharFunc = np.exp(-r * T) * ( B96_char_func(v, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta) / (alpha**2 + alpha - vo**2 + 1j * (2 * alpha + 1) * vo) ) else: # 虚值期权（OTM） alpha = 1.1 # 计算两个修正特征函数 v = (vo - 1j * alpha) - 1j modcharFunc1 = np.exp(-r * T) * ( 1 / (1 + 1j * (vo - 1j * alpha)) - np.exp(r * T) / (1j * (vo - 1j * alpha)) - B96_char_func( v, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta ) / ((vo - 1j * alpha) ** 2 - 1j * (vo - 1j * alpha)) ) v = (vo + 1j * alpha) - 1j modcharFunc2 = np.exp(-r * T) * ( 1 / (1 + 1j * (vo + 1j * alpha)) - np.exp(r * T) / (1j * (vo + 1j * alpha)) - B96_char_func( v, T, r, kappa_v, theta_v, sigma_v, rho, v0, lamb, mu, delta ) / ((vo + 1j * alpha) ** 2 - 1j * (vo + 1j * alpha)) ) # 数值FFT求解 delt = np.zeros(N) delt[0] = 1 j = np.arange(1, N + 1, 1) SimpsonW = (3 + (-1) ** j - delt) / 3 # 辛普森权重，提高积分精度 if S0 >= 0.95 * K: FFTFunc = np.exp(1j * b * vo) * modcharFunc * eta * SimpsonW payoff = (np.fft.fft(FFTFunc)).real CallValueM = np.exp(-alpha * k) / np.pi * payoff else: FFTFunc = ( np.exp(1j * b * vo) * (modcharFunc1 - modcharFunc2) * 0.5 * eta * SimpsonW ) payoff = (np.fft.fft(FFTFunc)).real CallValueM = payoff / (np.sinh(alpha * k) * np.pi) # 定位对应执行价的期权价格 pos = int((k + b) / eps) CallValue = CallValueM[pos] * S0 return CallValue
# 基于Lewis方法call_price_lewis = B96_call_value_Lewis(100, 100, 1, 0.05, 2.0, 0.04, 0.5, -0.5, 0.04, 1.0, -0.2, 0.1)# 基于FFT方法call_price_fft = B96_call_FFT(100, 100, 1, 0.05, 2.0, 0.04, 0.5, -0.5, 0.04, 1.0, -0.2, 0.1)print("Lewis方法定价：", call_price_lewis)print("FFT方法定价：", call_price_fft)
- 数据格式：
.npy是 NumPy 高效存储格式，适合期权定价等数值计算场景 - 核心难点：跳跃扩散模型的参数不稳定性，正则化是行业通用解决手段
- 校准优势：序贯校准（先随机波动率、再跳跃、最后全模型）能有效规避局部最优，提升参数稳定性
- 模型局限：Bates (1996) 假设无风险利率固定，未考虑随机利率
Heston、Merton、Bates 模型与 FFT 期权定价一、核心模型体系
1. 三大基础模型
| | | |
|---|
| Heston (1993) | | 刻画波动率的连续时变、聚类特征，捕捉资产收益率与波动率的负相关性 | |
| Merton (1976) | | 刻画资产价格的不连续突发跳变，拟合极端行情下的期权定价 | |
| Bates (1996) | | 同时捕捉连续波动率波动和不连续价格跳跃，兼顾两类风险 | |
二、核心数值方法：FFT 期权定价
1. FFT 核心作用
基于 Carr-Madan/Lewis 积分框架，将期权定价的解析解转化为快速傅里叶变换，大幅提升大规模期权定价的计算效率，是 Heston/Merton/Bates 模型的标准数值实现方式。
2. FFT 关键参数
| | |
|---|
| g（积分截断上限） | 控制数值积分的范围，保证被积函数充分收敛，提升定价精度 | g |
| α（阻尼因子） | | 核心决定条件：模型是否存在跳跃 |
| N（采样点数） | | 与eps（频域网格步长）强相关：eps=1/(g×N)，共同决定网格疏密 |
| 积分截断 | 理论积分区间为(−∞,+∞)，实际截断为有限值（如 50） | |
三、模型校准全流程
1. 校准方式
- 序贯校准（Sequential Calibration）：先校准 Heston 随机波动率参数，固定后再校准 Merton 跳跃参数，降低优化难度
- 完整校准（Full Calibration）：同时校准所有 Bates 模型参数，拟合效果最优
2. 校准核心参数与约束
| | |
|---|
| μ | −0.6<μ<0 | 符合股票市场 “下跌跳空更常见” 的实证规律，避免极端无意义值（如μ<−1，跌幅超 100%） |
| ρ | −1<ρ<0 | |
| tol | | 核心作用：筛选平值（ATM）期权，仅保留行权价在S0的±2%区间内的期权，保证校准数据的流动性与代表性 |
3. 拟合效果评估指标
- MSE（均方误差，Mean Squared Error）：核心评估指标，
纯 Merton（jump-only）校准：MSE 较高，拟合偏差大- Bates（full）完整校准：MSE 显著下降，拟合效果大幅提升
四、工程实现与工具
1. 数据存储：.npy文件格式
- 核心优势：Python/NumPy 读写速度极快，二进制原生存储，无需字符串解析，适配大规模金融数据处理
- 其他特性：文件体积更小，仅在 Python 生态内高效，不具备人类可读性
2. 校准质量可视化
- 标准方法：按到期期限（maturity）分组，绘制市场价格 vs 模型价格对比图，直观展示不同区间的拟合效果
- 辅助方法：误差直方图（统计误差分布）、隐含波动率曲面（展示整体拟合）
- Heston、Merton、Bates 是量化金融核心期权定价模型，FFT 是其关键数值实现方法
- 模型校准的核心是参数估计，需通过tol等参数筛选代表性数据，结合市场特征设定参数约束，保证估计的合理性与准确性
- 随机波动率、跳跃幅度、积分范围等参数直接决定模型拟合效果，MSE 是核心评估指标
- Bates 模型融合了 Heston 与 Merton 的优势，是当前市场拟合效果最优的经典期权定价模型之一
