简介这是一份面向量化金融与随机波动率研究者的C实现资源主题为rBergomi粗糙波动率模型下的欧式期权定价。代码采用蒙特卡洛模拟用混合方案生成分数布朗运动fBm目前实现恒定前向方差情形适合研究rough volatility、期权定价算法及高性能计算的进阶开发者参考或二次开发。资源包共185个文件压缩后约62.37MB其中52个cpp与46个头文件构成核心定价引擎29个R脚本和Notebook用于交叉验证与实验分析PDF/Tex文档提供算法与编译说明目录结构清晰便于定位。目前已有175人浏览/学习这份资源。通过源码可学习sobol序列生成、蒙特卡洛路径模拟、rBergomi模型参数处理等关键环节附带的测试程序展示了拟随机数在路径模拟中的具体用法配套的R脚本与Notebook便于交叉验证也可移植到R环境或扩展为其他随机波动率模型。代码需要C14、FFTW3与OpenMP支持适宜有一定C基础的读者。1. 初识rBergomi一个H0.5的随机波动率模型能改变什么做期权定价的同行应该都有过这种体验用Heston或者经典SABR去拟合短期限的隐含波动率微笑ATM附近总是差一口气曲率不够陡、翅膀抬不起来。传统随机波动率模型里波动率路径太“平滑”时间序列上看涨跌绵绵根本复现不了真实市场中那种短时间尺度上的剧烈抖动。这几年量化金融里被反复讨论的rough volatility方向核心就是把这个“波动率的波动率”从H0.5改成H0.5让波动率路径本身变成一条粗糙的分数布朗运动。rBergomi模型就是这个方向里最直接、也最常被拿到生产环境里试的那一套。本文要拆的是一份C实现覆盖了rBergomi模型的方差过程模拟、蒙特卡洛期权定价和参数校准流水线。适合刚接触rough volatility、想找一套能编译能跑的C代码的量化开发者也适合已经在用经典随机波动率模型、想对比H0.5到底改进了什么的工程师。2. rBergomi的离散化:从分数布朗运动到可以上手的矩阵卷积2.1 模型再读一遍:方差过程为什么是指数形式rBergomi模型里方差过程的长相是v_t ξ_0 · exp( η · W_t^H − ½ η² · t^{2H} )其中W_t^H是Riemann-Liouville型分数布朗运动定义为W_t^H √(2H) ∫_0^t (t−s)^{H−1/2} dW_s^v这里的H就是Hurst指数rough volatility的核心参数。当H0.5时这个式子退化成普通布朗运动驱动的几何随机波动率过程;当H取0.05~0.15时路径的粗糙度明显上升。指数形式保证了v_t恒为正不需要像Heston那样担心Feller条件这是它工程上讨喜的原因之一。资产价格沿用经典的几何扩散形式:dS_t / S_t √v_t · dW_t^S其中dW_t^S与dW_s^v的瞬时相关系数为ρ。因为方差过程是路径依赖的整个模型没有解析解定价基本只能靠蒙特卡洛模拟。所以C实现好不好用主要就看两件事:一是分数布朗运动路径生成得够不够快二是蒙特卡洛循环的向量化程度高不高。2.2 核矩阵构造:把分数积分变成一次矩阵乘法离散化W_t^H时一个常见做法是把积分核在每个小区间上做闭式平均把Volterra积分表示成一个下三角核矩阵K与高斯向量z的乘积。设时间网格为0t_0t_1...t_NTΔtT/N离散形式为W_{t_k}^H ≈ √(2HΔt) · Σ_{i0}^{k} [ (t_k−t_i)^{H0.5} − (t_k−t_{i1})^{H0.5} ] / ((H0.5)·Δt) · z_i这里的z_i是独立标准正态随机变量。为什么用区间平均而不是直接取左端点核值?因为在H0.5时核函数(t−s)^{H−0.5}在s→t时是奇异的直接用左端点会喂给计算机一堆无穷大。区间平均的闭式形式天然消除了这个奇点数值上稳定很多。// 构造下三角核矩阵 K:每行k对应一个时间点t_k每列i对应一个驱动噪声z_i Eigen::MatrixXd buildKernel(int N, double T, double H) { const double dt T / N; Eigen::MatrixXd K Eigen::MatrixXd::Zero(N, N); const double a H 0.5; for (int k 0; k N; k) { const double tk (k 1) * dt; for (int i 0; i k; i) { const double r tk - i * dt; const double l tk - (i 1) * dt; // 区间平均核闭式积分避免端点奇异 K(k, i) std::sqrt(2.0 * H * dt) * (std::pow(r, a) - std::pow(l, a)) / (a * dt); } } return K; }这段代码的输出是一个N×N的下三角矩阵。构建完成后W_t^H的整条路径就是K乘以一个N维高斯向量。需要留意的是这里给K的每一行分配了不同的时间点行索引越大包含的历史噪声越多这正是Volterra过程的记忆特性。内存上用Eigen的MatrixXd就够N到500时矩阵只有约2MB,完全可接受。2.3 方差路径与资产价格路径的生成顺序有了W_t^H后方差路径是一个逐元素的指数映射:// 输入高斯噪声矩阵Z(M×N)输出方差路径矩阵(M×N) Eigen::MatrixXd simulateVolPaths(const Eigen::MatrixXd Z, const Eigen::MatrixXd K, double xi0, double eta, double H, double dt) { // 先算分式布朗运动路径:M×N Eigen::MatrixXd fbm Z * K.transpose(); // 构造确定性时间项:t^{2H} Eigen::VectorXd tVec Eigen::VectorXd::LinSpaced(K.rows(), 1, K.rows()) * dt; Eigen::ArrayXd tPow tVec.array().pow(2.0 * H); // 指数映射得到每条路径每个时刻的v_t Eigen::MatrixXd v xi0 * (eta * fbm.array() - 0.5 * eta * eta * tPow.replicate(1, Z.rows()).transpose()).exp().matrix(); return v; }这段代码里Z是一个M×N矩阵每行是一条独立路径的高斯白噪声序列。Z * K.transpose()是完整的矩阵乘法Eigen会调用BLAS级别的优化比逐条路径循环快得多。eta控制了波动率的波动率eta越大v_t的路径起伏越剧烈;xi0是初始方差水平校准中通常由市场ATM隐含波动率平方给出。资产价格路径的模拟要用到与方差路径相同的驱动噪声:// 模拟到期日资产价格Z驱动方差U是与Z独立的补充噪声 Eigen::VectorXd simulateTerminalSpot(const Eigen::MatrixXd v, const Eigen::MatrixXd Z, const Eigen::MatrixXd U, double S0, double r, double rho, double dt) { const int M v.rows(); const int N v.cols(); Eigen::MatrixXd S Eigen::MatrixXd::Zero(M, N 1); S.col(0).setConstant(S0); for (int j 0; j N; j) { // 相关性分解:资产噪声与方差噪声共享同一个正交分量 Eigen::MatrixXd eps rho * Z.col(j) std::sqrt(1.0 - rho * rho) * U.col(j); // Euler-Maruyama离散用前一时点的瞬时方差 S.col(j 1) S.col(j).array() * ((r - 0.5 * v.col(j).array()) * dt (v.col(j).array() * dt).sqrt() * eps.array()).exp(); } return S.col(N); }这里最关键的细节是:eps里的Z.col(j)必须和驱动该时段方差路径的噪声保持一致。如果图省事先用Z生成了v_t再重新抽一组独立噪声去模拟资产价格那资产与方差之间的相关性就被扔掉了ρ设成多少都不会生效。这个错误在初版实现里很容易犯后面避坑章会再展开。3. C实现骨架:模拟器、定价器与一条蒙特卡洛流水线的搭法3.1 随机数引擎的选择:普通生产场景不必迷信SobolrBergomi模拟涉及两个噪声源:驱动方差的高斯噪声Z和资产价格额外的补充噪声U。最直接的方案是std::mt19937加Box-Muller但Box-Muller的三角函数在M和N都很大时吞吐量一般。我一般会直接用Eigen的机制封装一个正态生成器或者用标准库normal_distribution取决于编译器的实现质量。对于低差异序列这里多说一句:Sobol序列在粗糙波动率模型里的优势并没有想象中大。原因是rBergomi的方差路径对高维噪声很敏感而Sobol在高维时的边际分布均匀性会被Box-Muller的非线性变换破坏低维投影反而出现条带。普通生产场景下paranoia级别调好的梅森旋转加ziggurat算法性价比更高。// 随机数引擎封装:每线程独立种子避免不同路径共享同一条随机流 class RngPool { public: explicit RngPool(int numThreads) { for (int t 0; t numThreads; t) { std::seed_seq seq{t, 0x9E3779B9}; engines_.emplace_back(seq); } } // 返回第i条路径对应线程的引擎引用 std::mt19937 get(int threadId) { return engines_[threadId]; } private: std::vectorstd::mt19937 engines_; };这个封装的要点是每条线程持有独立的随机流而不是所有线程共用一个全局引擎加锁。否则蒙特卡洛跑到一半线程竞争会让随机数的相关性变得一团糟路径之间出现肉眼可见的相似性。种子用线程ID加一个固定常数派生保证不同线程的随机流不重叠。3.2 蒙特卡洛定价器:一次模拟多个行权价rBergomi定价的常规做法是一次模拟出M条到期日资产价格路径然后对同一批路径计算多个行权价的payoff。这样核矩阵构建和方差路径生成都只做一遍边际成本大幅下降。定价器结构大致如下:struct RBergomiParams { double H 0.07; // Hurst指数越小越粗糙 double eta 1.3; // 波动率的波动率 double rho -0.9; // 资产与方差瞬时相关系数 double xi0 0.04; // 初始远期方差 double S0 100.0; // 初始资产价格 double r 0.02; // 无风险利率 double T 1.0; // 期限 int N 100; // 时间步数 int M 200000; // 路径数 }; class RBergomiPricer { public: std::vectordouble priceCallSmile(const RBergomiParams p, const std::vectordouble strikes, int seed) { double dt p.T / p.N; auto K buildKernel(p.N, p.T, p.H); Eigen::MatrixXd Z randn(p.M, p.N, seed); Eigen::MatrixXd U randn(p.M, p.N, seed 1); Eigen::MatrixXd v simulateVolPaths(Z, K, p.xi0, p.eta, p.H, dt); Eigen::VectorXd ST simulateTerminalSpot(v, Z, U, p.S0, p.r, p.rho, dt); std::vectordouble prices; prices.reserve(strikes.size()); for (double strike : strikes) { Eigen::VectorXd payoff (ST.array() - strike).max(0.0); prices.push_back(std::exp(-p.r * p.T) * payoff.mean()); } return prices; } private: Eigen::MatrixXd randn(int rows, int cols, int seed); };这个设计里M200000路径、N100步时Z矩阵是200000×100的double数组大约160MB。加上U和v峰值内存接近500MB。如果内存吃紧可以把路径分块处理每块50000条分四批跑。实际上分块还有一个额外好处:可以用OpenMP把四个块的pricing丢给四个线程线性增长吞吐。# 编译时建议开启O3和march nativeEigen的向量化依赖编译器的自动矢量化 g -O3 -marchnative -fopenmp -I/path/to/eigen rBergomi_mc.cpp -o rBergomi_mc发布版编译务必开O3Eigen的表达式模板在未优化编译下性能会掉一个量级。N100、M200000、10个行权价的典型场景开O3后单线程跑一次全流程大概在几秒到十几秒之间校准的时候能接受。3.3 内存布局与性能陷阱:矩阵乘法不是唯一瓶颈用Eigen做批量模拟时一个容易被忽略的点是v.col(j)按列访问。Eigen默认列主序存储S矩阵逐列写入时缓存友好。但如果你把路径数放在列、时间步放行访问模式就完全反了性能会明显下滑。我习惯统一把路径数放在行维度时间放在列维度避免每次取列都要隔M个double跳一次内存。还有一个小优化:核矩阵K可以不重复构建。在校准循环里H和T不变的情况下K只依赖N和H、T。如果H没变K完全可以缓存下来复用。很多初版实现每跑一次MC就重建一次KN500时会有肉眼可见的浪费。把buildKernel的结果存成成员变量校验H和N没变就直接返回。4. 参数校准:从Hurst指数到相关性四个核心参数的调试路径4.1 校准的数据准备:剥掉IV微笑的时间维度rBergomi的参数校准分两步:先拟合波动率曲面再拟合模型参数。实际生产中市场给出的通常是不同期限、不同行权价的一堆隐含波动率。第一步要把这些IV报价转换成一组的模型输入:每个期限对应一组行权价乖离度(K/S_0)和implied vol。常见做法是用SSVI或SVI参数化先平滑一遍曲面得到每个期限的ATM波动率σ_ATM(t)同时提取出初始远期方差曲线ξ_0(t)σ_ATM(t)²。这一步不做直接拿原始报价喂给优化器噪声会大到校准完全跑不动。4.2 目标函数与优化器选型:网格搜索做初值局部优化做精调校准目标函数写起来很直接:对每个期限T_k用当前参数跑一批蒙特卡洛路径算出模型隐含波动率与市场IV的差平方求和。由于rBergomi模拟本身带蒙特卡洛噪声目标函数不平滑直接用梯度下降会来回震荡。我一般先用粗网格搜索出一组过得去的初值再用无导数方法局部精调。// 网格初选:在H和eta的二维网格上跑小路径数MC找损失最小点 double bestLoss 1e30; double bestH 0.1, bestEta 1.0; for (int i 0; i 12; i) { for (int j 0; j 10; j) { double H 0.02 i * 0.02; // 0.02~0.24 double eta 0.5 j * 0.2; // 0.5~2.3 RBergomiParams p baseParams; p.H H; p.eta eta; double loss smileLoss(p, marketSmile); if (loss bestLoss) { bestLoss loss; bestH H; bestEta eta; } } } // 用网格最优解作初值跑一轮Nelder-Mead精调 RBergomiParams p baseParams; p.H bestH; p.eta bestEta; nelderMead(smileLoss, p, { {0.02, 0.3}, {0.8, 2.5}, {-0.98, 0.0} });这里有个工程习惯必须养成:网格搜索阶段每次调用MC必须固定随机种子否则这次算出来的loss和上次的loss之间不仅有参数差异还混了一坨随机噪声根本无法比较。我的做法是把全局seed设为固定值每次调用smileLoss前用seed预算ID派生新种子保证同一参数组合不同次调用也能重复。4.3 参数敏感性与边界行为:H和rho的实际可辨识性校准中最常遇到的现象是H和rho在短期限数据上相互扯皮。粗糙度影响微笑曲率rho影响微笑倾斜。当期限很短、行权价范围又窄时曲率和倾斜的信息混在一起H和rho的联合置信区域会拉得很长。四个参数的影响范围大致如下:参数含义典型范围主要影响HHurst指数控制路径粗糙度0.03~0.30短期限IV微笑曲率H越小曲率越陡eta波动率的波动率0.8~2.5整体微笑高度与曲率幅度rho资产与方差相关系数-0.98~0.0ATM附近倾斜方向与斜率xi0初始远期方差由市场ATM vol平方决定IV水平几乎不影响形状一个实用的起点是:先用不同期限的ATM vol把xi0(t)钉死然后固定rho-0.9先只调H和eta把微笑形状大概对上最后再放开rho精调。经验上rho对短期限微笑的影响远大于长期限如果你要跨期限一起校准权重分配要偏向短期限否则最终模型的期限结构会被长期限带走。5. 避坑指南:五个真实项目中踩过的rBergomi实现坑5.1 负方差:把指数形式拆成SDE欧拉离散现象:模拟出的v_t路径出现负值期权价格里出现NaN或明显偏大的值。原因:rBergomi的v_t本来就是指数形式保证正性。但有些实现为了沿用Heston的代码先把v_t取对数再拆成SDE然后用Euler-Maruyama离散。对数过程的离散会引入额外的离散误差而H0.5时的粗糙路径让误差放大得更快偶尔就会把指数内部顶到很大的负数进而出现溢出。解决:不要对v_t做任何额外的SDE展开。直接用闭式指数表达式v_tξ_0·exp(...)生成路径。核矩阵生成和指数映射这两步完全够了中间不要插入额外的随机微分方程求解器。5.2 核函数端点奇异:处处NaN很难排查现象:生成的fbm路径前几步正常跑到后半段全部变成NaN或Inf。原因:直接用左端点核(t_k−t_i)^{H−0.5}时当ik也就是核的自变量取到0而H0.5导致0的负次幂是无穷大。矩阵乘法碰到这个无穷就直接污染整条路径。解决:改成2.2节里的区间平均闭式核。它把(t_k−s)^{H−0.5}在每个小区间上做积分平均端点自洽不会出现除零或无穷。如果代码里已经用了左端点核可以把ik时那一项人为设成Δt^{H−0.5}也能跑通但精度不如闭式平均核。5.3 资产与方差噪声脱钩:ρ设了等于没设现象:把rho从-0.9调到0.5发现期权价格几乎不变IV smile的倾斜也纹丝不动。原因:资产价格模拟里的ε_k和驱动方差路径的z_k来自两组独立的随机变量两者其实完全没有相关性。常见的错误写法是为v_t生成一份随机数为S_t又生成一份随机数最后靠某个全局数组把两者“装成”相关的样子。解决:必须把驱动W_t^H的z_k直接复用到资产价格的随机项里用ε_kρ·z_k√(1−ρ²)·u_k构造相关性其中u_k是一组与z_k独立的标准正态变量。这个结构在3.2节的代码里已经体现改的时候要确认z_k不是复制副本而是同一个矩阵的同一列。5.4 校准中H跑到边界:短期限微笑拟合不足现象:H每次校准最终都停在0.02的下边界或者0.3的上边界而且短期限的IV微笑拟合残差始终很大。原因:单因子rBergomi的参数化空间对短期限微笑形状的控制有限。H控制了粗糙度但欠拟合时优化器会把H推到边界来强行增加曲率代价是长期限的微笑形状被破坏。另一个常见诱因是目标函数里短期限和长期限的权重没配好短期残差被长期限稀释。解决:检查分期限的拟合残差。如果短期限始终差先考虑给短期限加权重或者换用双因子rBergomi而非继续死磕单因子。如果H真的一直贴着边界说明这个参数组合本身就不适合当前市场状态别硬调。5.5 Sobol序列在高维MC里出现条带现象:用Sobol序列做低差异采样期权价格在路径数增大时反而出现振荡收敛速度比梅森旋转还差。原因:Sobol在高维时低维投影的均匀性会被Box-Muller的非线性变换打散生成的高斯噪声出现明显的网格状聚集。rBergomi的时间步N100对应100个方差驱动维度Sobol在这个维度上并没有优势。解决:普通生产场景直接用梅森旋转加ziggurat或标准正态分布生成器。如果确实想用低差异序列可以配合随机移位(Cranley-Patterson rotation)来打破条带但收益有限不值得为了这点精度去增加排查复杂度。6. 验证与进阶:经验Hurst检验、矩对比与双因子扩展拿到一套rBergomi实现后第一件事不是直接跑定价而是验证模拟器本身没写歪。我有一套固定的三件套检查每个新参数组合都要过一遍。第一步是经验Hurst估计:从模拟出的v_t路径里取对数增量统计|Δlog v_t|在不同时间尺度Δt上的均值再做log-log回归斜率再换算回的H应当与设定值差在±0.02以内。这一步能直接暴露核矩阵构造错误。第二步是矩对比:模型的理论方差Var(log v_t)η²·t^{2H}均值E[v_t]ξ_0。用模拟样本算一算偏差超过几个百分点就说明路径数太少或离散有问题。第三步是微笑冒烟测试:保持其他参数不变把H从0.07调到0.03再看看短期限IV曲线的曲率是否明显变陡。如果H怎么调曲率都纹丝不动说明模拟器中方差路径的时间尺度没接对。三件套过了再做双因子扩展。单因子rBergomi有个天生短板:短期限和长期限的微笑形态由同一个H控制经常顾此失彼。双因子版本直接把方差过程改成两个不同粗糙度的因子叠加:v_t ξ_0 · exp( η_1·W_t^{H_1} η_2·W_t^{H_2} − ½(η_1²·t^{2H_1} η_2²·t^{2H_2}) )实现上改动很小:把两个核矩阵K_1、K_2横向拼接成一个N×2N的矩阵驱动噪声Z从M×N扩成M×2N前半列驱动第一个因子后半列驱动第二个因子。资产价格的噪声仍然只用前半列的z_k来构造相关性。这样一次矩阵乘法就能同时生成两个粗糙度因子的路径收益是校准时的参数空间明显变大短期限与长期限的微笑可以解耦。从那以后我每次拿到新的参数组合都强制把三件套检查跑一遍再进校准流程。这三步加起来不过几百行代码却帮我挡住了多次核矩阵写错、路径配对错乱的返工。希望帮到你。本文还有配套的精品资源点击获取