简介rBergomi粗糙波动率模型的C实现面向量化金融领域从事期权定价、随机波动率建模的研究人员与中高级C开发者。项目以蒙特卡洛方法为核心通过混合方案生成分数布朗运动fBm完成欧式期权定价并实现了恒定前向方差场景为后续扩展随机方差等功能提供基础。压缩包共185个文件、约62.37MB以52个cpp、46个h源码文件为主配合29个r语言脚本和部分笔记本、文档、测试文件便于阅读、测试与二次开发。代码依赖fftw3和OpenMP需在支持C14的g5及以上环境中编译运行项目本身包含清晰的工程目录和多个测试用例。已有173人浏览学习适合需要参考完整蒙特卡洛定价实现与fBm混合方案代码的量化研究者和开发者。 最近在做一个随机波动率模型的参数校准项目Python原型跑得挺顺但一进入网格搜索就原形毕露——一次校准要上万次蒙特卡洛模拟等Python跑完咖啡都快凉了。于是下定决心把核心定价器用C重写而这一版的主角就是rBergomi模型。rBergomiRough Bergomi是金融工程里近年很火的一个粗糙波动率模型核心特点是Hurst参数远小于0.5能自然复现短期隐含波动率微笑的陡峭形态。这篇文章我把从公式到C代码的完整落地过程拆开讲包括路径模拟、核权重预计算、蒙特卡洛定价以及C工程里踩过的坑适合正在做量化开发、想用高性能语言实现随机过程模型的读者。1. rBergomi模型到底在模拟什么粗糙不是形容词是数学前提1.1 H 0.5 让波动率路径自带锯齿很多做衍生品定价的人对Bergomi模型不陌生它是经典的马尔可夫两因子随机波动率模型。rBergomi的差别看起来只是把Hurst参数从H0.5改成H0.5但这个改动在数学性质上完全是另一个物种。Hurst参数刻画的是路径的粗糙程度。H0.5对应标准布朗运动路径是连续的但处处不可微H0.5时路径变得更加锯齿化相邻增量的负相关性变强。放到波动率过程上H≈0.08意味着波动率的变化呈现出一种刚涨完马上回落的微观结构这正是市场短期波动率微笑的典型特征。传统随机波动率模型普遍难以同时拟合短端微笑的形状和长期偏度衰减速度而rBergomi用一个额外的分数阶参数就把这个问题解决了大半。代价也很直白过程不再是马尔可夫的因为方差过程依赖于整个历史的分数布朗运动路径这给数值实现带来了本质性的困难。1.2 核心方程V_t ξ(t) exp(η I_t - η²t^(2H)/2)风险中性测度下rBergomi模型通常写成两个关联方程资产价格方程dS_t / S_t sqrt(V_t) dW_t远期方差过程V_t ξ(t) * exp(η * I_t - 0.5 * η² * t^(2H))其中驱动的核心量是I_t sqrt(2H) * ∫₀ᵗ (t-s)^(H-1/2) dW_s这里的积分I_t是Riemann-Liouville分数阶积分的核心结构权重函数(t-s)^(H-1/2)是个幂律核。当H0.5时指数H-1/2为负靠近积分上限时核是发散的但仍然是可积的这种弱奇异性质正是Heston类模型里完全遇不到的。注意V_t公式里减掉的项是η²·t^(2H)/2而不是通常的η²·t/2。原因是I_t的方差正好是t^(2H)Var(I_t) 2H * ∫₀ᵗ (t-s)^(2H-1) ds t^(2H)所以指数里减掉的是半方差加上这个归一化项之后远期方差曲线的期望值刚好回到ξ(t)。这个细节实现了过拟合风险也解释了为什么网上很多错误实现算出来的平均方差路径歪得离谱就是因为少了一个维度上的归一化。2. 路径模拟的核心用离散卷积替掉连续积分2.1 高斯噪声路径是唯一的随机源整个模型只有一个随机源驱动I_t的标准布朗运动W_s。资产价格的布朗运动W_t和它相关相关系数用ρ表示。因此蒙特卡洛模拟的第一步是生成时间网格上的独立标准正态随机数序列记为ε_0, ε_1, ..., ε_{N-1}。假设时间网格是均匀的步长Δt T/N节点t_i i·Δt。离散化分数阶积分I_t时把每个小区间[t_j, t_{j1}]里的布朗增量近似为ε_j·sqrt(Δt)然后解析地计算核函数在该小区间上的积分∫_{t_j}^{t_{j1}} (t_i-s)^(H-1/2) ds令k i - j令u (s - t_j)/Δt则t_i - s (k-u)·Δt积分变成∫₀¹ (k-u)^(H-1/2) du · (Δt)^(H1/2) [(k)^(H1/2) - (k-1)^(H1/2)] / (H1/2) · (Δt)^(H1/2)因此离散化公式为I_{t_i} ≈ sqrt(2H) · (Δt)^(H1/2) · Σ_{j0}^{i-1} w_{i-j} · ε_j其中权重序列w_k [k^(H1/2) - (k-1)^(H1/2)] / (H1/2)这个形式非常关键权重只依赖于下标差i-j与具体的i和j无关。这意味着整条I_t路径是一个截断的Toeplitz卷积权重数组可以提前一次性算好。2.2 核权重数组为什么能预计算权重的预计算在构造rBergomiPricer对象时完成只要N和H固定w数组就是固定的。下面是C实现std::vectordouble computeKernelWeights(int N, double H) { std::vectordouble w(N); const double alpha H 0.5; // 注意这里是H0.5不是H-0.5 for (int k 1; k N; k) { w[k-1] (std::pow(k, alpha) - std::pow(k - 1.0, alpha)) / alpha; } return w; }这里的alpha是H0.5落在(0.5, 1)区间。当k1时w_0是0到1之间的一个小数它代表最近一个小区间的核积分贡献k越大权重越小因为核函数随距离衰减。实际模拟I_t路径时不是每次都递归调用pow而是用一个累积型循环std::vectordouble simulateConvolvedNoise(const std::vectordouble eps, const std::vectordouble w, int N, double H, double dt) { std::vectordouble I(N, 0.0); const double scale std::sqrt(2.0 * H) * std::pow(dt, H 0.5); for (int i 1; i N; i) { double acc 0.0; for (int j 0; j i; j) { acc w[i - j - 1] * eps[j]; } I[i] scale * acc; } return I; }这个双层循环是O(N²)的。N1024时单条路径的卷积约50万次乘加看起来不多但蒙特卡洛跑几十万条路径就变成几百亿次操作了所以后面必须优化。优化的思路我可以先剧透既然权重只依赖差这就是卷积结构FFT可以直接把单路径复杂度降到O(N log N)更激进的方案是用指数核近似变成O(M·N)M只需二三十后面详细说。2.3 这样实现出来的I_t协方差对吗用离散卷积模拟I_t最让人担心的就是它和理论I_t的统计性质是否一致答案是离散化带来的误差是O(Δt^(H1/2))量级网格越细越接近真实过程而且由于权重是按小区间解析积分算出来的比直接取中点值的欧拉离散要精确得多。验证方法很简单直接算模拟路径的自协方差函数和理论值对比。I_t和I_s的理论协方差为Cov(I_t, I_s) 2H · ∫₀^{min(t,s)} (t-u)^(H-1/2) · (s-u)^(H-1/2) du这个积分在网格上每个时间点都可以用数值积分算出来后面在验证章节我会给出具体的对比方式。3. C实现拆解一条完整的蒙特卡洛定价流水线3.1 随机数生成器别在细节上翻车蒙特卡洛定价最容易被轻视的就是随机数发生器。我的做法是每个线程持有独立的std::mt19937_64引擎并用std::normal_distribution生成标准正态随机数。多线程并行时每个线程独立构造引擎并配不同的种子序列避免数据竞争导致结果错误。一个常见的坑是老版本的编译器或平台对std::normal_distribution的实现方式不同同一套种子在不同平台生成的数据不一样。如果你需要跨平台复现比如回测框架在Windows开发、Linux部署强烈建议自己实现一个固定的Box-Muller或Ziggurat生成器并在测试里锁定一批随机数做回归对比。我自己的生产代码里用一个简单的Box-Muller实现额外缓存第二个高斯值避免浪费class NormalRNG { public: explicit NormalRNG(uint64_t seed) : engine_(seed) {} double next() { if (has_spare_) { has_spare_ false; return spare_; } double u1 uniform_(); double u2 uniform_(); double r std::sqrt(-2.0 * std::log(u1)); double theta 2.0 * M_PI * u2; spare_ r * std::sin(theta); has_spare_ true; return r * std::cos(theta); } private: std::mt19937_64 engine_; std::uniform_real_distributiondouble uniform_{0.0, 1.0}; double spare_ 0.0; bool has_spare_ false; };注意std::log和std::sqrt在每次调用时都有开销但这个开销在蒙特卡洛中的占比远小于卷积计算不用在这里过度优化。3.2 Forward Variance来自市场数据的输入模型里的ξ(t)是远期方差曲线表达的是期限结构信息。最简单的设定是常数ξ(t)V₀即初始瞬时方差处处相等。更贴近市场的做法是从SSVI或SABR参数化拟合出来的方差期限结构插值。我在代码里把ξ(t)抽象成一个接口class ForwardVarianceCurve { public: virtual ~ForwardVarianceCurve() default; virtual double value(double t) const 0; }; class ConstantForwardVariance : public ForwardVarianceCurve { public: explicit ConstantForwardVariance(double v0) : v0_(v0) {} double value(double /*t*/) const override { return v0_; } private: double v0_; };这样后续想换成插值曲线不需要改动模拟主流程。3.3 一条路径从头走到尾方差过程的递推与存储有了随机数和核权重模拟一条完整路径的流程就清晰了。核心是一个模拟引擎类struct RoughBergomiPath { std::vectordouble variance; std::vectordouble spot; double terminalSpot 0.0; }; class RoughBergomiEngine { public: RoughBergomiEngine(int N, double T, double H, double eta, double rho, const ForwardVarianceCurve fwdVar) : N_(N), T_(T), H_(H), eta_(eta), rho_(rho), dt_(T / N), fwdVar_(fwdVar), weights_(computeKernelWeights(N, H)) {} RoughBergomiPath simulatePath(NormalRNG rng) const { std::vectordouble eps(N_); for (int i 0; i N_; i) { eps[i] rng.next(); } // 1. 卷积得到 I_t std::vectordouble I convolve(eps); // 2. 计算方差路径和资产路径 RoughBergomiPath path; path.variance.resize(N_); path.spot.resize(N_); path.spot[0] S0_; double s S0_; for (int i 1; i N_; i) { double t i * dt_; double var fwdVar_.value(t) * std::exp(eta_ * I[i] - 0.5 * eta_ * eta_ * std::pow(t, 2.0 * H_)); path.variance[i] var; double zS rho_ * eps[i-1] std::sqrt(1.0 - rho_ * rho_) * rng.next(); s * std::exp((riskFreeRate_ - 0.5 * var) * dt_ std::sqrt(var * dt_) * zS); path.spot[i] s; } path.terminalSpot s; return path; } private: int N_; double T_, H_, eta_, rho_, dt_; double S0_ 100.0; double riskFreeRate_ 0.0; const ForwardVarianceCurve fwdVar_; std::vectordouble weights_; };注意资产路径的随机项zS必须复用第i-1个eps[i-1]来保证相关系数ρ这里就是教科书之外最容易写错的地方。很多人先把I_t完整模拟出来再用一个全新独立的高斯变量驱动资产结果相关系数完全失效。3.4 蒙特卡洛定价器与置信区间定价器只需要不断调用模拟引擎、统计贴现后的收益。单条路径的内存开销是O(N)路径之间完全独立天然适合并行double priceEuropeanCall(const RoughBergomiEngine engine, double strike, int numPaths, int numThreads) { std::atomicdouble sumPayoff{0.0}; std::atomicdouble sumPayoffSq{0.0}; #pragma omp parallel num_threads(numThreads) { uint64_t seed 0x9E3779B97F4A7C15ULL ^ static_castuint64_t(omp_get_thread_num() 1); NormalRNG rng(seed); double localSum 0.0; double localSq 0.0; #pragma omp for schedule(static) for (int p 0; p numPaths; p) { auto path engine.simulatePath(rng); double payoff std::max(path.terminalSpot - strike, 0.0); localSum payoff; localSq payoff * payoff; } sumPayoff localSum; sumPayoffSq localSq; } double mean sumPayoff / numPaths; double sd std::sqrt((sumPayoffSq / numPaths - mean * mean) / (numPaths - 1)); return mean; // 这里可以同时输出sd/sqrt(numPaths)作为标准误 }这种实现方式的好处是大规模并行时每个线程只维护自己的局部和最后用原子操作合并避免频繁同步。同时每条路径内部用流式方式模拟不需要把整批路径都存在内存里对缓存非常友好。4. 数值验证三件套自协方差、退化测试、文献对照4.1 自协方差与理论值的对比任何随机过程模拟代码写完后第一件事不是看期权价格而是验证路径的统计分布是否正确。对我而言最有效的验证是模拟10万条路径计算I_t在固定时间点t0.5, 1.0上的方差和协方差与理论公式对比。理论方差是Var(I_t) t^(2H)。取H0.08, t1.0时理论方差1^(0.16)1.0H0.08, t0.5时理论方差0.5^(0.16)≈0.895。模拟结果的偏差应在1%以内才算通过。协方差验证稍微复杂一点我用自适应Simpson积分算出理论值再对比。如果偏差超过几个百分点优先检查权重数组的索引是否偏移了。4.2 H→0.5退化测试rBergomi在H0.5时应该退化成标准Bergomi模型或者说马尔可夫版本。这是最廉价的回归测试把H设成0.5跑一组固定路径算出来的价格应该和一个用标准欧拉离散的Heston模型当相关参数对齐时几乎一致。这个测试能抓出一大类实现错误。比如指数归一化项如果把t^(2H)误写成tH0.5时恰好没区别但在H0.05时就会产生巨大的偏差。反过来如果H0.5时结果就对不上说明核权重数组的基本公式就有问题。4.3 复现论文里的期权价格表最后一步是拿公开文献里的数值结果做对标。Bayer, Friz, Gatheral在2016年的论文里给了不同Hurst参数、不同到期日下的期权价格表选了H0.07, 0.1等参数组合。复现时不需要完全一致因为他们的网格设置和随机数种子不同但价格差异应该在蒙特卡洛标准误的三倍以内。我第一次跑复现时ATM期权价格对上了但OTM期权差了将近2%后来发现是资产路径的离散化用了对数值还是原值的问题——使用几何布朗运动的解析条件就得用对数形式直接对S做欧拉离散在波动率高时会引入显著偏差。5. 为什么我不建议一上来用Cholesky以及生产环境的加速方案5.1 Cholesky方法的三重麻烦看到I_t的协方差矩阵很多人第一反应是用Cholesky分解来精确生成I_t向量计算N×N协方差矩阵C分解成CLL^T然后ILz。这个方法在N128时完全没问题但实际定价时N至少要1024此时协方差矩阵是100万个数Cholesky分解的计算量是O(N³)光分解就要几十秒。更大的问题是Cholesky生成的是I_t的边际联合分布它丢失了与驱动布朗运动W_s路径增量之间的对应关系。而资产价格路径恰恰需要ρ·ε_i这层耦合你没法从I_t反推出和它线性相关的ε序列。当然也可以把ε^H和I_t拼成一个2N维的联合高斯向量再做Cholesky但2N维矩阵分解的内存和耗时就更不现实了。直接卷积法之所以是事实标准就是因为它天然保留了完整耦合结构同时O(N²)的计算量在N1024时可以接受。所以我的结论是教学、验证阶段用Cholesky可以生产环境第一版就用卷积法。5.2 从O(N²)到O(M·N)指数核近似生产环境性能不够时最先考虑的不是换语言而是换算法。既然权重函数(t-s)^(H-1/2)是幂律核而幂律核可以用指数函数和逼近(t-s)^(H-1/2) ≈ Σ_{k1}^{M} a_k · exp(-b_k·(t-s))这项技术来自Gatheral等人在2018年提出的Markovian approximation思想。每个指数项对应一个一阶线性SDEdX_k(t) -b_k·X_k(t)·dt a_k·dW(t)在均匀网格上可以直接递推x_k * std::exp(-b_k * dt); x_k a_k * eps[i] * std::sqrt(dt);M取20到30项时就能达到很高的逼近精度。这样生成整条I_t路径的复杂度从O(N²)降到了O(M·N)在N2048时大致有60到100倍的加速而且不用每条路径都存权重矩阵。代价是要先把a_k和b_k算出来这通常离线完成。网上有开源的预计算参数表也可以自己用最小二乘拟合。理论上这套近似不是精确的分数布朗运动但蒙特卡洛的统计误差通常远大于近似误差所以工程上完全够用。5.3 与Python原型的性能对比我自己的对比经验是Python用NumPy向量化实现同样的卷积方案模拟10万条路径N256大约需要几十秒C单线程版本大概1到2秒加上OpenMP八线程并行后能压到0.2到0.3秒。再配上指数核近似整体可以做到比原始Python快两个量级以上。这个提升对参数校准是质变原来一个网格点要等一分钟现在几秒就出结果整个校准流程从过夜任务变成午休任务。6. 编译和工程上的四个坑6.1 浮点幂运算的隐藏成本在模拟循环里最容易被忽略的昂贵操作就是std::pow。计算V_t时每次都要调std::pow(t, 2*H)H又是小数这个调用比普通的加减乘除慢几十倍。一条路径N2048就要调2048次pow几十万条路径就是几亿次pow。解法很简单在路径模拟前预计算每个网格点的pow(t_i, 2*H)数组查询就行。核权重数组也是一样的道理N固定的话提前算好不要在每条路径里重复计算。6.2 OpenMP并行时的随机数种子陷阱很多人并行蒙特卡洛时图省事所有线程共用一个随机数引擎加锁或者用固定的种子加上线程号。前者的性能会暴跌后者可能因为两个线程的噪声序列相关导致方差低估。我用的方案是主线程生成一个64位种子每个工作线程用seed ^ (threadId * 大质数)做种子。注意不同线程之间不能复用完全相同的种子。另外std::normal_distribution里可能有内部状态每个线程必须持有一份独立的分布对象。6.3-ffast-math会悄悄破坏你的金融计算很多追求性能的C开发者在Release编译时习惯加-ffast-math。这个选项会允许编译器做不安全的浮点重排、忽略NaN和Inf的检查对蒙特卡洛算价格来说问题不大因为少量数值误差会被统计平均掉。但如果你在做协方差验证、或者要用价格梯度做校准-ffast-math可能带来不可控的精度损失。我的做法是开发版不开任何激进浮点优化生产版只开-O3 -marchnative保留IEEE浮点语义。实测性能差异在5%以内远小于排查错误的时间成本。6.4 内存布局路径流式比存储整批快得多刚开始写并行版本时我图省事把模拟出的所有路径都存进一个大数组最后统一算价格。结果N2048、路径数20万时内存占用超过了3GB而且缓存命中率很差。改成每路径内部直接算payoff、只保留一条路径状态的流式方案后内存占用降到几MB速度反而因为缓存友好而提升了。只有当你要用对偶变量法或重要性抽样等方差缩减技术时才有必要保留完整的路径信息。7. 我自己的工程体会rBergomi模型的C实现真正难的地方不在C语法而在两处一是把分数阶积分和弱奇异核的数学性质吃透二是把蒙特卡洛的性能瓶颈分析清楚。数学理解不到位写出来的代码看起来在模拟rBergomi实际统计性质全是错的性能不做算法层面的优化换什么语言都只能能用但难用。如果你准备在量化项目里落地rBergomi我建议按这个顺序动手先实现朴素卷积版把验证三件套跑通再考虑指数核近似或FFT加速最后做OpenMP并行和编译优化。每一步都有明确的验证手段不会出现改完性能上去了但不知道结果对不对的尴尬局面。本文还有配套的精品资源点击获取 rBergomi随机波动率模型的C++高性能实现与蒙特卡洛定价
简介rBergomi粗糙波动率模型的C实现面向量化金融领域从事期权定价、随机波动率建模的研究人员与中高级C开发者。项目以蒙特卡洛方法为核心通过混合方案生成分数布朗运动fBm完成欧式期权定价并实现了恒定前向方差场景为后续扩展随机方差等功能提供基础。压缩包共185个文件、约62.37MB以52个cpp、46个h源码文件为主配合29个r语言脚本和部分笔记本、文档、测试文件便于阅读、测试与二次开发。代码依赖fftw3和OpenMP需在支持C14的g5及以上环境中编译运行项目本身包含清晰的工程目录和多个测试用例。已有173人浏览学习适合需要参考完整蒙特卡洛定价实现与fBm混合方案代码的量化研究者和开发者。 最近在做一个随机波动率模型的参数校准项目Python原型跑得挺顺但一进入网格搜索就原形毕露——一次校准要上万次蒙特卡洛模拟等Python跑完咖啡都快凉了。于是下定决心把核心定价器用C重写而这一版的主角就是rBergomi模型。rBergomiRough Bergomi是金融工程里近年很火的一个粗糙波动率模型核心特点是Hurst参数远小于0.5能自然复现短期隐含波动率微笑的陡峭形态。这篇文章我把从公式到C代码的完整落地过程拆开讲包括路径模拟、核权重预计算、蒙特卡洛定价以及C工程里踩过的坑适合正在做量化开发、想用高性能语言实现随机过程模型的读者。1. rBergomi模型到底在模拟什么粗糙不是形容词是数学前提1.1 H 0.5 让波动率路径自带锯齿很多做衍生品定价的人对Bergomi模型不陌生它是经典的马尔可夫两因子随机波动率模型。rBergomi的差别看起来只是把Hurst参数从H0.5改成H0.5但这个改动在数学性质上完全是另一个物种。Hurst参数刻画的是路径的粗糙程度。H0.5对应标准布朗运动路径是连续的但处处不可微H0.5时路径变得更加锯齿化相邻增量的负相关性变强。放到波动率过程上H≈0.08意味着波动率的变化呈现出一种刚涨完马上回落的微观结构这正是市场短期波动率微笑的典型特征。传统随机波动率模型普遍难以同时拟合短端微笑的形状和长期偏度衰减速度而rBergomi用一个额外的分数阶参数就把这个问题解决了大半。代价也很直白过程不再是马尔可夫的因为方差过程依赖于整个历史的分数布朗运动路径这给数值实现带来了本质性的困难。1.2 核心方程V_t ξ(t) exp(η I_t - η²t^(2H)/2)风险中性测度下rBergomi模型通常写成两个关联方程资产价格方程dS_t / S_t sqrt(V_t) dW_t远期方差过程V_t ξ(t) * exp(η * I_t - 0.5 * η² * t^(2H))其中驱动的核心量是I_t sqrt(2H) * ∫₀ᵗ (t-s)^(H-1/2) dW_s这里的积分I_t是Riemann-Liouville分数阶积分的核心结构权重函数(t-s)^(H-1/2)是个幂律核。当H0.5时指数H-1/2为负靠近积分上限时核是发散的但仍然是可积的这种弱奇异性质正是Heston类模型里完全遇不到的。注意V_t公式里减掉的项是η²·t^(2H)/2而不是通常的η²·t/2。原因是I_t的方差正好是t^(2H)Var(I_t) 2H * ∫₀ᵗ (t-s)^(2H-1) ds t^(2H)所以指数里减掉的是半方差加上这个归一化项之后远期方差曲线的期望值刚好回到ξ(t)。这个细节实现了过拟合风险也解释了为什么网上很多错误实现算出来的平均方差路径歪得离谱就是因为少了一个维度上的归一化。2. 路径模拟的核心用离散卷积替掉连续积分2.1 高斯噪声路径是唯一的随机源整个模型只有一个随机源驱动I_t的标准布朗运动W_s。资产价格的布朗运动W_t和它相关相关系数用ρ表示。因此蒙特卡洛模拟的第一步是生成时间网格上的独立标准正态随机数序列记为ε_0, ε_1, ..., ε_{N-1}。假设时间网格是均匀的步长Δt T/N节点t_i i·Δt。离散化分数阶积分I_t时把每个小区间[t_j, t_{j1}]里的布朗增量近似为ε_j·sqrt(Δt)然后解析地计算核函数在该小区间上的积分∫_{t_j}^{t_{j1}} (t_i-s)^(H-1/2) ds令k i - j令u (s - t_j)/Δt则t_i - s (k-u)·Δt积分变成∫₀¹ (k-u)^(H-1/2) du · (Δt)^(H1/2) [(k)^(H1/2) - (k-1)^(H1/2)] / (H1/2) · (Δt)^(H1/2)因此离散化公式为I_{t_i} ≈ sqrt(2H) · (Δt)^(H1/2) · Σ_{j0}^{i-1} w_{i-j} · ε_j其中权重序列w_k [k^(H1/2) - (k-1)^(H1/2)] / (H1/2)这个形式非常关键权重只依赖于下标差i-j与具体的i和j无关。这意味着整条I_t路径是一个截断的Toeplitz卷积权重数组可以提前一次性算好。2.2 核权重数组为什么能预计算权重的预计算在构造rBergomiPricer对象时完成只要N和H固定w数组就是固定的。下面是C实现std::vectordouble computeKernelWeights(int N, double H) { std::vectordouble w(N); const double alpha H 0.5; // 注意这里是H0.5不是H-0.5 for (int k 1; k N; k) { w[k-1] (std::pow(k, alpha) - std::pow(k - 1.0, alpha)) / alpha; } return w; }这里的alpha是H0.5落在(0.5, 1)区间。当k1时w_0是0到1之间的一个小数它代表最近一个小区间的核积分贡献k越大权重越小因为核函数随距离衰减。实际模拟I_t路径时不是每次都递归调用pow而是用一个累积型循环std::vectordouble simulateConvolvedNoise(const std::vectordouble eps, const std::vectordouble w, int N, double H, double dt) { std::vectordouble I(N, 0.0); const double scale std::sqrt(2.0 * H) * std::pow(dt, H 0.5); for (int i 1; i N; i) { double acc 0.0; for (int j 0; j i; j) { acc w[i - j - 1] * eps[j]; } I[i] scale * acc; } return I; }这个双层循环是O(N²)的。N1024时单条路径的卷积约50万次乘加看起来不多但蒙特卡洛跑几十万条路径就变成几百亿次操作了所以后面必须优化。优化的思路我可以先剧透既然权重只依赖差这就是卷积结构FFT可以直接把单路径复杂度降到O(N log N)更激进的方案是用指数核近似变成O(M·N)M只需二三十后面详细说。2.3 这样实现出来的I_t协方差对吗用离散卷积模拟I_t最让人担心的就是它和理论I_t的统计性质是否一致答案是离散化带来的误差是O(Δt^(H1/2))量级网格越细越接近真实过程而且由于权重是按小区间解析积分算出来的比直接取中点值的欧拉离散要精确得多。验证方法很简单直接算模拟路径的自协方差函数和理论值对比。I_t和I_s的理论协方差为Cov(I_t, I_s) 2H · ∫₀^{min(t,s)} (t-u)^(H-1/2) · (s-u)^(H-1/2) du这个积分在网格上每个时间点都可以用数值积分算出来后面在验证章节我会给出具体的对比方式。3. C实现拆解一条完整的蒙特卡洛定价流水线3.1 随机数生成器别在细节上翻车蒙特卡洛定价最容易被轻视的就是随机数发生器。我的做法是每个线程持有独立的std::mt19937_64引擎并用std::normal_distribution生成标准正态随机数。多线程并行时每个线程独立构造引擎并配不同的种子序列避免数据竞争导致结果错误。一个常见的坑是老版本的编译器或平台对std::normal_distribution的实现方式不同同一套种子在不同平台生成的数据不一样。如果你需要跨平台复现比如回测框架在Windows开发、Linux部署强烈建议自己实现一个固定的Box-Muller或Ziggurat生成器并在测试里锁定一批随机数做回归对比。我自己的生产代码里用一个简单的Box-Muller实现额外缓存第二个高斯值避免浪费class NormalRNG { public: explicit NormalRNG(uint64_t seed) : engine_(seed) {} double next() { if (has_spare_) { has_spare_ false; return spare_; } double u1 uniform_(); double u2 uniform_(); double r std::sqrt(-2.0 * std::log(u1)); double theta 2.0 * M_PI * u2; spare_ r * std::sin(theta); has_spare_ true; return r * std::cos(theta); } private: std::mt19937_64 engine_; std::uniform_real_distributiondouble uniform_{0.0, 1.0}; double spare_ 0.0; bool has_spare_ false; };注意std::log和std::sqrt在每次调用时都有开销但这个开销在蒙特卡洛中的占比远小于卷积计算不用在这里过度优化。3.2 Forward Variance来自市场数据的输入模型里的ξ(t)是远期方差曲线表达的是期限结构信息。最简单的设定是常数ξ(t)V₀即初始瞬时方差处处相等。更贴近市场的做法是从SSVI或SABR参数化拟合出来的方差期限结构插值。我在代码里把ξ(t)抽象成一个接口class ForwardVarianceCurve { public: virtual ~ForwardVarianceCurve() default; virtual double value(double t) const 0; }; class ConstantForwardVariance : public ForwardVarianceCurve { public: explicit ConstantForwardVariance(double v0) : v0_(v0) {} double value(double /*t*/) const override { return v0_; } private: double v0_; };这样后续想换成插值曲线不需要改动模拟主流程。3.3 一条路径从头走到尾方差过程的递推与存储有了随机数和核权重模拟一条完整路径的流程就清晰了。核心是一个模拟引擎类struct RoughBergomiPath { std::vectordouble variance; std::vectordouble spot; double terminalSpot 0.0; }; class RoughBergomiEngine { public: RoughBergomiEngine(int N, double T, double H, double eta, double rho, const ForwardVarianceCurve fwdVar) : N_(N), T_(T), H_(H), eta_(eta), rho_(rho), dt_(T / N), fwdVar_(fwdVar), weights_(computeKernelWeights(N, H)) {} RoughBergomiPath simulatePath(NormalRNG rng) const { std::vectordouble eps(N_); for (int i 0; i N_; i) { eps[i] rng.next(); } // 1. 卷积得到 I_t std::vectordouble I convolve(eps); // 2. 计算方差路径和资产路径 RoughBergomiPath path; path.variance.resize(N_); path.spot.resize(N_); path.spot[0] S0_; double s S0_; for (int i 1; i N_; i) { double t i * dt_; double var fwdVar_.value(t) * std::exp(eta_ * I[i] - 0.5 * eta_ * eta_ * std::pow(t, 2.0 * H_)); path.variance[i] var; double zS rho_ * eps[i-1] std::sqrt(1.0 - rho_ * rho_) * rng.next(); s * std::exp((riskFreeRate_ - 0.5 * var) * dt_ std::sqrt(var * dt_) * zS); path.spot[i] s; } path.terminalSpot s; return path; } private: int N_; double T_, H_, eta_, rho_, dt_; double S0_ 100.0; double riskFreeRate_ 0.0; const ForwardVarianceCurve fwdVar_; std::vectordouble weights_; };注意资产路径的随机项zS必须复用第i-1个eps[i-1]来保证相关系数ρ这里就是教科书之外最容易写错的地方。很多人先把I_t完整模拟出来再用一个全新独立的高斯变量驱动资产结果相关系数完全失效。3.4 蒙特卡洛定价器与置信区间定价器只需要不断调用模拟引擎、统计贴现后的收益。单条路径的内存开销是O(N)路径之间完全独立天然适合并行double priceEuropeanCall(const RoughBergomiEngine engine, double strike, int numPaths, int numThreads) { std::atomicdouble sumPayoff{0.0}; std::atomicdouble sumPayoffSq{0.0}; #pragma omp parallel num_threads(numThreads) { uint64_t seed 0x9E3779B97F4A7C15ULL ^ static_castuint64_t(omp_get_thread_num() 1); NormalRNG rng(seed); double localSum 0.0; double localSq 0.0; #pragma omp for schedule(static) for (int p 0; p numPaths; p) { auto path engine.simulatePath(rng); double payoff std::max(path.terminalSpot - strike, 0.0); localSum payoff; localSq payoff * payoff; } sumPayoff localSum; sumPayoffSq localSq; } double mean sumPayoff / numPaths; double sd std::sqrt((sumPayoffSq / numPaths - mean * mean) / (numPaths - 1)); return mean; // 这里可以同时输出sd/sqrt(numPaths)作为标准误 }这种实现方式的好处是大规模并行时每个线程只维护自己的局部和最后用原子操作合并避免频繁同步。同时每条路径内部用流式方式模拟不需要把整批路径都存在内存里对缓存非常友好。4. 数值验证三件套自协方差、退化测试、文献对照4.1 自协方差与理论值的对比任何随机过程模拟代码写完后第一件事不是看期权价格而是验证路径的统计分布是否正确。对我而言最有效的验证是模拟10万条路径计算I_t在固定时间点t0.5, 1.0上的方差和协方差与理论公式对比。理论方差是Var(I_t) t^(2H)。取H0.08, t1.0时理论方差1^(0.16)1.0H0.08, t0.5时理论方差0.5^(0.16)≈0.895。模拟结果的偏差应在1%以内才算通过。协方差验证稍微复杂一点我用自适应Simpson积分算出理论值再对比。如果偏差超过几个百分点优先检查权重数组的索引是否偏移了。4.2 H→0.5退化测试rBergomi在H0.5时应该退化成标准Bergomi模型或者说马尔可夫版本。这是最廉价的回归测试把H设成0.5跑一组固定路径算出来的价格应该和一个用标准欧拉离散的Heston模型当相关参数对齐时几乎一致。这个测试能抓出一大类实现错误。比如指数归一化项如果把t^(2H)误写成tH0.5时恰好没区别但在H0.05时就会产生巨大的偏差。反过来如果H0.5时结果就对不上说明核权重数组的基本公式就有问题。4.3 复现论文里的期权价格表最后一步是拿公开文献里的数值结果做对标。Bayer, Friz, Gatheral在2016年的论文里给了不同Hurst参数、不同到期日下的期权价格表选了H0.07, 0.1等参数组合。复现时不需要完全一致因为他们的网格设置和随机数种子不同但价格差异应该在蒙特卡洛标准误的三倍以内。我第一次跑复现时ATM期权价格对上了但OTM期权差了将近2%后来发现是资产路径的离散化用了对数值还是原值的问题——使用几何布朗运动的解析条件就得用对数形式直接对S做欧拉离散在波动率高时会引入显著偏差。5. 为什么我不建议一上来用Cholesky以及生产环境的加速方案5.1 Cholesky方法的三重麻烦看到I_t的协方差矩阵很多人第一反应是用Cholesky分解来精确生成I_t向量计算N×N协方差矩阵C分解成CLL^T然后ILz。这个方法在N128时完全没问题但实际定价时N至少要1024此时协方差矩阵是100万个数Cholesky分解的计算量是O(N³)光分解就要几十秒。更大的问题是Cholesky生成的是I_t的边际联合分布它丢失了与驱动布朗运动W_s路径增量之间的对应关系。而资产价格路径恰恰需要ρ·ε_i这层耦合你没法从I_t反推出和它线性相关的ε序列。当然也可以把ε^H和I_t拼成一个2N维的联合高斯向量再做Cholesky但2N维矩阵分解的内存和耗时就更不现实了。直接卷积法之所以是事实标准就是因为它天然保留了完整耦合结构同时O(N²)的计算量在N1024时可以接受。所以我的结论是教学、验证阶段用Cholesky可以生产环境第一版就用卷积法。5.2 从O(N²)到O(M·N)指数核近似生产环境性能不够时最先考虑的不是换语言而是换算法。既然权重函数(t-s)^(H-1/2)是幂律核而幂律核可以用指数函数和逼近(t-s)^(H-1/2) ≈ Σ_{k1}^{M} a_k · exp(-b_k·(t-s))这项技术来自Gatheral等人在2018年提出的Markovian approximation思想。每个指数项对应一个一阶线性SDEdX_k(t) -b_k·X_k(t)·dt a_k·dW(t)在均匀网格上可以直接递推x_k * std::exp(-b_k * dt); x_k a_k * eps[i] * std::sqrt(dt);M取20到30项时就能达到很高的逼近精度。这样生成整条I_t路径的复杂度从O(N²)降到了O(M·N)在N2048时大致有60到100倍的加速而且不用每条路径都存权重矩阵。代价是要先把a_k和b_k算出来这通常离线完成。网上有开源的预计算参数表也可以自己用最小二乘拟合。理论上这套近似不是精确的分数布朗运动但蒙特卡洛的统计误差通常远大于近似误差所以工程上完全够用。5.3 与Python原型的性能对比我自己的对比经验是Python用NumPy向量化实现同样的卷积方案模拟10万条路径N256大约需要几十秒C单线程版本大概1到2秒加上OpenMP八线程并行后能压到0.2到0.3秒。再配上指数核近似整体可以做到比原始Python快两个量级以上。这个提升对参数校准是质变原来一个网格点要等一分钟现在几秒就出结果整个校准流程从过夜任务变成午休任务。6. 编译和工程上的四个坑6.1 浮点幂运算的隐藏成本在模拟循环里最容易被忽略的昂贵操作就是std::pow。计算V_t时每次都要调std::pow(t, 2*H)H又是小数这个调用比普通的加减乘除慢几十倍。一条路径N2048就要调2048次pow几十万条路径就是几亿次pow。解法很简单在路径模拟前预计算每个网格点的pow(t_i, 2*H)数组查询就行。核权重数组也是一样的道理N固定的话提前算好不要在每条路径里重复计算。6.2 OpenMP并行时的随机数种子陷阱很多人并行蒙特卡洛时图省事所有线程共用一个随机数引擎加锁或者用固定的种子加上线程号。前者的性能会暴跌后者可能因为两个线程的噪声序列相关导致方差低估。我用的方案是主线程生成一个64位种子每个工作线程用seed ^ (threadId * 大质数)做种子。注意不同线程之间不能复用完全相同的种子。另外std::normal_distribution里可能有内部状态每个线程必须持有一份独立的分布对象。6.3-ffast-math会悄悄破坏你的金融计算很多追求性能的C开发者在Release编译时习惯加-ffast-math。这个选项会允许编译器做不安全的浮点重排、忽略NaN和Inf的检查对蒙特卡洛算价格来说问题不大因为少量数值误差会被统计平均掉。但如果你在做协方差验证、或者要用价格梯度做校准-ffast-math可能带来不可控的精度损失。我的做法是开发版不开任何激进浮点优化生产版只开-O3 -marchnative保留IEEE浮点语义。实测性能差异在5%以内远小于排查错误的时间成本。6.4 内存布局路径流式比存储整批快得多刚开始写并行版本时我图省事把模拟出的所有路径都存进一个大数组最后统一算价格。结果N2048、路径数20万时内存占用超过了3GB而且缓存命中率很差。改成每路径内部直接算payoff、只保留一条路径状态的流式方案后内存占用降到几MB速度反而因为缓存友好而提升了。只有当你要用对偶变量法或重要性抽样等方差缩减技术时才有必要保留完整的路径信息。7. 我自己的工程体会rBergomi模型的C实现真正难的地方不在C语法而在两处一是把分数阶积分和弱奇异核的数学性质吃透二是把蒙特卡洛的性能瓶颈分析清楚。数学理解不到位写出来的代码看起来在模拟rBergomi实际统计性质全是错的性能不做算法层面的优化换什么语言都只能能用但难用。如果你准备在量化项目里落地rBergomi我建议按这个顺序动手先实现朴素卷积版把验证三件套跑通再考虑指数核近似或FFT加速最后做OpenMP并行和编译优化。每一步都有明确的验证手段不会出现改完性能上去了但不知道结果对不对的尴尬局面。本文还有配套的精品资源点击获取