企业官网建设流程全解析

简介：这是一份面向量化金融与期权定价方向的C++工程实现，聚焦rBergomi粗糙波动率模型下的欧式期权定价。定价核心采用蒙特卡洛模拟，分数布朗运动（fBm）由混合方案生成，当前版本实现了恒定前向方差假设，适合具备C++14基础、希望研究粗糙波动率建模与蒙特卡洛定价的开发者与研究者参考。资源包共185个文件，约62.37MB，以cpp与h源码为主体，辅以R脚本、hpp头文件、ipynb实验笔记、tex与pdf文档，以及Makevars、Namespace等构建配置，整体结构覆盖源码、测试与文档多个层面。编译依赖fftw3库与OpenMP支持，需C++14环境，在g++ 5.x及以上版本可正常编译运行。已有175人学习下载。读者可从中获得完整的模型实现代码、蒙特卡洛定价流程、fBm混合生成方案及跨平台编译排错思路，便于在此基础上扩展前向方差结构或对比其他定价方法。

1. rBergomi 的 C++ 实现：从蒙特卡洛定价到路径模拟的工程落地

如果你做过随机波动率模型的衍生品定价，大概率绕不开两个名字：Heston 和 Bergomi。Heston 的好处是有半解析解，校准快，但它的波动率是 CIR 过程，隐含波动率期限结构偏平，短端 skew 拟合经常翻车。rBergomi（rough Bergomi）换了个思路，用分数布朗运动驱动方差过程，Hurst 指数 H 通常取 0.1 左右，能同时把短端 skew 和长端期限结构压得比较准。代价是：没有解析解，只能靠蒙特卡洛，而且路径模拟的协方差结构比标准布朗运动复杂得多。这篇笔记讲的是怎么用 C++ 把 rBergomi 的路径模拟和期权定价跑通，适合已经写过蒙特卡洛定价、想上手 rough volatility 的量化开发。核心问题就三个：方差过程怎么离散、协方差怎么构造、C++ 里怎么把性能压到可接受范围。

2. rBergomi 模型拆解：方差过程、Hurst 指数与定价公式

2.1 模型定义与参数含义

rBergomi 的标准形式在风险中性测度下写成两行：

S_t = S_0 * exp( ∫_0^t sqrt(v_u) dW_u - 1/2 ∫_0^t v_u du ) v_t = ξ_0(t) * exp( η * W_t^H - 1/2 * η^2 * t^{2H} )

其中 W 和 W^H 是两个相关的布朗运动，相关系数 ρ 通常取负值（股票和波动率负相关，这是杠杆效应的来源）。ξ_0(t) 是初始方差曲线，一般从市场隐含波动率校准得到，工程上常用幂函数形式 ξ_0(t) = ξ_0 * t^{2H} 或者直接插值市场 variance swap 报价。η 是波动率的波动率（vol-of-vol），H 是 Hurst 指数，控制路径的粗糙程度。

这里最容易搞混的是 H 的取值。标准布朗运动 H=0.5，路径是半鞅；rBergomi 里 H<0.5，路径不是半鞅，这也是它被称为“rough”的原因。H 越小，短端 skew 越陡，但路径模拟的协方差矩阵条件数越差，数值上越难处理。实操中 H 一般取 0.05 到 0.2，0.1 是最常见的起点。

2.2 为什么不能直接用欧拉离散

很多人第一反应是把 v_t 的指数形式直接欧拉离散：

// 错误示范：直接对 v_t 做欧拉离散 v[i+1] = v[i] + eta * v[i] * dW_H[i] - 0.5 * eta * eta * v[i] * dt;

这样写有两个问题。第一，v_t 的解析形式里 W_t^H 是分数布朗运动，它的增量不是独立的，协方差结构是 E[W_t^H W_s^H] = 1/2 (t^{2H} + s^{2H} - |t-s|^{2H})，直接当独立增量处理会丢掉长记忆性，短端 skew 会明显偏平。第二，指数形式在 η 较大时数值不稳定，容易爆掉。

正确做法是先模拟 W^H 的路径，再代入 v_t 的解析表达式。这样 v_t 永远为正，不需要额外的截断处理。

2.3 协方差矩阵的构造与 Cholesky 分解

模拟 W^H 的标准方法是协方差矩阵法。给定时间网格 0 = t_0 < t_1 < ... < t_n = T，W^H 在各点上的联合分布是均值为零的高斯分布，协方差矩阵 Σ 的元素为：

Σ_ij = 1/2 ( t_i^{2H} + t_j^{2H} - |t_i - t_j|^{2H} )

然后对 Σ 做 Cholesky 分解 Σ = L L^T，用独立标准正态向量 Z 生成路径：W^H = L Z。

// 构造协方差矩阵并做 Cholesky 分解 #include <Eigen/Dense> Eigen::MatrixXd buildCovMatrix(const std::vector<double>& t, double H) { int n = t.size(); Eigen::MatrixXd Sigma(n, n); for (int i = 0; i < n; ++i) { for (int j = 0; j <= i; ++j) { double ti = t[i], tj = t[j]; double val = 0.5 * (std::pow(ti, 2*H) + std::pow(tj, 2*H) - std::pow(std::abs(ti - tj), 2*H)); Sigma(i, j) = val; Sigma(j, i) = val; } } // 加微小对角扰动，防止数值上非正定 Sigma += 1e-10 * Eigen::MatrixXd::Identity(n, n); return Sigma; } Eigen::MatrixXd chol = buildCovMatrix(t, H).llt().matrixL();

这里有几个参数要注意。时间网格 t 的步长不一定要均匀，但均匀网格下 Cholesky 分解的复杂度是 O(n^3)，n=1000 时大概几秒，n=5000 就上分钟了。实操中常用 hybrid scheme 或者 FFT 方法来加速，后面第 5 章会讲。1e-10 的对角扰动是血泪经验，H 小于 0.1 时协方差矩阵条件数能到 1e12 以上，不加扰动 Cholesky 直接报非正定。

2.4 股票价格路径的联合模拟

有了 W^H 的路径，v_t 就确定了。接下来模拟股票价格 S_t，需要处理 W 和 W^H 的相关性。标准做法是把 W 分解成：

W = ρ W^H + sqrt(1 - ρ^2) W_perp

其中 W_perp 是和 W^H 独立的布朗运动。这样每一步的增量 dW 就可以用 dW^H 和独立的 dW_perp 组合出来。

// 联合模拟 v_t 和 S_t double rho = -0.9; // 相关系数 double eta = 1.5; // vol-of-vol double xi0 = 0.04; // 初始方差 double S0 = 100.0; Eigen::VectorXd Z = Eigen::VectorXd::Random(n); // 标准正态 Eigen::VectorXd WH = chol * Z; // 分数布朗运动路径 Eigen::VectorXd v(n), S(n); S[0] = S0; for (int i = 0; i < n; ++i) { double t_i = t[i]; v[i] = xi0 * std::exp(eta * WH[i] - 0.5 * eta * eta * std::pow(t_i, 2*H)); } // 模拟 S 路径，用对数欧拉格式 Eigen::VectorXd Z_perp = Eigen::VectorXd::Random(n); double logS = std::log(S0); for (int i = 1; i < n; ++i) { double dt = t[i] - t[i-1]; double dWH = WH[i] - WH[i-1]; double dW_perp = std::sqrt(dt) * Z_perp[i]; double dW = rho * dWH + std::sqrt(1 - rho*rho) * dW_perp; logS += -0.5 * v[i-1] * dt + std::sqrt(v[i-1]) * dW; S[i] = std::exp(logS); }

对数欧拉格式比直接欧拉稳定，因为 log S 的漂移项和扩散项都是线性的。注意 v[i-1] 用的是上一步的方差，这是欧拉格式的标准做法，误差是 O(dt)。如果想更准，可以用 Milstein 或者 predictor-corrector，但 rBergomi 的主要误差来源是 W^H 的协方差近似，不是 S 的离散误差，所以对数欧拉够用。

3. C++ 工程实现：从单路径到批量蒙特卡洛的完整链路

3.1 项目结构与依赖选择

一个能跑的 rBergomi C++ 实现，目录结构我一般这么组织：

rbergomi/ ├── include/ │ ├── covariance.h // 协方差矩阵构造与 Cholesky │ ├── path_simulator.h // 路径模拟 │ └── pricer.h // 期权定价 ├── src/ │ ├── covariance.cpp │ ├── path_simulator.cpp │ └── pricer.cpp ├── tests/ │ └── test_black_scholes.cpp // 用 BS 验证 MC 框架 └── CMakeLists.txt

依赖就两个：Eigen 做线性代数，OpenMP 做并行。不建议引入太重的框架，rBergomi 的计算核心就是矩阵乘法和随机数生成，自己写反而可控。CMakeLists 里开 -O3 -march=native，Eigen 开 -DEIGEN_USE_MKL_ALL 如果有 MKL 的话。

3.2 协方差矩阵的缓存与复用

蒙特卡洛定价要跑几万到几十万条路径，但协方差矩阵的 Cholesky 分解只跟时间网格和 H 有关，跟路径无关。所以一定要把 chol 矩阵缓存起来，不要每条路径都重新分解。

class PathSimulator { public: PathSimulator(const std::vector<double>& t, double H) : t_(t), H_(H) { Eigen::MatrixXd Sigma = buildCovMatrix(t, H); chol_ = Sigma.llt().matrixL(); } void simulate(double eta, double rho, double xi0, double S0, Eigen::VectorXd& S_out, Eigen::VectorXd& v_out) { int n = t_.size(); Eigen::VectorXd Z = Eigen::VectorXd::Random(n); Eigen::VectorXd WH = chol_ * Z; // ... 后续模拟 } private: std::vector<double> t_; double H_; Eigen::MatrixXd chol_; };

这个缓存策略在 n=500、路径数 10 万的时候，能把总时间从几小时压到几分钟。Cholesky 分解本身是 O(n^3)，但只做一次；每条路径的矩阵向量乘是 O(n^2)，这是主要开销。

3.3 随机数生成器的选择与种子管理

Eigen 的 Random() 用的是 std::rand，质量和速度都一般。实操中我一般换成 PCG64 或者 Mersenne Twister，并且每条路径用独立的种子流，方便复现和并行。

#include <random> class RNG { public: RNG(uint64_t seed) : gen_(seed) {} double normal() { return dist_(gen_); } private: std::mt19937_64 gen_; std::normal_distribution<double> dist_{0.0, 1.0}; }; // 并行时每条线程独立 RNG #pragma omp parallel { int tid = omp_get_thread_num(); RNG rng(12345 + tid * 1000000); #pragma omp for for (int i = 0; i < num_paths; ++i) { // 用 rng 生成路径 } }

种子管理有个坑：如果所有线程用同一个种子，生成的路径会完全一样，蒙特卡洛误差估计直接失效。每条线程的种子要错开至少 1e6，避免 mt19937_64 的状态空间重叠。

3.4 期权定价的蒙特卡洛循环与方差缩减

定价欧式看涨期权，payoff 是 max(S_T - K, 0)，MC 估计是 payoff 的样本均值。

double priceCall(const PathSimulator& sim, double K, double r, double T, int num_paths) { double sum = 0.0, sum_sq = 0.0; #pragma omp parallel for reduction(+:sum, sum_sq) for (int i = 0; i < num_paths; ++i) { Eigen::VectorXd S, v; sim.simulate(eta, rho, xi0, S0, S, v); double payoff = std::max(S[S.size()-1] - K, 0.0); sum += payoff; sum_sq += payoff * payoff; } double mean = sum / num_paths; double var = sum_sq / num_paths - mean * mean; double stderr = std::sqrt(var / num_paths); return std::exp(-r * T) * mean; }

方差缩减方面，对偶变量法（antithetic variates）最简单有效：每条路径用 Z 和 -Z 各跑一次，payoff 取平均。rBergomi 里对偶变量能降 30% 到 50% 的方差，因为 payoff 对 Z 的依赖近似单调。控制变量法也可以用，用 Black-Scholes 价格作为控制变量，但 rBergomi 的隐含波动率跟 BS 差得远，控制变量效果一般。

4. 避坑与排查：rBergomi 实现里最容易翻车的五个点

4.1 协方差矩阵非正定导致 Cholesky 失败

现象：H 小于 0.1 时，Eigen 的 llt() 返回 NaN 或者直接断言失败。

原因：分数布朗运动的协方差矩阵在 H 很小时条件数极大，浮点误差累积导致矩阵在数值上非正定。

解决：加对角扰动 Sigma += 1e-10 * I，或者用 LDLT 分解代替 LLT。如果还不行，把 H 的下限卡在 0.05，再小的话建议换 hybrid scheme。

4.2 路径数够但标准误不降

现象：跑 10 万条路径，标准误还是 0.1 以上，置信区间宽得没法用。

原因：大概率是随机数种子没管好，多条路径用了相同的随机数流，样本之间不独立。

解决：检查每条线程的 RNG 种子是否独立，用 omp_get_thread_num() 错开。另外确认 Eigen 的 Random() 没有被多线程同时调用，它不是线程安全的。

4.3 短端 skew 拟合偏平

现象：校准出来的 H 接近 0.5，模型退化成了普通随机波动率，短端 skew 明显不够。

原因：时间网格太粗，或者协方差矩阵用了近似（比如把 W^H 当独立增量处理）。

解决：短端（一周到一个月）的时间步长要细，至少每天一个点。协方差矩阵必须用精确的分数布朗运动协方差，不能用近似。

4.4 内存爆掉

现象：n=5000、路径数 10 万时，程序被 OOM killer 干掉。

原因：chol 矩阵是 n x n 的稠密矩阵，n=5000 时单精度就要 100MB，双精度 200MB。如果每条路径还存完整的 S 和 v 向量，内存直接爆炸。

解决：chol 矩阵用 float 存（精度损失可接受），路径模拟用原地更新，不存完整路径，只存最后一步的 S_T。定价只需要终端值，中间路径没必要留。

4.5 并行加速比不线性

现象：单线程 10 分钟，开 8 线程只降到 3 分钟，加速比不到 3。

原因：chol 矩阵的矩阵向量乘是内存带宽瓶颈，不是计算瓶颈。多线程同时读 chol 矩阵，内存带宽被打满。

解决：把 chol 矩阵转成行优先存储，提高缓存命中率。或者用分块矩阵乘法，每次只加载一部分 chol 到缓存。再不行就上 GPU，但那是另一个话题了。

5. 进阶技巧：Hybrid Scheme 与 GPU 加速的取舍

协方差矩阵法的瓶颈在 O(n^3) 的 Cholesky 分解和 O(n^2) 的矩阵向量乘。n=1000 时还能忍，n=5000 以上就必须换方法。Hybrid scheme 是现在比较主流的替代方案，核心思路是把分数布朗运动的协方差核近似成几项指数函数的和，然后用 Ornstein-Uhlenbeck 过程递推模拟，复杂度降到 O(n)。

具体做法是把 t^{2H} 的协方差核用 Laplace 变换展开，取前几项做主成分近似。实操中取 5 到 10 项就能把误差控制在 1e-4 以内。这样每条路径的模拟从 O(n^2) 降到 O(n * k)，k 是展开项数，n=5000 时能快 100 倍以上。

// Hybrid scheme 的递推核心 // 把 W^H 近似成 k 个 OU 过程的线性组合 for (int i = 1; i < n; ++i) { double dt = t[i] - t[i-1]; double dWH = 0.0; for (int j = 0; j < k; ++j) { // OU 过程递推 ou[j] = ou[j] * std::exp(-lambda[j] * dt) + std::sqrt((1 - std::exp(-2*lambda[j]*dt)) / (2*lambda[j])) * std::sqrt(dt) * Z[j]; dWH += w[j] * ou[j]; } WH[i] = WH[i-1] + dWH; }

lambda 和 w 是近似系数，需要预先算好。这套方法的精度在 H=0.1、k=8 时，和精确协方差矩阵法的价格差异在 0.1% 以内，但速度快两个数量级。

GPU 加速是另一个方向。chol 矩阵乘法和路径模拟都是天然并行的，用 CUDA 或者 SYCL 能把 10 万条路径的时间从分钟级压到秒级。但 GPU 的显存有限，n=5000 时 chol 矩阵就要 200MB，加上路径数据，显存很容易不够。实操中 GPU 方案适合 n 小于 2000 的场景，n 大了还是得上 hybrid scheme。

我自己的习惯是：n 小于 500 用协方差矩阵法，代码简单不容易错；n 在 500 到 2000 之间用 hybrid scheme，速度和精度平衡得最好；n 大于 2000 就上 hybrid scheme 加 GPU，但这时候校准的复杂度已经很高了，要慎重评估值不值得。rBergomi 的坑主要在数值稳定性上，模型本身不复杂，但 H 小的时候协方差矩阵的条件数能到 1e12，双精度都不一定扛得住。我一般会在代码里加个断言，检查 chol 矩阵的对角元是否都为正，一旦出现负数就自动加扰动重试。这个习惯帮我省了很多调试时间。希望帮到你。

本文还有配套的精品资源，点击获取