1. 项目概述与COS方法核心思想在量化金融和衍生品定价的日常工作中我们经常需要面对一个核心问题如何快速且精确地为复杂模型下的期权定价。传统的蒙特卡洛模拟虽然通用但计算成本高昂而解析解往往只存在于少数理想化的模型中。当遇到像Heston随机波动率模型这类能更真实刻画市场“波动率微笑”现象但又没有简单闭式解的模型时一个高效的数值方法就显得至关重要。COS方法Fourier-Cosine Series Expansion Method正是为解决这类问题而生的利器。我第一次接触这个方法时就被其优雅的数学形式和惊人的计算效率所吸引——它巧妙地将概率密度函数的特征函数与傅里叶余弦级数展开相结合把定价问题转化为一系列余弦系数的加权求和从而绕开了直接求解复杂积分或偏微分方程的难题。简单来说COS方法的核心思想是任何概率密度函数只要我们知道其特征函数即其傅里叶变换就可以在一个有限的区间[a, b]内用余弦级数来逼近。而欧式期权的价格本质上就是收益函数关于标的资产到期日价格分布的期望贴现。因此我们可以先将收益函数也用同样的余弦基展开那么期权价格的积分就神奇地转化为两个余弦级数系数序列的点积计算量骤减。对于Heston模型其对数收益的特征函数有已知的半解析形式这使得COS方法成为了一个“天作之合”的选择——我们无需模拟成千上万条路径只需计算有限项级数就能得到高精度的价格。本文将手把手带你用R语言实现这一过程。无论你是刚刚踏入量化领域的新手还是希望优化现有定价库的老兵这篇实践指南都将从理论到代码完整呈现如何将一篇论文中的公式比如Junike (2024)关于截断区间和项数选择的最新成果转化为可运行、可调试、可复现的工业级代码。我们会深入每个函数背后的金融数学原理讨论关键参数如截断范围L和展开项数N选择的“艺术”与“科学”并分享我在实现过程中踩过的坑和总结的调试技巧。最终你将获得一套可以直接嵌入你策略回测或风险管理系统中的Heston模型COS定价模块。2. Heston模型与特征函数解析2.1 Heston模型参数的经济含义在实现COS方法之前我们必须先理解定价的“对象”——Heston模型。它是对经典Black-Scholes模型的重要扩展引入了随机波动率。模型由一对随机微分方程描述标的资产价格过程dS_t μ S_t dt √v_t S_t dW_t^S波动率过程dv_t κ (θ - v_t) dt ξ √v_t dW_t^v其中dW_t^S和dW_t^v是相关的布朗运动相关系数为ρ。模型包含五个核心参数每一个都有明确的金融含义κ (kappa) - 均值回归速度它衡量波动率向其长期平均水平回归的快慢。κ值越大波动率在受到冲击后回归到θ的速度越快。在实际市场中波动率通常表现出均值回归特性即不会长期停留在极高或极低的水平。θ (theta) - 长期平均方差波动率过程v_t在长期内趋向的平均水平。注意这里是方差所以长期平均波动率是√θ。ξ (xi) - 波动率的波动率这个参数描述了波动率本身的不确定性或“波动”。ξ越大波动率过程越“狂野”更容易产生极端的波动率值从而在期权市场上生成更显著的“波动率微笑”。ρ (rho) - 相关系数资产收益与波动率变化之间的相关性。在股票市场中我们通常观察到ρ 0即市场下跌时负收益波动率往往上升这被称为“杠杆效应”。这个参数对于生成不对称的微笑曲线至关重要。v0 - 初始方差在定价时刻观察到的或隐含的初始波动率水平。理解这些参数是后续校准和模型应用的基础。在代码中我们将用一个长度为5的向量params c(kappa, theta, xi, rho, v0)来统一管理它们。2.2 Heston模型特征函数的推导与实现COS方法的基石是特征函数。对于对数资产价格x_T log(S_T)其特征函数定义为ψ(u) E[e^(i u x_T)]。幸运的是对于Heston模型这个期望有半解析形式的解参见Heston (1993) 或 Schoutens et al. (2004)。虽然公式看起来复杂但我们可以将其分解为几个有明确含义的部分。下面是我们将在R中实现的特征函数psiLogST_Heston。我强烈建议你在阅读代码时对照下面的分解说明# 特征函数对数资产价格 ln(S_T) 的特征函数 # u: 傅里叶变换变量 mat: 到期时间 params: 模型参数向量 S0: 标的现价 r: 无风险利率 psiLogST_Heston function(u, mat, params, S0, r){ kappa params[1] # 均值回归速度 theta params[2] # 长期平均方差 xi params[3] # 波动率的波动率 rho params[4] # 相关系数 v0 params[5] # 初始方差 # 1. 计算辅助变量 d它与复平面上的特征根有关 d sqrt((rho * xi * u * 1i - kappa)^2 - xi^2 * (-1i * u - u^2)) # 2. 计算中间变量 mytmp 和 g mytmp kappa - rho * xi * u * 1i g (mytmp - d) / (mytmp d) # 与边界条件相关的比率 # 3. 计算时间衰减项 expdmat exp(-d * mat) # 4. 组合特征函数的三个部分 # tmp0: 与初始价格和漂移项相关的部分 tmp0 1i * u * (log(S0) r * mat) # tmp1 tmp2: 与波动率过程均值回归部分相关的贡献 tmp1 (mytmp - d) * mat - 2 * log((1 - g * expdmat) / (1 - g)) tmp2 theta * kappa * xi^(-2) * tmp1 # tmp3: 与初始波动率 v0 相关的贡献 tmp3 v0 * xi^(-2) * (mytmp - d) * (1 - expdmat) / (1 - g * expdmat) # 5. 返回特征函数值 exp(tmp0 tmp2 tmp3) }注意代码中的1i是R语言中虚数单位i的表示。整个计算都在复数域中进行。虽然公式复杂但现代计算机处理复数运算速度很快我们无需手动分离实部和虚部。为什么需要中心化直接使用psiLogST_Heston得到的是log(S_T)的特征函数。在COS方法中为了优化级数展开的收敛性我们通常会对随机变量进行中心化处理即考虑y log(S_T) - E[log(S_T)]。中心化后的变量y均值为0其概率质量更集中所需的截断区间[a, b]可以更小从而用更少的展开项N达到相同的精度。因此我们需要计算中心化对数收益的特征函数φ(u) E[e^(i u y)] ψ(u) * e^(-i u * μ)其中μ E[log(S_T)]。μ可以通过特征函数在原点处的导数求得μ E[log(S_T)] -i * ψ(0)。在R中我们可以用Deriv包进行符号求导但这里有个重要的实操心得对于生产环境频繁调用符号求导函数Deriv来计算μ效率很低因为每次定价都要重新求导。更好的做法是手动推导出μ的解析表达式并编码。根据Heston模型的性质可以推导出μ log(S0) r * mat - 0.5 * theta * mat 0.5 * (theta - v0) * (1 - exp(-kappa * mat)) / kappa注意这个表达式依赖于具体的特征函数形式可能与上述代码对应的特征函数略有不同需要根据你的psi函数形式进行验证或重新推导。在初步实现和验证阶段我们可以先用Deriv但在性能优化时替换为解析解是必须的步骤。3. COS方法的核心算法拆解3.1 从连续积分到离散求和数学变换COS方法的魔力在于它将一个复杂的积分问题简化了。欧式看跌期权的价格看涨期可通过看跌-看涨平价关系得到为P e^(-rT) * E[ max(K - S_T, 0) ] e^(-rT) ∫_a^b max(K - e^x, 0) f(x) dx其中f(x)是xlog(S_T)的概率密度函数[a, b]是积分截断区间。COS方法的核心步骤如下密度函数展开在区间[a, b]上用余弦级数逼近密度函数f(x)f(x) ≈ Σ_{k0}^{N-1} A_k * cos(kπ (x-a)/(b-a))其中Σ表示求和第一项权重为1/2A_k是展开系数。系数关联特征函数关键的一步是这些余弦系数A_k可以通过密度函数的特征函数φ(u)即中心化后的特征函数直接计算而无需知道f(x)的具体形式A_k (2/(b-a)) * Re{ φ(kπ/(b-a)) * exp(-i kπ a/(b-a)) }收益函数展开同样将看跌期权的收益函数g(x)max(K-e^x, 0)在同一个区间上用相同的余弦基展开得到系数V_k。价格近似将两个展开式代入积分利用余弦函数的正交性积分奇迹般地简化为系数点积P ≈ e^(-rT) * Σ_{k0}^{N-1} A_k * V_k至此问题转化为计算两个序列{A_k}和{V_k}并求它们的加权和。计算A_k需要特征函数φ计算V_k有解析公式见下文。3.2 关键参数截断区间[a, b]与项数N精度和效率的权衡就体现在这两个参数上。区间[a, b]选得太窄会截断密度函数的尾部导致误差选得太宽则需要更大的N来捕捉函数细节降低效率。项数N直接决定了计算量N越大精度越高但速度越慢。早期COS方法论文如Fang Oosterlee, 2009常用基于随机变量前几阶矩的经验公式来确定[a, b]例如a c1 - L√c2,b c1 L√c2其中c1, c2是中心化变量的前两阶矩L是一个常数如10。这种方法简单但不够严谨L的选择比较随意。本文代码采用了Junike (2024) 提出的更严谨的方法。它基于对截断误差和级数截断误差的严格数学界反向推导出给定误差容忍度eps下最优的L这里区间定义为[-L, L]因为变量已中心化和N。计算L区间半宽L (2 * K * e^(-rT) * μ_n / eps)^(1/4)其中μ_n是中心化对数收益的四阶绝对矩E[|y|^4]。这个公式的直觉是收益函数的尾部行为由矩控制和允许的误差共同决定了需要覆盖的区间范围。四阶矩可以通过特征函数的四阶导数在0点的值计算μ_n |φ^{(4)}(0)|。计算N展开项数 公式更为复杂见代码它涉及对特征函数高阶导数衰减性的积分估计boundDeriv。简单理解特征函数φ(u)在u很大时衰减得越快即分布越平滑所需的N就越小。这个计算需要数值积分是预处理中计算成本最高的部分但好消息是对于固定的参数集(params, mat, S0, r)它只需要计算一次。我们可以将结果保存下来如save(phi4, file“phi4.RData”)后续定价直接调用。实操心得在实际应用中如果需要对同一标的、同一到期日但不同行权价K的多个期权定价如计算整个波动率曲面利用看跌-看涨平价关系我们只需要计算一次看跌期权价格。但注意上述L的公式依赖于K。一个常见的简化是对所有K使用一个统一的、足够大的L例如基于最大或最小行权价计算虽然可能略微降低单个期权的计算效率但避免了为每个K重复计算L和N在批量定价时反而可能更快。这需要根据具体应用场景进行权衡。4. R语言实现全流程与代码精讲4.1 环境准备与辅助函数实现首先确保安装了必要的R包。我们主要依赖Deriv包进行符号求导来计算高阶矩。对于生产环境建议将求导结果硬编码为解析函数以提升速度。# 加载符号求导包用于初始开发和验证 if (!requireNamespace(Deriv, quietly TRUE)) install.packages(Deriv) library(Deriv) # 1. 定义Heston模型对数价格的特征函数 (已在前文给出) psiLogST_Heston - function(u, mat, params, S0, r) { # ... 函数体同上此处省略 ... } # 2. 计算 E[log(S_T)] μ # 方法1使用符号求导方便但慢用于验证 psi_prime_at_0 - Deriv(psiLogST_Heston, u) mu_via_deriv - function(mat, params, S0, r) { Re(-1i * psi_prime_at_0(0, mat, params, S0, r)) } # 方法2解析解推荐用于生产 # 注意此解析解需要与你使用的psiLogST_Heston函数形式严格匹配。 # 下面是一个常见形式的Heston特征函数对应的μ解析解示例请务必根据你的公式推导验证。 mu_analytic - function(mat, params, S0, r) { kappa - params[1]; theta - params[2]; xi - params[3]; rho - params[4]; v0 - params[5] # 假设psi函数形式导致 E[log(S_T)] log(S0) r*mat - 0.5 * E[∫0^T v_t dt] # E[∫0^T v_t dt] 在Heston模型下有解析解 expected_integrated_var - theta * mat (v0 - theta) * (1 - exp(-kappa * mat)) / kappa return(log(S0) r * mat - 0.5 * expected_integrated_var) } # 在后续代码中我们将使用mu_analytic并假设其正确。实际使用时应用psi_prime_at_0在测试用例上验证mu_analytic的正确性。 # 3. 中心化对数收益的特征函数 phi - function(u, mat, params, S0, r) { # 使用解析解μ以提高速度 mu_val - mu_analytic(mat, params, S0, r) return(psiLogST_Heston(u, mat, params, S0, r) * exp(-1i * u * mu_val)) }4.2 余弦系数A_k与收益系数V_k的计算这是COS方法的核心计算步骤。# 4. 计算密度函数的余弦展开系数 A_k (代码中的ck函数) # L: 中心化区间半宽 [-L, L], N: 展开项数计算0,...,N ck - function(L, mat, N, params, S0, r) { k - 0:N # 注意公式中的 a -L, b L, 所以 (b-a)2L # A_k (2/(2L)) * Re { φ(kπ/(2L)) * exp(-i kπ * a/(2L)) } (1/L) * Re { φ(kπ/(2L)) * exp(i kπ/2) } # 因为 a -L, 所以 -i * kπ * a / (2L) i * kπ/2 u_vals - k * pi / (2 * L) return( (1/L) * Re( phi(u_vals, mat, params, S0, r) * exp(1i * k * pi / 2) ) ) } # 5. 计算看跌期权收益函数的余弦展开系数 V_k (代码中的vk函数) # 收益函数g(x) max(K - exp(x), 0), 其中 x log(S_T) # 经过推导参见Junike and Pankrashkin (2022)附录V_k有闭式解。 vk_put - function(K, L, mat, N, params, S0, r) { mu_val - mu_analytic(mat, params, S0, r) # 中心化前的均值 d - min(log(K) - mu_val, L) # 收益函数非零区间的右端点因为看跌期权在价格低于K时才有收益 if (d -L) { # 如果整个区间[-L, L]上收益函数都为0即K极小则所有系数为0 return(rep(0, N 1)) } k - 0:N # 计算两个辅助函数 Ψ0 和 Ψ1 # Ψ0 对应收益函数中常数K部分的积分 psi0 - (2 * L / (k * pi)) * sin(k * pi * (d L) / (2 * L)) psi0[1] - d L # k0时的特例因为sin(0)/0需要单独处理 # Ψ1 对应收益函数中 -exp(x) 部分的积分 tmp1 - (k * pi / (2 * L)) * sin(k * pi * (d L) / (2 * L)) tmp2 - cos(k * pi * (d L) / (2 * L)) tmp3 - 1 (k * pi / (2 * L))^2 psi1 - (exp(d) * (tmp1 tmp2) - exp(-L)) / tmp3 # 组合成V_k系数并贴现 return(exp(-r * mat) * (K * psi0 - exp(mu_val) * psi1)) }4.3 定价函数与看跌-看涨平价有了系数定价就是点积求和。# 6. COS方法近似计算欧式看跌期权价格 put_COS - function(K, L, mat, N, params, S0, r) { A_k - ck(L, mat, N, params, S0, r) V_k - vk_put(K, L, mat, N, params, S0, r) # 求和第一项k0权重为1/2 sum_val - sum(A_k * V_k) sum_val - sum_val - 0.5 * A_k[1] * V_k[1] # 先全加再减去一半的第一项 # 更清晰的写法sum_val 0.5*A_0*V_0 Σ_{k1}^{N} A_k * V_k # 但上面利用向量化运算更高效 return(sum_val) } # 7. 利用看跌-看涨平价关系计算欧式看涨期权价格 # 看跌-看涨平价C K*e^{-rT} P S0 call_COS - function(K, L, mat, N, params, S0, r) { put_price - put_COS(K, L, mat, N, params, S0, r) return(put_price S0 - K * exp(-r * mat)) }4.4 自动确定L与N的完整流程现在我们将Junike (2024)的方法集成进来实现给定误差容忍度下的参数自动选择。# 8. 高阶导数计算预处理只需算一次 # 计算中心化特征函数的1-4阶导数用于计算矩和误差界 phi1 - Deriv(phi, u) phi2 - Deriv(phi1, u) phi3 - Deriv(phi2, u) phi4 - Deriv(phi3, u) # 计算四阶导数用于求四阶矩 # 9. 根据目标精度自动计算L和N的函数 calculate_L_N - function(K, mat, params, S0, r, eps 1e-6, s 20) { # 计算四阶绝对矩 μ_n |φ^{(4)}(0)| mu_n - abs(phi4(0, mat, params, S0, r)) # 计算截断区间半宽 L公式来自 Junike (2024, Eq. (3.10)) L - (2 * K * exp(-r * mat) * mu_n / eps)^(1/4) # 计算误差界积分用于确定N integrand - function(u) { (1 / (2 * pi)) * (abs(u)^(s 1)) * abs(phi(u, mat, params, S0, r)) } # 数值积分注意积分区间为整个实数轴。实践中可截断到足够大的范围。 # 使用 integrate 函数设置适当的下限和上限。 boundDeriv - integrate(integrand, lower -100, upper 100, subdivisions 1000)$value # 注意积分限应足够大使得phi(u)衰减到接近0。对于Heston模型|u|很大时|phi(u)|衰减很快。 # 计算所需项数 N公式来自 Junike (2024, Sec. 6.1) numerator - 2^(s 5/2) * boundDeriv * L^(s 2) * 12 * K * exp(-r * mat) denominator - s * pi^(s 1) * eps N - ceiling((numerator / denominator)^(1 / s)) return(list(L L, N N, mu_n mu_n)) } # 10. 封装一个用户友好的定价函数 price_option_COS - function(K, mat, params, S0, r, type call, eps 1e-6) { # 步骤1计算该参数集下的最优L和N LN_info - calculate_L_N(K, mat, params, S0, r, eps) L - LN_info$L N - LN_info$N # 步骤2计算价格 if (type put) { price - put_COS(K, L, mat, N, params, S0, r) } else if (type call) { price - call_COS(K, L, mat, N, params, S0, r) } else { stop(Option type must be either call or put.) } # 返回价格和使用的参数用于调试和验证 return(list(price price, L L, N N, mu_n LN_info$mu_n)) }5. 实战测试、常见问题与性能优化5.1 完整测试用例与结果验证让我们用文献中常见的参数来测试我们的实现。这里我们使用与输入材料中相同的参数。# 设置测试参数 S0 - 100 # 标的现价 K - 90 # 行权价 r - 0.1 # 无风险利率 mat - 0.7 # 到期时间年 params - c( kappa 0.6067, # 均值回归速度 theta 0.0707, # 长期平均方差 xi 0.2928, # 波动率的波动率 rho -0.7571, # 相关系数 v0 0.0654 # 初始方差 ) # 设置误差容忍度 eps - 1e-6 # 计算看跌期权价格 result_put - price_option_COS(K, mat, params, S0, r, type put, eps eps) cat(sprintf(欧式看跌期权价格: %.6f\n, result_put$price)) cat(sprintf(使用的参数: L %.4f, N %d\n, result_put$L, result_put$N)) # 利用看跌-看涨平价计算看涨期权价格 call_price_by_parity - result_put$price S0 - K * exp(-r * mat) # 直接用COS方法计算看涨期权价格 result_call - price_option_COS(K, mat, params, S0, r, type call, eps eps) cat(sprintf(欧式看涨期权价格 (COS): %.6f\n, result_call$price)) cat(sprintf(欧式看涨期权价格 (平价): %.6f\n, call_price_by_parity)) cat(sprintf(两者差异: %.10f\n, abs(result_call$price - call_price_by_parity)))运行这段代码我们应该得到看跌期权价格约为2.773954这与输入材料中给出的结果一致。看涨期权的COS计算结果与通过看跌-看涨平价计算的结果应该只有极微小的数值误差远小于eps这验证了我们代码实现的一致性。5.2 常见问题排查与调试技巧在实现和运行过程中你可能会遇到以下问题结果为NaN或Inf检查特征函数首先验证psiLogST_Heston函数。对于较大的虚数参数u复数运算可能导致溢出。可以尝试在计算exp(d * mat)等项时检查中间变量d的实部是否为负确保衰减。Heston模型的特征函数在参数满足Feller条件2κθ ξ^2时通常表现良好。我们的测试参数满足该条件。检查积分区间在calculate_L_N函数中数值积分integrate的上下限(-100, 100)可能不够大。可以尝试增加上限如200, 500并观察积分值是否收敛。也可以先画出integrand(u)的函数图像看看它在多大范围内有显著非零值。计算速度慢瓶颈分析最耗时的部分是calculate_L_N中的数值积分和phi4的高阶符号求导。对于固定(params, mat, S0, r)的批量定价务必预计算并保存phi4函数和boundDeriv积分结果。向量化我们的ck和vk_put函数已经对k0:N进行了向量化运算这比用循环快得多。简化μ的计算如前所述用解析解mu_analytic替代通过Deriv得到的mu_via_deriv。精度不达标增大N如果结果与参考值偏差较大可以手动设置一个更大的N比如128, 256进行测试看看价格是否收敛。如果收敛到一个不同的值可能是L计算有误。检查L手动增大L比如乘以1.5或2如果价格发生显著变化说明初始L可能截断了密度函数的重要尾部。需要检查四阶矩mu_n的计算是否正确。验证特征函数与已知的Heston特征函数实现进行交叉验证例如与QuantLib或其他可靠库比较。可以计算psiLogST_Heston在一些u值上的输出。看跌-看涨平价不成立如果COS计算的看涨和看跌价格不满足平价关系问题几乎肯定出在vk_put函数或ck函数上。首先检查mu_analytic的计算是否正确与mu_via_deriv对比。然后用简单的数值积分方法如对密度函数采样计算一个基准价格来孤立问题。5.3 性能优化与生产环境建议对于需要高频或批量定价的生产环境可以考虑以下优化预计算与缓存这是最重要的优化。创建一个缓存系统为每个不同的参数组合(params, mat, S0, r)预计算并存储以下信息中心化均值mu_val四阶矩mu_n误差界积分值boundDeriv甚至可以为常用的(L, N)组合预计算好ck系数因为它不依赖于K。替换符号求导完全避免使用Deriv包。手动推导出phi(u)的1-4阶导数公式并实现为R函数。这需要一些复杂的微积分工作但能带来巨大的速度提升。使用更快的数值积分器integrate函数可能不是最快的。可以尝试使用更专业的数值积分库或者对于衰减很快的函数使用自适应高斯-埃尔米特求积法。并行计算如果需要对大量不同行权价K的期权定价由于每个K的计算是独立的可以很容易地用parallel包进行并行计算。代码移植对于极致性能要求可以考虑用C重写核心计算部分通过Rcpp包集成到R中特别是特征函数和系数计算的循环部分5.4 扩展与应用隐含波动率计算与模型校准一旦我们有了快速准确的定价器两个最直接的应用就是计算隐含波动率和模型校准。计算隐含波动率给定市场观察到的期权价格C_market我们需要找到使得Heston模型价格C_model(σ)等于市场价格的波动率参数这里其他Heston参数固定。由于COS定价器很快我们可以使用标准的求根算法如uniroot来反解。# 示例计算隐含波动率 (这里假设只变动v0其他参数固定) # 注意Heston模型的“隐含波动率”通常指与BS模型波动率σ的映射这里我们反解的是v0。 # 更常见的校准是反解所有参数。 calc_implied_v0 - function(market_price, K, mat, params_guess, S0, r, typecall) { # params_guess是包含v0初始猜测的参数向量 price_diff - function(v0_trial) { params_trial - params_guess params_trial[5] - v0_trial # 替换v0 model_price - price_option_COS(K, mat, params_trial, S0, r, typetype, eps1e-6)$price return(model_price - market_price) } # 使用二分法求根 root_result - uniroot(price_diff, interval c(0.001, 1.0)) # 给v0一个合理的搜索区间 return(root_result$root) }模型校准这是量化金融中的核心任务。给定一组不同行权价和到期日的期权市场价格我们寻找一组Heston模型参数(κ, θ, ξ, ρ, v0)使得模型价格与市场价格的总体误差如均方误差最小。这通常是一个高维非凸优化问题可以使用optim或nlminb等R优化函数结合我们的COS定价器来完成。# 简化的校准框架示例 calibrate_heston - function(market_data, params_init, S0, r) { # market_data: data.frame with columns: K, mat, market_price, type objective_function - function(params) { total_error - 0 for (i in 1:nrow(market_data)) { row - market_data[i, ] model_price - price_option_COS(row$K, row$mat, params, S0, r, typerow$type, eps1e-6)$price total_error - total_error (model_price - row$market_price)^2 } return(total_error) } # 添加参数约束如κ, θ, ξ 0, |ρ| 1, v0 0 lower_bounds - c(1e-3, 1e-3, 1e-3, -0.999, 1e-3) upper_bounds - c(10, 1, 5, 0.999, 1) result - optim(params_init, objective_function, method L-BFGS-B, lower lower_bounds, upper upper_bounds) return(result$par) }在校准过程中定价速度至关重要。因此之前提到的预计算和性能优化在这里会得到丰厚的回报。
R语言实现Heston模型COS期权定价：从傅里叶变换到高效数值计算
1. 项目概述与COS方法核心思想在量化金融和衍生品定价的日常工作中我们经常需要面对一个核心问题如何快速且精确地为复杂模型下的期权定价。传统的蒙特卡洛模拟虽然通用但计算成本高昂而解析解往往只存在于少数理想化的模型中。当遇到像Heston随机波动率模型这类能更真实刻画市场“波动率微笑”现象但又没有简单闭式解的模型时一个高效的数值方法就显得至关重要。COS方法Fourier-Cosine Series Expansion Method正是为解决这类问题而生的利器。我第一次接触这个方法时就被其优雅的数学形式和惊人的计算效率所吸引——它巧妙地将概率密度函数的特征函数与傅里叶余弦级数展开相结合把定价问题转化为一系列余弦系数的加权求和从而绕开了直接求解复杂积分或偏微分方程的难题。简单来说COS方法的核心思想是任何概率密度函数只要我们知道其特征函数即其傅里叶变换就可以在一个有限的区间[a, b]内用余弦级数来逼近。而欧式期权的价格本质上就是收益函数关于标的资产到期日价格分布的期望贴现。因此我们可以先将收益函数也用同样的余弦基展开那么期权价格的积分就神奇地转化为两个余弦级数系数序列的点积计算量骤减。对于Heston模型其对数收益的特征函数有已知的半解析形式这使得COS方法成为了一个“天作之合”的选择——我们无需模拟成千上万条路径只需计算有限项级数就能得到高精度的价格。本文将手把手带你用R语言实现这一过程。无论你是刚刚踏入量化领域的新手还是希望优化现有定价库的老兵这篇实践指南都将从理论到代码完整呈现如何将一篇论文中的公式比如Junike (2024)关于截断区间和项数选择的最新成果转化为可运行、可调试、可复现的工业级代码。我们会深入每个函数背后的金融数学原理讨论关键参数如截断范围L和展开项数N选择的“艺术”与“科学”并分享我在实现过程中踩过的坑和总结的调试技巧。最终你将获得一套可以直接嵌入你策略回测或风险管理系统中的Heston模型COS定价模块。2. Heston模型与特征函数解析2.1 Heston模型参数的经济含义在实现COS方法之前我们必须先理解定价的“对象”——Heston模型。它是对经典Black-Scholes模型的重要扩展引入了随机波动率。模型由一对随机微分方程描述标的资产价格过程dS_t μ S_t dt √v_t S_t dW_t^S波动率过程dv_t κ (θ - v_t) dt ξ √v_t dW_t^v其中dW_t^S和dW_t^v是相关的布朗运动相关系数为ρ。模型包含五个核心参数每一个都有明确的金融含义κ (kappa) - 均值回归速度它衡量波动率向其长期平均水平回归的快慢。κ值越大波动率在受到冲击后回归到θ的速度越快。在实际市场中波动率通常表现出均值回归特性即不会长期停留在极高或极低的水平。θ (theta) - 长期平均方差波动率过程v_t在长期内趋向的平均水平。注意这里是方差所以长期平均波动率是√θ。ξ (xi) - 波动率的波动率这个参数描述了波动率本身的不确定性或“波动”。ξ越大波动率过程越“狂野”更容易产生极端的波动率值从而在期权市场上生成更显著的“波动率微笑”。ρ (rho) - 相关系数资产收益与波动率变化之间的相关性。在股票市场中我们通常观察到ρ 0即市场下跌时负收益波动率往往上升这被称为“杠杆效应”。这个参数对于生成不对称的微笑曲线至关重要。v0 - 初始方差在定价时刻观察到的或隐含的初始波动率水平。理解这些参数是后续校准和模型应用的基础。在代码中我们将用一个长度为5的向量params c(kappa, theta, xi, rho, v0)来统一管理它们。2.2 Heston模型特征函数的推导与实现COS方法的基石是特征函数。对于对数资产价格x_T log(S_T)其特征函数定义为ψ(u) E[e^(i u x_T)]。幸运的是对于Heston模型这个期望有半解析形式的解参见Heston (1993) 或 Schoutens et al. (2004)。虽然公式看起来复杂但我们可以将其分解为几个有明确含义的部分。下面是我们将在R中实现的特征函数psiLogST_Heston。我强烈建议你在阅读代码时对照下面的分解说明# 特征函数对数资产价格 ln(S_T) 的特征函数 # u: 傅里叶变换变量 mat: 到期时间 params: 模型参数向量 S0: 标的现价 r: 无风险利率 psiLogST_Heston function(u, mat, params, S0, r){ kappa params[1] # 均值回归速度 theta params[2] # 长期平均方差 xi params[3] # 波动率的波动率 rho params[4] # 相关系数 v0 params[5] # 初始方差 # 1. 计算辅助变量 d它与复平面上的特征根有关 d sqrt((rho * xi * u * 1i - kappa)^2 - xi^2 * (-1i * u - u^2)) # 2. 计算中间变量 mytmp 和 g mytmp kappa - rho * xi * u * 1i g (mytmp - d) / (mytmp d) # 与边界条件相关的比率 # 3. 计算时间衰减项 expdmat exp(-d * mat) # 4. 组合特征函数的三个部分 # tmp0: 与初始价格和漂移项相关的部分 tmp0 1i * u * (log(S0) r * mat) # tmp1 tmp2: 与波动率过程均值回归部分相关的贡献 tmp1 (mytmp - d) * mat - 2 * log((1 - g * expdmat) / (1 - g)) tmp2 theta * kappa * xi^(-2) * tmp1 # tmp3: 与初始波动率 v0 相关的贡献 tmp3 v0 * xi^(-2) * (mytmp - d) * (1 - expdmat) / (1 - g * expdmat) # 5. 返回特征函数值 exp(tmp0 tmp2 tmp3) }注意代码中的1i是R语言中虚数单位i的表示。整个计算都在复数域中进行。虽然公式复杂但现代计算机处理复数运算速度很快我们无需手动分离实部和虚部。为什么需要中心化直接使用psiLogST_Heston得到的是log(S_T)的特征函数。在COS方法中为了优化级数展开的收敛性我们通常会对随机变量进行中心化处理即考虑y log(S_T) - E[log(S_T)]。中心化后的变量y均值为0其概率质量更集中所需的截断区间[a, b]可以更小从而用更少的展开项N达到相同的精度。因此我们需要计算中心化对数收益的特征函数φ(u) E[e^(i u y)] ψ(u) * e^(-i u * μ)其中μ E[log(S_T)]。μ可以通过特征函数在原点处的导数求得μ E[log(S_T)] -i * ψ(0)。在R中我们可以用Deriv包进行符号求导但这里有个重要的实操心得对于生产环境频繁调用符号求导函数Deriv来计算μ效率很低因为每次定价都要重新求导。更好的做法是手动推导出μ的解析表达式并编码。根据Heston模型的性质可以推导出μ log(S0) r * mat - 0.5 * theta * mat 0.5 * (theta - v0) * (1 - exp(-kappa * mat)) / kappa注意这个表达式依赖于具体的特征函数形式可能与上述代码对应的特征函数略有不同需要根据你的psi函数形式进行验证或重新推导。在初步实现和验证阶段我们可以先用Deriv但在性能优化时替换为解析解是必须的步骤。3. COS方法的核心算法拆解3.1 从连续积分到离散求和数学变换COS方法的魔力在于它将一个复杂的积分问题简化了。欧式看跌期权的价格看涨期可通过看跌-看涨平价关系得到为P e^(-rT) * E[ max(K - S_T, 0) ] e^(-rT) ∫_a^b max(K - e^x, 0) f(x) dx其中f(x)是xlog(S_T)的概率密度函数[a, b]是积分截断区间。COS方法的核心步骤如下密度函数展开在区间[a, b]上用余弦级数逼近密度函数f(x)f(x) ≈ Σ_{k0}^{N-1} A_k * cos(kπ (x-a)/(b-a))其中Σ表示求和第一项权重为1/2A_k是展开系数。系数关联特征函数关键的一步是这些余弦系数A_k可以通过密度函数的特征函数φ(u)即中心化后的特征函数直接计算而无需知道f(x)的具体形式A_k (2/(b-a)) * Re{ φ(kπ/(b-a)) * exp(-i kπ a/(b-a)) }收益函数展开同样将看跌期权的收益函数g(x)max(K-e^x, 0)在同一个区间上用相同的余弦基展开得到系数V_k。价格近似将两个展开式代入积分利用余弦函数的正交性积分奇迹般地简化为系数点积P ≈ e^(-rT) * Σ_{k0}^{N-1} A_k * V_k至此问题转化为计算两个序列{A_k}和{V_k}并求它们的加权和。计算A_k需要特征函数φ计算V_k有解析公式见下文。3.2 关键参数截断区间[a, b]与项数N精度和效率的权衡就体现在这两个参数上。区间[a, b]选得太窄会截断密度函数的尾部导致误差选得太宽则需要更大的N来捕捉函数细节降低效率。项数N直接决定了计算量N越大精度越高但速度越慢。早期COS方法论文如Fang Oosterlee, 2009常用基于随机变量前几阶矩的经验公式来确定[a, b]例如a c1 - L√c2,b c1 L√c2其中c1, c2是中心化变量的前两阶矩L是一个常数如10。这种方法简单但不够严谨L的选择比较随意。本文代码采用了Junike (2024) 提出的更严谨的方法。它基于对截断误差和级数截断误差的严格数学界反向推导出给定误差容忍度eps下最优的L这里区间定义为[-L, L]因为变量已中心化和N。计算L区间半宽L (2 * K * e^(-rT) * μ_n / eps)^(1/4)其中μ_n是中心化对数收益的四阶绝对矩E[|y|^4]。这个公式的直觉是收益函数的尾部行为由矩控制和允许的误差共同决定了需要覆盖的区间范围。四阶矩可以通过特征函数的四阶导数在0点的值计算μ_n |φ^{(4)}(0)|。计算N展开项数 公式更为复杂见代码它涉及对特征函数高阶导数衰减性的积分估计boundDeriv。简单理解特征函数φ(u)在u很大时衰减得越快即分布越平滑所需的N就越小。这个计算需要数值积分是预处理中计算成本最高的部分但好消息是对于固定的参数集(params, mat, S0, r)它只需要计算一次。我们可以将结果保存下来如save(phi4, file“phi4.RData”)后续定价直接调用。实操心得在实际应用中如果需要对同一标的、同一到期日但不同行权价K的多个期权定价如计算整个波动率曲面利用看跌-看涨平价关系我们只需要计算一次看跌期权价格。但注意上述L的公式依赖于K。一个常见的简化是对所有K使用一个统一的、足够大的L例如基于最大或最小行权价计算虽然可能略微降低单个期权的计算效率但避免了为每个K重复计算L和N在批量定价时反而可能更快。这需要根据具体应用场景进行权衡。4. R语言实现全流程与代码精讲4.1 环境准备与辅助函数实现首先确保安装了必要的R包。我们主要依赖Deriv包进行符号求导来计算高阶矩。对于生产环境建议将求导结果硬编码为解析函数以提升速度。# 加载符号求导包用于初始开发和验证 if (!requireNamespace(Deriv, quietly TRUE)) install.packages(Deriv) library(Deriv) # 1. 定义Heston模型对数价格的特征函数 (已在前文给出) psiLogST_Heston - function(u, mat, params, S0, r) { # ... 函数体同上此处省略 ... } # 2. 计算 E[log(S_T)] μ # 方法1使用符号求导方便但慢用于验证 psi_prime_at_0 - Deriv(psiLogST_Heston, u) mu_via_deriv - function(mat, params, S0, r) { Re(-1i * psi_prime_at_0(0, mat, params, S0, r)) } # 方法2解析解推荐用于生产 # 注意此解析解需要与你使用的psiLogST_Heston函数形式严格匹配。 # 下面是一个常见形式的Heston特征函数对应的μ解析解示例请务必根据你的公式推导验证。 mu_analytic - function(mat, params, S0, r) { kappa - params[1]; theta - params[2]; xi - params[3]; rho - params[4]; v0 - params[5] # 假设psi函数形式导致 E[log(S_T)] log(S0) r*mat - 0.5 * E[∫0^T v_t dt] # E[∫0^T v_t dt] 在Heston模型下有解析解 expected_integrated_var - theta * mat (v0 - theta) * (1 - exp(-kappa * mat)) / kappa return(log(S0) r * mat - 0.5 * expected_integrated_var) } # 在后续代码中我们将使用mu_analytic并假设其正确。实际使用时应用psi_prime_at_0在测试用例上验证mu_analytic的正确性。 # 3. 中心化对数收益的特征函数 phi - function(u, mat, params, S0, r) { # 使用解析解μ以提高速度 mu_val - mu_analytic(mat, params, S0, r) return(psiLogST_Heston(u, mat, params, S0, r) * exp(-1i * u * mu_val)) }4.2 余弦系数A_k与收益系数V_k的计算这是COS方法的核心计算步骤。# 4. 计算密度函数的余弦展开系数 A_k (代码中的ck函数) # L: 中心化区间半宽 [-L, L], N: 展开项数计算0,...,N ck - function(L, mat, N, params, S0, r) { k - 0:N # 注意公式中的 a -L, b L, 所以 (b-a)2L # A_k (2/(2L)) * Re { φ(kπ/(2L)) * exp(-i kπ * a/(2L)) } (1/L) * Re { φ(kπ/(2L)) * exp(i kπ/2) } # 因为 a -L, 所以 -i * kπ * a / (2L) i * kπ/2 u_vals - k * pi / (2 * L) return( (1/L) * Re( phi(u_vals, mat, params, S0, r) * exp(1i * k * pi / 2) ) ) } # 5. 计算看跌期权收益函数的余弦展开系数 V_k (代码中的vk函数) # 收益函数g(x) max(K - exp(x), 0), 其中 x log(S_T) # 经过推导参见Junike and Pankrashkin (2022)附录V_k有闭式解。 vk_put - function(K, L, mat, N, params, S0, r) { mu_val - mu_analytic(mat, params, S0, r) # 中心化前的均值 d - min(log(K) - mu_val, L) # 收益函数非零区间的右端点因为看跌期权在价格低于K时才有收益 if (d -L) { # 如果整个区间[-L, L]上收益函数都为0即K极小则所有系数为0 return(rep(0, N 1)) } k - 0:N # 计算两个辅助函数 Ψ0 和 Ψ1 # Ψ0 对应收益函数中常数K部分的积分 psi0 - (2 * L / (k * pi)) * sin(k * pi * (d L) / (2 * L)) psi0[1] - d L # k0时的特例因为sin(0)/0需要单独处理 # Ψ1 对应收益函数中 -exp(x) 部分的积分 tmp1 - (k * pi / (2 * L)) * sin(k * pi * (d L) / (2 * L)) tmp2 - cos(k * pi * (d L) / (2 * L)) tmp3 - 1 (k * pi / (2 * L))^2 psi1 - (exp(d) * (tmp1 tmp2) - exp(-L)) / tmp3 # 组合成V_k系数并贴现 return(exp(-r * mat) * (K * psi0 - exp(mu_val) * psi1)) }4.3 定价函数与看跌-看涨平价有了系数定价就是点积求和。# 6. COS方法近似计算欧式看跌期权价格 put_COS - function(K, L, mat, N, params, S0, r) { A_k - ck(L, mat, N, params, S0, r) V_k - vk_put(K, L, mat, N, params, S0, r) # 求和第一项k0权重为1/2 sum_val - sum(A_k * V_k) sum_val - sum_val - 0.5 * A_k[1] * V_k[1] # 先全加再减去一半的第一项 # 更清晰的写法sum_val 0.5*A_0*V_0 Σ_{k1}^{N} A_k * V_k # 但上面利用向量化运算更高效 return(sum_val) } # 7. 利用看跌-看涨平价关系计算欧式看涨期权价格 # 看跌-看涨平价C K*e^{-rT} P S0 call_COS - function(K, L, mat, N, params, S0, r) { put_price - put_COS(K, L, mat, N, params, S0, r) return(put_price S0 - K * exp(-r * mat)) }4.4 自动确定L与N的完整流程现在我们将Junike (2024)的方法集成进来实现给定误差容忍度下的参数自动选择。# 8. 高阶导数计算预处理只需算一次 # 计算中心化特征函数的1-4阶导数用于计算矩和误差界 phi1 - Deriv(phi, u) phi2 - Deriv(phi1, u) phi3 - Deriv(phi2, u) phi4 - Deriv(phi3, u) # 计算四阶导数用于求四阶矩 # 9. 根据目标精度自动计算L和N的函数 calculate_L_N - function(K, mat, params, S0, r, eps 1e-6, s 20) { # 计算四阶绝对矩 μ_n |φ^{(4)}(0)| mu_n - abs(phi4(0, mat, params, S0, r)) # 计算截断区间半宽 L公式来自 Junike (2024, Eq. (3.10)) L - (2 * K * exp(-r * mat) * mu_n / eps)^(1/4) # 计算误差界积分用于确定N integrand - function(u) { (1 / (2 * pi)) * (abs(u)^(s 1)) * abs(phi(u, mat, params, S0, r)) } # 数值积分注意积分区间为整个实数轴。实践中可截断到足够大的范围。 # 使用 integrate 函数设置适当的下限和上限。 boundDeriv - integrate(integrand, lower -100, upper 100, subdivisions 1000)$value # 注意积分限应足够大使得phi(u)衰减到接近0。对于Heston模型|u|很大时|phi(u)|衰减很快。 # 计算所需项数 N公式来自 Junike (2024, Sec. 6.1) numerator - 2^(s 5/2) * boundDeriv * L^(s 2) * 12 * K * exp(-r * mat) denominator - s * pi^(s 1) * eps N - ceiling((numerator / denominator)^(1 / s)) return(list(L L, N N, mu_n mu_n)) } # 10. 封装一个用户友好的定价函数 price_option_COS - function(K, mat, params, S0, r, type call, eps 1e-6) { # 步骤1计算该参数集下的最优L和N LN_info - calculate_L_N(K, mat, params, S0, r, eps) L - LN_info$L N - LN_info$N # 步骤2计算价格 if (type put) { price - put_COS(K, L, mat, N, params, S0, r) } else if (type call) { price - call_COS(K, L, mat, N, params, S0, r) } else { stop(Option type must be either call or put.) } # 返回价格和使用的参数用于调试和验证 return(list(price price, L L, N N, mu_n LN_info$mu_n)) }5. 实战测试、常见问题与性能优化5.1 完整测试用例与结果验证让我们用文献中常见的参数来测试我们的实现。这里我们使用与输入材料中相同的参数。# 设置测试参数 S0 - 100 # 标的现价 K - 90 # 行权价 r - 0.1 # 无风险利率 mat - 0.7 # 到期时间年 params - c( kappa 0.6067, # 均值回归速度 theta 0.0707, # 长期平均方差 xi 0.2928, # 波动率的波动率 rho -0.7571, # 相关系数 v0 0.0654 # 初始方差 ) # 设置误差容忍度 eps - 1e-6 # 计算看跌期权价格 result_put - price_option_COS(K, mat, params, S0, r, type put, eps eps) cat(sprintf(欧式看跌期权价格: %.6f\n, result_put$price)) cat(sprintf(使用的参数: L %.4f, N %d\n, result_put$L, result_put$N)) # 利用看跌-看涨平价计算看涨期权价格 call_price_by_parity - result_put$price S0 - K * exp(-r * mat) # 直接用COS方法计算看涨期权价格 result_call - price_option_COS(K, mat, params, S0, r, type call, eps eps) cat(sprintf(欧式看涨期权价格 (COS): %.6f\n, result_call$price)) cat(sprintf(欧式看涨期权价格 (平价): %.6f\n, call_price_by_parity)) cat(sprintf(两者差异: %.10f\n, abs(result_call$price - call_price_by_parity)))运行这段代码我们应该得到看跌期权价格约为2.773954这与输入材料中给出的结果一致。看涨期权的COS计算结果与通过看跌-看涨平价计算的结果应该只有极微小的数值误差远小于eps这验证了我们代码实现的一致性。5.2 常见问题排查与调试技巧在实现和运行过程中你可能会遇到以下问题结果为NaN或Inf检查特征函数首先验证psiLogST_Heston函数。对于较大的虚数参数u复数运算可能导致溢出。可以尝试在计算exp(d * mat)等项时检查中间变量d的实部是否为负确保衰减。Heston模型的特征函数在参数满足Feller条件2κθ ξ^2时通常表现良好。我们的测试参数满足该条件。检查积分区间在calculate_L_N函数中数值积分integrate的上下限(-100, 100)可能不够大。可以尝试增加上限如200, 500并观察积分值是否收敛。也可以先画出integrand(u)的函数图像看看它在多大范围内有显著非零值。计算速度慢瓶颈分析最耗时的部分是calculate_L_N中的数值积分和phi4的高阶符号求导。对于固定(params, mat, S0, r)的批量定价务必预计算并保存phi4函数和boundDeriv积分结果。向量化我们的ck和vk_put函数已经对k0:N进行了向量化运算这比用循环快得多。简化μ的计算如前所述用解析解mu_analytic替代通过Deriv得到的mu_via_deriv。精度不达标增大N如果结果与参考值偏差较大可以手动设置一个更大的N比如128, 256进行测试看看价格是否收敛。如果收敛到一个不同的值可能是L计算有误。检查L手动增大L比如乘以1.5或2如果价格发生显著变化说明初始L可能截断了密度函数的重要尾部。需要检查四阶矩mu_n的计算是否正确。验证特征函数与已知的Heston特征函数实现进行交叉验证例如与QuantLib或其他可靠库比较。可以计算psiLogST_Heston在一些u值上的输出。看跌-看涨平价不成立如果COS计算的看涨和看跌价格不满足平价关系问题几乎肯定出在vk_put函数或ck函数上。首先检查mu_analytic的计算是否正确与mu_via_deriv对比。然后用简单的数值积分方法如对密度函数采样计算一个基准价格来孤立问题。5.3 性能优化与生产环境建议对于需要高频或批量定价的生产环境可以考虑以下优化预计算与缓存这是最重要的优化。创建一个缓存系统为每个不同的参数组合(params, mat, S0, r)预计算并存储以下信息中心化均值mu_val四阶矩mu_n误差界积分值boundDeriv甚至可以为常用的(L, N)组合预计算好ck系数因为它不依赖于K。替换符号求导完全避免使用Deriv包。手动推导出phi(u)的1-4阶导数公式并实现为R函数。这需要一些复杂的微积分工作但能带来巨大的速度提升。使用更快的数值积分器integrate函数可能不是最快的。可以尝试使用更专业的数值积分库或者对于衰减很快的函数使用自适应高斯-埃尔米特求积法。并行计算如果需要对大量不同行权价K的期权定价由于每个K的计算是独立的可以很容易地用parallel包进行并行计算。代码移植对于极致性能要求可以考虑用C重写核心计算部分通过Rcpp包集成到R中特别是特征函数和系数计算的循环部分5.4 扩展与应用隐含波动率计算与模型校准一旦我们有了快速准确的定价器两个最直接的应用就是计算隐含波动率和模型校准。计算隐含波动率给定市场观察到的期权价格C_market我们需要找到使得Heston模型价格C_model(σ)等于市场价格的波动率参数这里其他Heston参数固定。由于COS定价器很快我们可以使用标准的求根算法如uniroot来反解。# 示例计算隐含波动率 (这里假设只变动v0其他参数固定) # 注意Heston模型的“隐含波动率”通常指与BS模型波动率σ的映射这里我们反解的是v0。 # 更常见的校准是反解所有参数。 calc_implied_v0 - function(market_price, K, mat, params_guess, S0, r, typecall) { # params_guess是包含v0初始猜测的参数向量 price_diff - function(v0_trial) { params_trial - params_guess params_trial[5] - v0_trial # 替换v0 model_price - price_option_COS(K, mat, params_trial, S0, r, typetype, eps1e-6)$price return(model_price - market_price) } # 使用二分法求根 root_result - uniroot(price_diff, interval c(0.001, 1.0)) # 给v0一个合理的搜索区间 return(root_result$root) }模型校准这是量化金融中的核心任务。给定一组不同行权价和到期日的期权市场价格我们寻找一组Heston模型参数(κ, θ, ξ, ρ, v0)使得模型价格与市场价格的总体误差如均方误差最小。这通常是一个高维非凸优化问题可以使用optim或nlminb等R优化函数结合我们的COS定价器来完成。# 简化的校准框架示例 calibrate_heston - function(market_data, params_init, S0, r) { # market_data: data.frame with columns: K, mat, market_price, type objective_function - function(params) { total_error - 0 for (i in 1:nrow(market_data)) { row - market_data[i, ] model_price - price_option_COS(row$K, row$mat, params, S0, r, typerow$type, eps1e-6)$price total_error - total_error (model_price - row$market_price)^2 } return(total_error) } # 添加参数约束如κ, θ, ξ 0, |ρ| 1, v0 0 lower_bounds - c(1e-3, 1e-3, 1e-3, -0.999, 1e-3) upper_bounds - c(10, 1, 5, 0.999, 1) result - optim(params_init, objective_function, method L-BFGS-B, lower lower_bounds, upper upper_bounds) return(result$par) }在校准过程中定价速度至关重要。因此之前提到的预计算和性能优化在这里会得到丰厚的回报。
