拓冰网站定制 
Lévy 过程莱维过程这个词我第一次真切记住它是从一次做市亏钱的教训里来的。那阵子我盯的期权组合短期限隐含波动率曲线两端翘得厉害用连续扩散那一套怎么调参都对不上来回折腾两周才认账财报日、突发公告那种一夜跳空 8% 的行情模型里根本就没有位置参数再掰也是硬掰。后来朋友甩给我一篇讲跳扩散的论文我才算正式走进 Lévy 过程这扇门。它做的事情说穿了很朴素——把漂移、连续扩散、跳跃这三件事装进同一个数学框架用一个三元组 (b, σ², ν) 讲清楚一个随机过程到底是怎么走的。做量化定价的、做风险计量的、写保险精算模型的、研究排队系统的甚至做粒子输运和反常扩散研究的最后都会撞上它。这篇文章我不打算从测度论开始堆定义而是按我踩坑的顺序把定义、三元组、从属过程、参数标定、定价落地和模拟误差这几块拆开讲能直接抄的代码我都贴上。1. 为什么连续无跳跃的假设会翻车1.1 布朗运动的两个隐藏代价布朗运动是我们在入门阶段最熟悉的连续时间模型它的路径连续、增量独立、增量服从正态分布。这三条看起来很温和实际上藏着两个致命代价只要把它拿去拟合真实的价格序列问题立刻暴露。第一个代价是路径连续性。连续路径意味着在极小的时间尺度上位移量级只有 O(√Δt)。换句话说无论你把采样频率提到多高单步位移都会随步长缩小而平稳地缩小永远不可能出现隔夜跳 8%这种断崖。可现实里这种事每年都要发生几次。更麻烦的是连续性假设下你没法用跳跃强度这个概念去描述事件到达的密集程度只能被动地靠提高波动率去吸收结果就是平时波动率被高估、关键时刻又被低估。第二个代价来自正态分布的尾部。正态分布的尾部衰减速度是指数平方级的5 倍标准差的事件概率被压到小数点后好多位。但真实数据里日频收益的超额峰度往往在 3 到 10 之间也就是同样的标准差下极端事件的频率比正态预测高出好几个数量级。这两件事叠加起来还会导致第三个副作用方差随时间线性增长年化波动率变成常数任意时间尺度上的分布形状完全一样。可实际数据里期限越短的分布越尖、偏斜越重这跟形状不随时间尺度变化的性质正好相反。1.2 独立增量与平稳增量定义只有四行Lévy 过程的定义其实短得惊人。一个实值随机过程 X {X_t}{t≥0}如果满足下面四条就叫 Lévy 过程X_0 0 几乎必然成立任意 0 ≤ t_1 … t_n各个增量 X{t_i} - X_{t_{i-1}} 相互独立增量平稳也就是 X_{ts} - X_s 的分布只依赖 t 而不依赖 s以及随机连续性对任意 ε 0当 s 趋近 t 时 P(|X_s - X_t| ε) 趋近 0。严格来说第四条在数学上可以由前三条推出来这是 Lévy 本人的一个定理但把它直接写进定义里省事也符合直觉。前两条的含义值得多说一句独立增量意味着过程没有记忆未来怎么走跟过去走成什么样无关平稳增量意味着时间齐次随机机制的规则不随时间改变。这两条合起来就保证了 X_t 的分布一定是单步分布 Q 的 t 次卷积反过来也就推出 Q 必须是无穷可分分布。Lévy 过程和无穷可分分布之间是一一对应的这也是它在极限定理里能扮演高斯分布之于中心极限定理那个角色的根本原因。平稳性这一条在实际建模里需要留个心眼。波动率有日内周期开盘半小时和收盘半小时完全不是一个物种也有区制切换平静期和危机期的机制差别巨大。严格讲日内不同时段的增量分布并不相同直接拿裸数据去拟合会得到一堆自相矛盾的参数。我的一般做法是先做时间去季节化把每个时间点的收益按该时段的历史波动率标准化或者干脆只在同质窗口内做估计。1.3 两个极端之间的一整条光谱布朗运动和泊松过程可以看成两个极端。前者只有无穷多的小位移路径连续方差有限后者只有正跳、幅度恒为定值完全没有连续成分。真正的 Lévy 过程世界是这两个极端之间的一整条光谱而这条光谱上散落着大量在工程上非常好用的模型。过程三元组 (b, σ², ν) 的关键特征路径特点典型场景布朗运动(μ, σ², 0)连续正态增量扩散近似、基准模型泊松过程(0, 0, λδ_1)纯跳跳幅固定事故计数、订单到达复合泊松(0, 0, λF)跳幅服从分布 F单日索赔总额、批量订单Gamma 过程σ² 0ν 集中于正半轴单调不减只有正跳累积降雨量、信用损失累积方差 Gammaσ² 0ν 为两侧指数衰减的 1/|x| 型无穷多小跳 有限大跳日频对数收益α-稳定过程σ² 0ν(dx) ∝ |x|^{-1-α}跳跃测度在零点爆炸方差无穷极端尾部建模Merton 跳扩散(b, σ², λ·N(μ_J, σ_J²))连续扩散 高斯跳权益期权定价CGMYσ² 0ν 为指数阻尼的幂律跳跃活跃度可精细调节高频精细校准这张表里我最想强调的是 α-稳定过程那一行。它是把小跳跃强度无穷大这件事推到极致的结果ν 在零点附近按 |x|^{-1-α} 爆炸导致过程连二阶矩都不存在。这听起来很吓人但恰恰就是很多重尾数据的真实写照。而广义中心极限定理告诉我们独立同分布随机变量和在适当归一化下的极限要么是高斯分布要么是 α-稳定分布——换句话说只要你的数据生成机制允许无穷方差最后收敛到的就是稳定过程。Lévy 过程之所以值得认真学这个定理是最大的理论靠山。2. 三元组把 Lévy-Khintchine 公式读明白2.1 特征函数是最好用的那把尺子Lévy 过程有一个极其漂亮的性质它的特征函数有封闭形式。对任意 t ≥ 0 和 u ∈ RE[e^{iuX_t}] e^{tψ(u)}其中 ψ 叫 Lévy 指数形式是ψ(u) ibu - σ²u²/2 ∫_{R{0}} (e^{iux} - 1 - iux·1_{|x|1}) ν(dx)为什么特征函数能写成 t 次方的指数形式这完全是独立平稳增量的直接后果X_t 是 t 个独立同分布增量的和特征函数自然就是单步特征函数的 t 次方。这个公式最实用的地方在于它告诉你只要知道 (b, σ², ν) 这三个东西过程的全部有限维分布就确定了一整套定价、模拟、估计的工作都可以围绕这三个量展开。那为什么不直接用概率密度而要用特征函数原因很现实绝大多数 Lévy 过程的密度函数没有初等表达式。方差 Gamma 的密度要用第二类贝塞尔函数写CGMY 连贝塞尔函数都写不出来只能用无穷级数或者数值反演。但它们的特征函数几乎总是初等函数形式干净、计算便宜、求导方便。所以整条工程链条——参数估计、期权定价、风险度量、随机数生成——都是围绕特征函数转的。我刚开始做这块的时候总想着先把密度搞出来再说结果在特殊函数里绕了大半个月后来老老实实改走特征函数通道效率提升了不止一个量级。2.2 跳跃测度 ν 和那个不得不加的减号跳跃测度 ν 的读法是ν(A) 表示单位时间内跳跃幅度落进集合 A 的期望次数也就是 E[N((0,t]×A)] t·ν(A)。它不是一个概率测度而是一个强度测度所以 ν(R{0}) 完全可以是无穷大。事实上对于方差 Gamma 和 CGMY 这类模型零点附近的小跳跃强度就是无穷的一秒钟里发生无限多次极其微小的跳跃。但有个约束必须满足∫ min(1, x²) ν(dx) ∞。这条把跳跃测度在无穷远处的增长和在零点附近的增长同时管住了。现在解释那个看起来很别扭的减号 -iux·1_{|x|1}。道理在于幅度大于 1 的大跳跃其总强度 ν({|x|≥1}) 是有限的这部分可以直接处理但幅度小于 1 的小跳跃强度可能是无穷的如果直接写 ∫ x N(dt,dx) 这一项它的期望可能根本不存在或者发散。减掉 iux 相当于对每一层小跳跃先做中心化让积分收敛。生活化的类比是这样的你每天要交无数笔小额手续费单笔金额都很小但笔数无穷多直接把所有笔加起来是发散的必须先对每一笔减去它的期望值这个总偏差才是收敛的、有意义的东西。这个减号的副作用必须记住漂移参数 b 依赖于截断函数的选择。如果你把 1_{|x|1} 换成 1_{|x|c}b 的数值就会变而 σ² 和 ν 保持不变。所以写论文或者写接口文档的时候一定要标注清楚用的是哪种截断约定否则别人照着复现一定对不上。我自己就栽过这个坑照着某篇文献的参数往代码里填模拟出来的路径均值怎么都偏查了整整两天最后发现是文献用的截断和我用的不一致b 差了 ∫_{|x|1} x ν(dx) 那一项。顺带给一个能直接用的一致性关系。在 1_{|x|1} 的约定下拆成漂移 扩散 补偿小跳 大跳之后期望满足 E[X_t] t·(b ∫_{|x|≥1} x ν(dx))因为补偿泊松积分是鞅、均值为零。写完模型先拿这个式子做一次数值校验能挡掉不少低级错误。2.3 Lévy-Itô 分解把过程切成四块如果说特征函数是看 Lévy 过程的方式那 Lévy-Itô 分解就是造 Lévy 过程的方式。它把任意 Lévy 过程唯一地拆成四块X_t bt σW_t ∫_0^t∫_{|x|1} x·Ñ(ds,dx) ∫_0^t∫_{|x|≥1} x·N(ds,dx)四块分别是确定性的线性漂移、标准的连续布朗部分、无穷多次被补偿的小跳跃、有限次不做补偿的大跳跃其中 Ñ N - ν(dx)ds 是补偿泊松随机测度。这四块之间相互独立这个独立性在模拟和蒙特卡洛里非常关键——意味着你可以分开采样再叠加不会引入相关性误差。这个分解最直接的价值是给了一切实操的起点。想模拟分别处理这四项。想推导定价方程把生成元写出来Lf(x) b·f(x) (σ²/2)·f(x) ∫ [f(xy) - f(x) - y·f(x)·1_{|y|1}] ν(dy)有了生成元套用 Feynman-Kac 就能得到跳扩散框架下的定价方程。以权益期权为例价格函数 V(t,S) 满足的是一个偏积分微分方程PIDE比 Black-Scholes 方程多的就是那个积分项∂V/∂t (r - λκ)S·∂V/∂S (σ²S²/2)·∂²V/∂S² ∫ [V(t, S·e^y) - V(t,S)] ν(dy) - rV 0这个积分项是求解时的头号麻烦。有限差分本来能生成稀疏的三对角矩阵加上积分项之后矩阵立刻变稠密直接求逆的代价从 O(N) 涨到 O(N³)。工程上的标准做法是把积分项当成卷积用 FFT 以 O(N log N) 算出来然后配合隐式处理微分算子、显式处理积分算子的 IMEX 时间离散格式。这套组合拳我在做美式期权 PDE 求解的时候用得最多实测稳定性比全隐式好速度比全显式快。3. 从属过程与时间变换最实用的一个技巧3.1 用随机时钟换掉日历时间从属过程subordinator是一类特殊的 Lévy 过程它单调不减只有正跳跳跃测度完全集中在正半轴并且满足比一般情形更强的约束 ∫ min(1, x) ν(dx) ∞。它的拉普拉斯变换有简洁形式 E[e^{-λS_t}] e^{-tΨ(λ)}其中 Ψ(λ) γλ ∫_0^∞ (1 - e^{-λx}) ν(dx) 叫拉普拉斯指数。从属过程最有价值的用法是做时间变换。取一个标准的布朗运动 W再取一个独立的从属过程 S构造 X_t μ·S_t σ·W_{S_t}。这个构造的直觉非常贴切布朗运动本身描述的是每一步都是独立小扰动的扩散但真实的成交和信息到达并不匀速——开盘、数据发布、突发事件的时候市场活跃度暴涨午间流动性枯竭的时候几乎不动。如果用一只有随机指针的时钟去给布朗运动计时活跃时段指针走得快、平静时段走得慢那么最终观察到的日历时间里的价格过程自然就带上了尖峰厚尾。这个物理含义在工程上有个更接地气的版本拿交易笔数或者成交量当业务时钟去替换日历时间很多统计性质会立刻变好。这不是玄学而是因为从属过程本质上就是把事件到达速率随机化而成交量恰好是事件到达速率的一个可观测量。我在处理高频数据的时候用成交量时钟重采样之后的收益率正态性检验的 p 值经常能从 10 的负若干次方提升到 0.05 以上。3.2 方差 Gamma 与 NIG把布朗运动拖慢的两个经典做法最简单的时间变换是拿 Gamma 过程当 S。Gamma 过程 S_t 服从形状参数 t/ν、尺度参数 ν 的 Gamma 分布因此 E[S_t] t、Var[S_t] ν·t。把它代入得到 X_t θ·S_t σ·W_{S_t}这就是方差 Gamma 过程它的特征函数是φ(u) (1 - iθνu σ²νu²/2)^{-t/ν}三个参数的分工非常清楚θ 控制偏度正负都能做ν 控制尾部厚度和超额峰度σ 是整体的扩散尺度。特别要注意 ν 趋近 0 的极限此时 Gamma 时钟退化成确定性的钟方差 Gamma 过程收敛到布朗运动。这个极限正好解释了后面要讲的参数识别噩梦。换成逆高斯从属过程得到的就是 NIG 过程正态逆高斯它和方差 Gamma 是同一个正态方差混合家族的两个成员区别只在于时钟的分布。两者都是纯跳过程没有连续成分但小跳跃强度无穷大所以在有限的分辨率下看起来几乎是连续的。单位时间内的累积量可以拿来做参数含义的定量判断参数表达式说明κ₁θ均值直接给出漂移方向κ₂σ² θ²ν方差扩散与跳动的叠加κ₃3θσ²ν 2θ³ν²三阶累积量偏度的来源κ₄3σ⁴ν 12θ²σ²ν² 3θ⁴ν³四阶累积量峰度的来源超额峰度等于 κ₄ / κ₂²你能从这张表直接看出想让数据显得更尖要么加大 ν要么加大 |θ|。这一步的推演在实际调参时非常省时间比盲试参数高效得多。3.3 方差 Gamma 的模拟代码直接按时间步采样 Gamma 增量是最省事的做法。注意增量独立性让你可以逐段生成然后累加。import numpy as np def simulate_vg(n_steps, T, theta, sigma, nu, n_paths1, seedNone): 模拟方差 Gamma 过程 X_t theta * S_t sigma * W_{S_t} S_t 是 Gamma 从属过程S_t ~ Gamma(shapet/nu, scalenu) 参数: n_steps: 时间步数 T: 总期限 theta: 偏度参数 sigma: 扩散尺度 nu: 尾部厚度参数 (nu - 0 退化为布朗运动) 返回: shape (n_paths, n_steps 1) 的路径数组首列为 0 rng np.random.default_rng(seed) dt T / n_steps if dt / nu 0.05: raise ValueError(dt/nu 过小Gamma 采样会大量下溢请增大步长或调小 nu) shape dt / nu # Gamma 增量的形状参数 scale nu # Gamma 增量的尺度参数 dS rng.gamma(shape, scale, size(n_paths, n_steps)) # 条件在 Gamma 时钟上布朗部分增量的方差是 sigma^2 * dS dW rng.standard_normal((n_paths, n_steps)) * np.sqrt(dS) dX theta * dS sigma * dW X np.cumsum(dX, axis1) return np.concatenate([np.zeros((n_paths, 1)), X], axis1)有一个坑必须提醒当 dt/ν 很小比如 ν 0.1、日频数据 dt 1/250商只有 0.04时Gamma 采样会产生大量接近零的值浮点数下溢之后再开方、再乘正态会引入明显的数值偏差你会看到路径方差不匹配。我一般把商控制在 0.05 以上做不到就改用步长自适应或者用 Gamma 桥在已知终点值的条件下做桥式插值来生成路径。实测下来Gamma 桥在 ν 小于 0.05 的场景里几乎是唯一可靠的选择代价是实现复杂度明显上升。4. 参数标定从数据走到模型4.1 用经验特征函数做拟合既然特征函数是初等形式最顺理成章的估计法就是最小化理论特征函数和经验特征函数之间的距离。这个方法对纯跳过程特别友好因为它们的密度要么难算要么根本不存在。import numpy as np from scipy.optimize import minimize def vg_cf(u, theta, sigma, nu): 方差 Gamma 过程单位时间的特征函数 phi(u) E[exp(i u X_1)] return (1 - 1j * theta * nu * u 0.5 * sigma**2 * nu * u**2) ** (-1.0 / nu) def fit_vg_by_cf(log_returns, u_max50.0, n_u80): 用经验特征函数与理论特征函数的加权距离估计 VG 参数。 log_returns: 单位时间如日频的对数收益序列 r np.asarray(log_returns, dtypefloat) r r - r.mean() # 先去掉均值逐渐比较形状 u_grid np.linspace(0.2, u_max, n_u) def ecf(u): return np.mean(np.exp(1j * u * r)) ecf_vals np.array([ecf(u) for u in u_grid]) # 权重相位信息集中在 |ECF| 较大处用 n/|ECF|^2 做最优加权 n len(r) w n / np.maximum(np.abs(ecf_vals) ** 2, 1e-8) def loss(p): theta, sigma, nu p if sigma 1e-6 or nu 1e-4: return 1e10 diff np.array([vg_cf(u, theta, sigma, nu) - e for u, e in zip(u_grid, ecf_vals)]) return float(np.sum(w * np.abs(diff) ** 2) / n_u) p0 np.array([0.0, r.std(), 0.2]) best None for nu0 in (0.05, 0.2, 0.5, 1.0): # 多初值避开局部极小 res minimize(loss, [p0[0], p0[1], nu0], methodNelder-Mead, options{xatol: 1e-8, fatol: 1e-10, maxiter: 4000}) if best is None or res.fun best.fun: best res theta, sigma, nu best.x # 还原均值到 theta 上单位时间的漂移 return theta r.mean(), sigma, nu, best.fun这段代码有三个地方值得展开讲。第一先减掉均值再拟合形状参数是因为均值和偏度参数 θ 存在直接耦合一起拟合会让目标函数出现长长的浅谷收敛极慢先固定均值、只拟合尺度与尾部最后再把均值加回 θ稳定性好很多。第二权重用 n/|φ̂(u)|² 而不是均匀权重这是有理论依据的渐近最优权重原因是高频部分的经验特征函数被噪声主导其估计偏差量级是 O(1/√n)不降权就会把噪声当信号拟合进去。第三多初值是必须的因为 CGMY 和 VG 这类模型的目标函数经常是多峰的单初值的 Nelder-Mead 在 ν 方向容易卡住不动。u 网格的上限选择也有讲究。日频对数收益的尺度在 0.01 左右主要相位信息集中在 u ≤ 50 这个区间u 再往上理论特征函数的模已经衰减到千分之一以下经验估计的噪声完全盖过信号。我一般的做法是画一条理论 ECF 与经验 ECF 的模随 u 变化的曲线找一个交叉点把 u_max 定在交叉点之前。4.2 累积量匹配快手但脆弱的做法如果你手上只有样本的前几阶矩用累积量匹配能很快得到一组粗略参数特别适合作为上面数值优化的初值。以方差 Gamma 为例记单位时间的样本累积量为 c₁、c₂、c₃、c₄则 θ c₁ν 满足一元二次方程θ³ν² - 3θc₂ν c₃ 0取哪个根要配合约束条件判断必须同时满足 σ² c₂ - θ²ν 0以及判别式 9c₂² ≥ 4θc₃。选好根之后 σ² 直接从 c₂ 里减出来c₄ 留着做一致性检验——如果代入 c₄ 的残差大得离谱说明数据里还有 VG 描述不了的结构比如波动率聚集这时候硬拟合就是自欺欺人。这套方法的短板在于它对高阶矩极其敏感。c₃ 和 c₄ 的样本估计标准误很大几个异常值就能把参数带飞。所以我从来不把它当最终答案只当数值优化的起跳点。真正靠谱的还是 4.1 节那种用完整特征函数曲线做最小化的方式。4.3 从数据里把跳跃抠出来模型参数估计完之后还有一个反向验证的问题数据里到底有没有跳跃这里的经典工具是双幂变差。它的出发点很聪明——已实现方差 RV Σr_i² 在存在跳跃时会包含跳跃变差而双幂变差 BV (π/2)·Σ|r_i||r_{i-1}| 对跳跃稳健只捕捉连续部分的积分方差。两者之差就是跳跃变差。import numpy as np def jump_variation(r, k1): 用双幂变差把已实现方差拆成连续变差和跳跃变差。 r: 高频对数收益已去均值 k: 滞后阶数通常取 1 r np.asarray(r, dtypefloat) absr np.abs(r) RV float(np.sum(r ** 2)) BV float((np.pi / 2.0) * np.sum(absr[k:] * absr[:-k])) JV max(RV - BV, 0.0) ratio JV / RV if RV 0 else 0.0 return {RV: RV, BV: BV, JV: JV, jump_share: ratio}我在用这个指标时有两条血泪经验。第一微观结构噪声会同时抬高 RV 和 BV在采样频率高于几秒的时候噪声项会彻底压过跳跃信号所以必须先做子采样或者噪声校正否则你会检测出一大堆根本不存在的跳跃。第二jump_share 在不同交易日之间波动极大单日的结果没有意义必须看滚动窗口的均值和分布再配合跨资产比较才能判断跳跃是不是一个稳定的模型特征。如果某个资产的 jump_share 长年在 5% 以下那给它配一个纯跳模型多半是过度设计老老实实上跳扩散甚至纯扩散更划算。5. 落到业务里定价、风控与排队5.1 特征函数定价一条公式吃遍所有模型Lévy 过程在衍生品定价里最爽的一点在于不管你选哪个具体模型定价流程是完全一样的因为大家都有特征函数。欧式看涨期权价格可以写成C S₀·Π₁ - K·e^{-rT}·Π₂其中 Π₂ Q(S_T K) 是风险中性测度下的行权概率Π₁ Q^S(S_T K) 是份额测度下的行权概率两者都能用 Gil-Pelaez 反演公式从特征函数直接算出来Π_j 1/2 (1/π)∫₀^∞ Re[e^{-iu·log K}·φ_j(u)/(iu)] du。份额测度的特征函数通过测度变换得到φ^S(u) φ(u - i)/φ(-i)。import numpy as np def call_price_by_cf(S0, K, r, T, cf, u_max200.0, n_u2**15): 用 Gil-Pelaez 反演计算欧式看涨期权价格。 cf: 特征函数 phi(u) E[exp(i u X_T)]X_T log(S_T / S0) 必须已经做过鞅修正满足 cf(-1j) exp(r*T) k np.log(K / S0) u np.linspace(1e-8, u_max, n_u) phi_u cf(u) # 风险中性测度 phi_ui cf(u - 1j) # 份额测度 phi_mi cf(-1j) # 应等于 exp(r*T) integ2 np.real(np.exp(-1j * u * k) * phi_u / (1j * u)) integ1 np.real(np.exp(-1j * u * k) * phi_ui / (1j * u * phi_mi)) Pi1 0.5 np.trapezoid(integ1, u) / np.pi Pi2 0.5 np.trapezoid(integ2, u) / np.pi return S0 * Pi1 - K * np.exp(-r * T) * Pi2鞅修正是这一段最关键也最容易漏的一步。以方差 Gamma 为例模型原始形式给出的对数收益特征函数对应的期望不是无风险利率你必须加一个漂移补偿量 ω满足ω r (1/ν)·ln(1 - θν σ²ν/2)这里的约束 1 - θν σ²ν/2 0 必须成立否则 E[S_T] 发散整个定价框架失效。我第一次做这块的时候忘了加 ω算出来的看涨期权价格比看跌还便宜对着屏幕愣了半天才反应过来。所以接口里最好加一句断言直接检查 cf(-1j) 是不是等于 exp(r*T)事后再检查价格合理性就晚了。u_max 和 n_u 的取值决定精度。深虚值期权对 u_max 极其敏感u_max 太小会把高频信息截掉价格错得离谱。判断方法很简单把 u_max 翻倍再算一遍看价格变化量是否落在你能容忍的范围内比如 0.01%。如果振荡明显就说明被积函数衰减太慢这时候需要上阻尼技巧——把看涨价格乘上 e^{αk}α 是正阻尼参数让被积函数在复平面上衰减得更快再配合 FFT 一次性算出一整条行权价曲线。做市系统的报价表基本都是这么生成的单次 FFT 算出几百个行权价比逐个反演快很多。5.2 保险精算里的破产概率Lévy 过程在保险里的用法甚至比金融领域更直接。经典的 Cramér-Lundberg 模型把保险公司的盈余写成 U_t u ct - Σ_{i1}^{N_t} Y_i其中 u 是初始资本c 是保费率N 是强度 λ 的泊松过程Y_i 是独立同分布的索赔额。这个 U 本身就是个 Lévy 过程——线性漂移加上一个复合泊松跳跃没有扩散项。在这个框架下安全负荷率 ρ 决定了保费率 c (1 ρ)λE[Y]调节系数 R 由方程 λ(M_Y(R) - 1) cR 唯一确定其中 M_Y 是索赔额的矩母函数。Lundberg 不等式给出破产概率的上界 ψ(u) ≤ e^{-Ru}指数衰减的速度直接由 R 控制。推广方向有两个都很自然。一是把复合泊松换成一般的从属过程得到更灵活的 Lévy 保险风险模型能描述索赔频率本身也随机变化的情形二是往盈余过程里加一个布朗扰动项得到扩散扰动风险模型用来描述保费收入的连续波动和投资收益。这两条推广在文献里统称Lévy 风险过程好处是所有分析工具都从特征函数和拉普拉斯指数直接继承过来不需要重新推一遍。有意思的是这里的 Pollaczek-Khinchine 公式和排队论里的稳态等待时间公式本质上是同一个式子。盈余过程的破产概率、M/G/1 队列的稳态工作量、反射 Lévy 过程在零点的局部时三者在数学上是一套东西的三个投影。我个人觉得理清这层关系比背十个公式有用得多。5.3 反射过程与排队反射 Lévy 过程的标准构造是 Q_t X_t L_t其中 L 是 Q 在零点处的局部时负责把过程推回非负半轴。它的离散版本就是 Lindley 递归 W_{n1} max(W_n X_n, 0)几乎每个写过排队模拟的人都用过。连续版本的价值在于稳态工作量的拉普拉斯变换可以用 Lévy 过程的谱理论写成显式形式从而避开长时间的瞬态模拟。这一块我印象最深的是数值稳定性问题。反射过程在接近零点时局部时会迅速累积如果模拟时步长太粗你会看到路径在零附近反复穿越局部时被严重高估稳态分布的右尾被压扁。我的做法是做步长收敛测试把步长减半两次比较稳态分布的 95% 分位数变化超过百分之一就继续细化。这个测试花不了几分钟但能省掉后面大量的返工。6. 模拟误差控制真正花钱的地方6.1 小跳跃截断的偏差从哪来对于跳跃测度在零点附近无界的模型α-稳定、CGMY 的 Y 0、方差 Gamma你不可能真的把每一次微小跳跃都模拟出来。标准做法是设一个截断阈值 ε只模拟 |x| ε 的跳跃把 |x| ≤ ε 的部分近似成扩散加漂移。具体地补偿漂移用小跳跃的一阶矩 b_ε ∫_{|x|≤ε} x ν(dx)扩散系数用二阶矩 σ_ε² ∫_{|x|≤ε} x² ν(dx)。这个替换带来两层误差。一层是矩的不匹配被替换掉的真实小跳跃有非零的三阶及以上累积量而替代的高斯扩散三阶以上累积量为零所以分布的形状偏了。另一层是随机性结构的差异一阶和二阶矩匹配了但跳跃在时间上的聚集特性和扩散完全不同。α-稳定过程下截断偏差随 ε 的收敛速度比较慢通常是 O(ε^{2-α}) 量级这意味着 α 越接近 2收敛越慢需要越小的 ε。工程上的做法很简单把 ε 依次取 1e-3、5e-4、2.5e-4算三次价格看这个序列的模式。如果差值按预期比例缩小说明已经在渐近区间内如果两次之间跳动很大而且方向不定说明还没进渐近区需要更小的 ε 或者检查代码是不是有 bug。我一般把这个测试写成一个 pytest 用例每次改动模拟器都跑一遍防止某次重构悄悄把截断逻辑改坏。6.2 方差缩减三个我常用的招第一个招是条件蒙特卡洛。以方差 Gamma 为例条件在 Gamma 时钟的整条路径 S 上资产价格的对数收益是正态的欧式期权的条件期望可以直接用 Black-Scholes 公式解析算出蒙特卡洛只需要处理从属过程这一低维随机源。这一步能把方差降低一到两个数量级是所有技巧里性价比最高的。对路径依赖的期权只能做到部分条件化但即使只条件化终点附近的时段方差下降也很明显。第二个招是重要性抽样。给跳跃测度做指数倾斜Esscher 变换把稀有的大跳跃拉到分布中心附近再抽样同时给每条路径乘上对应的似然比权重。这个技巧在计算深虚值期权和极端分位数风险度量比如 99.9% 违约阈值的时候几乎是必需的因为裸蒙特卡洛在这种场景下的有效样本量可能小到个位数。有意思的是Esscher 变换既是定价时的测度变换工具又是模拟时的方差缩减工具同一套数学两个用途。第三个招是控制变量。用连续部分把跳跃全部关掉的贴现收益作为控制变量因为它的解析期望已知而且和全模型的收益高度相关。实测下来对短期期权能把方差降低 50% 到 70%实现成本还很低属于白捡的收益。7. 常见问题与排查速查7.1 参数识别为什么会拟合出假的σ这是我在 Lévy 过程建模里遇到频率最高的问题。症状很典型拟合出来的 ν 一路往零跑σ 稳定在样本标准差附近θ 稳定在样本均值附近最终结果退化成布朗运动。原因不是优化器坏了而是数据本身没法区分无穷多极小跳跃和连续扩散。有限样本下两者对特征函数在这段区间的贡献几乎不可分目标函数在 ν 方向是一条长长的平底谷。诊断方法很简单固定 θ 和 σ把 ν 从 0 到 1 扫一遍画轮廓目标函数曲线。如果曲线在 ν 小于某个值之后基本是平的就说明这个参数在这个数据上不可识别。处理手段有三条按我实际用过的优先级排列换成更高频的数据小跳跃在细时间尺度上更突出识别力提升明显把 ν 固定到文献里的典型量级只估计剩下的参数或者上贝叶斯方法用弱信息先验把参数约束在一个合理范围让后验分布而不是点估计来回答问题。最后这条在样本量小的时候特别好用因为后验分布会把这个参数其实不太确定这件事诚实地表现出来而不是给你一个看起来很精确的假数字。7.2 数值反演的三个典型故障第一个故障是深虚值期权价格明显失真原因一般是 u_max 设得太小。判断方法是把 u_max 翻倍重算如果价格变化超过阈值就说明之前截断了有效信息。千万不要凭感觉设 u_max一定要做收敛测试。第二个故障是价格曲线出现上下振荡。它来源于被积函数随 u 衰减太慢离散化时的截断误差在反演积分里形成了类似 Gibbs 现象的振荡。解法是引入 Carr-Madan 的阻尼因子把被积函数乘上一个在虚轴方向增长、在实轴方向快速衰减的函数让尾部行为变好同时把积分网格加密配合 Simpson 权重而不是简单的梯形权重。第三个故障是模型隐含的远期价格对不上市场远期。这几乎总是鞅修正漏了或者做错了。记住那个检查点cf(-1j) 必须等于 e^{rT}。如果分红率、外汇双币种利率这些都在里面就得把对应项补全否则整个定价链路会在最基础的一致性上崩掉。7.3 常见问题速查表现象可能原因我的处理方式ν 一路趋近 0σ ≈ 样本标准差参数不识别小跳跃与扩散不可分换高频数据、固定 ν、或改用贝叶斯方法深虚值价格失真反演积分上限 u_max 太小翻倍 u_max 做收敛测试价格曲线上下振荡被积函数衰减慢、网格太粗Carr-Madan 阻尼因子 加密网格模拟路径方差与理论不符时间步长过粗、Gamma 采样下溢增大步长或改用 Gamma 桥E[S_T] 与 e^{rT} 对不上漏了漂移补偿 ω按 ω r ln(1-θνσ²ν/2)/ν 修正跳跃检测结果异常多微观结构噪声被当成跳跃子采样、噪声校正、提高阈值CGMY 拟合不收敛C 与 Y 强相关固定 Y 做网格搜索再优化其余参数截断阈值 ε 缩小时价格不收敛还没进渐近区或代码有 bug依次减半 ε 观察收敛模式这张表里的每一条我都真实踩过尤其是C 与 Y 强相关那条让一个模型调了两天才发现不是优化器的问题而是参数本身在数据长度不够的时候根本分不开。经验是CGMY 这种四参数模型日频数据至少要有三年以上、并且包含至少一次明显的跳跃事件否则不要指望能把四个参数都稳定估出来。最后分享一个我自己的习惯。每次上线一个新模型之前我一定做三件事拿布朗运动做退化工况测试把跳跃参数推向极限看代码能否平滑退化到 Black-Scholes 解析解拿一条人工构造的已知参数路径做回环测试生成数据、再拟合、看能不能把参数找回来以及拿市场上真实报价做一次反算隐含参数的一致性检查。这三步花掉的时间和事后排查一个生产事故相比便宜得不像话。踩过几次坑之后我现在看到任何拟合结果特别漂亮的模型第一反应都是先怀疑参数识别和数据质量而不是先庆祝。 Lévy过程与跳扩散建模：三元组、特征函数、参数标定与定价落地
Lévy 过程莱维过程这个词我第一次真切记住它是从一次做市亏钱的教训里来的。那阵子我盯的期权组合短期限隐含波动率曲线两端翘得厉害用连续扩散那一套怎么调参都对不上来回折腾两周才认账财报日、突发公告那种一夜跳空 8% 的行情模型里根本就没有位置参数再掰也是硬掰。后来朋友甩给我一篇讲跳扩散的论文我才算正式走进 Lévy 过程这扇门。它做的事情说穿了很朴素——把漂移、连续扩散、跳跃这三件事装进同一个数学框架用一个三元组 (b, σ², ν) 讲清楚一个随机过程到底是怎么走的。做量化定价的、做风险计量的、写保险精算模型的、研究排队系统的甚至做粒子输运和反常扩散研究的最后都会撞上它。这篇文章我不打算从测度论开始堆定义而是按我踩坑的顺序把定义、三元组、从属过程、参数标定、定价落地和模拟误差这几块拆开讲能直接抄的代码我都贴上。1. 为什么连续无跳跃的假设会翻车1.1 布朗运动的两个隐藏代价布朗运动是我们在入门阶段最熟悉的连续时间模型它的路径连续、增量独立、增量服从正态分布。这三条看起来很温和实际上藏着两个致命代价只要把它拿去拟合真实的价格序列问题立刻暴露。第一个代价是路径连续性。连续路径意味着在极小的时间尺度上位移量级只有 O(√Δt)。换句话说无论你把采样频率提到多高单步位移都会随步长缩小而平稳地缩小永远不可能出现隔夜跳 8%这种断崖。可现实里这种事每年都要发生几次。更麻烦的是连续性假设下你没法用跳跃强度这个概念去描述事件到达的密集程度只能被动地靠提高波动率去吸收结果就是平时波动率被高估、关键时刻又被低估。第二个代价来自正态分布的尾部。正态分布的尾部衰减速度是指数平方级的5 倍标准差的事件概率被压到小数点后好多位。但真实数据里日频收益的超额峰度往往在 3 到 10 之间也就是同样的标准差下极端事件的频率比正态预测高出好几个数量级。这两件事叠加起来还会导致第三个副作用方差随时间线性增长年化波动率变成常数任意时间尺度上的分布形状完全一样。可实际数据里期限越短的分布越尖、偏斜越重这跟形状不随时间尺度变化的性质正好相反。1.2 独立增量与平稳增量定义只有四行Lévy 过程的定义其实短得惊人。一个实值随机过程 X {X_t}{t≥0}如果满足下面四条就叫 Lévy 过程X_0 0 几乎必然成立任意 0 ≤ t_1 … t_n各个增量 X{t_i} - X_{t_{i-1}} 相互独立增量平稳也就是 X_{ts} - X_s 的分布只依赖 t 而不依赖 s以及随机连续性对任意 ε 0当 s 趋近 t 时 P(|X_s - X_t| ε) 趋近 0。严格来说第四条在数学上可以由前三条推出来这是 Lévy 本人的一个定理但把它直接写进定义里省事也符合直觉。前两条的含义值得多说一句独立增量意味着过程没有记忆未来怎么走跟过去走成什么样无关平稳增量意味着时间齐次随机机制的规则不随时间改变。这两条合起来就保证了 X_t 的分布一定是单步分布 Q 的 t 次卷积反过来也就推出 Q 必须是无穷可分分布。Lévy 过程和无穷可分分布之间是一一对应的这也是它在极限定理里能扮演高斯分布之于中心极限定理那个角色的根本原因。平稳性这一条在实际建模里需要留个心眼。波动率有日内周期开盘半小时和收盘半小时完全不是一个物种也有区制切换平静期和危机期的机制差别巨大。严格讲日内不同时段的增量分布并不相同直接拿裸数据去拟合会得到一堆自相矛盾的参数。我的一般做法是先做时间去季节化把每个时间点的收益按该时段的历史波动率标准化或者干脆只在同质窗口内做估计。1.3 两个极端之间的一整条光谱布朗运动和泊松过程可以看成两个极端。前者只有无穷多的小位移路径连续方差有限后者只有正跳、幅度恒为定值完全没有连续成分。真正的 Lévy 过程世界是这两个极端之间的一整条光谱而这条光谱上散落着大量在工程上非常好用的模型。过程三元组 (b, σ², ν) 的关键特征路径特点典型场景布朗运动(μ, σ², 0)连续正态增量扩散近似、基准模型泊松过程(0, 0, λδ_1)纯跳跳幅固定事故计数、订单到达复合泊松(0, 0, λF)跳幅服从分布 F单日索赔总额、批量订单Gamma 过程σ² 0ν 集中于正半轴单调不减只有正跳累积降雨量、信用损失累积方差 Gammaσ² 0ν 为两侧指数衰减的 1/|x| 型无穷多小跳 有限大跳日频对数收益α-稳定过程σ² 0ν(dx) ∝ |x|^{-1-α}跳跃测度在零点爆炸方差无穷极端尾部建模Merton 跳扩散(b, σ², λ·N(μ_J, σ_J²))连续扩散 高斯跳权益期权定价CGMYσ² 0ν 为指数阻尼的幂律跳跃活跃度可精细调节高频精细校准这张表里我最想强调的是 α-稳定过程那一行。它是把小跳跃强度无穷大这件事推到极致的结果ν 在零点附近按 |x|^{-1-α} 爆炸导致过程连二阶矩都不存在。这听起来很吓人但恰恰就是很多重尾数据的真实写照。而广义中心极限定理告诉我们独立同分布随机变量和在适当归一化下的极限要么是高斯分布要么是 α-稳定分布——换句话说只要你的数据生成机制允许无穷方差最后收敛到的就是稳定过程。Lévy 过程之所以值得认真学这个定理是最大的理论靠山。2. 三元组把 Lévy-Khintchine 公式读明白2.1 特征函数是最好用的那把尺子Lévy 过程有一个极其漂亮的性质它的特征函数有封闭形式。对任意 t ≥ 0 和 u ∈ RE[e^{iuX_t}] e^{tψ(u)}其中 ψ 叫 Lévy 指数形式是ψ(u) ibu - σ²u²/2 ∫_{R{0}} (e^{iux} - 1 - iux·1_{|x|1}) ν(dx)为什么特征函数能写成 t 次方的指数形式这完全是独立平稳增量的直接后果X_t 是 t 个独立同分布增量的和特征函数自然就是单步特征函数的 t 次方。这个公式最实用的地方在于它告诉你只要知道 (b, σ², ν) 这三个东西过程的全部有限维分布就确定了一整套定价、模拟、估计的工作都可以围绕这三个量展开。那为什么不直接用概率密度而要用特征函数原因很现实绝大多数 Lévy 过程的密度函数没有初等表达式。方差 Gamma 的密度要用第二类贝塞尔函数写CGMY 连贝塞尔函数都写不出来只能用无穷级数或者数值反演。但它们的特征函数几乎总是初等函数形式干净、计算便宜、求导方便。所以整条工程链条——参数估计、期权定价、风险度量、随机数生成——都是围绕特征函数转的。我刚开始做这块的时候总想着先把密度搞出来再说结果在特殊函数里绕了大半个月后来老老实实改走特征函数通道效率提升了不止一个量级。2.2 跳跃测度 ν 和那个不得不加的减号跳跃测度 ν 的读法是ν(A) 表示单位时间内跳跃幅度落进集合 A 的期望次数也就是 E[N((0,t]×A)] t·ν(A)。它不是一个概率测度而是一个强度测度所以 ν(R{0}) 完全可以是无穷大。事实上对于方差 Gamma 和 CGMY 这类模型零点附近的小跳跃强度就是无穷的一秒钟里发生无限多次极其微小的跳跃。但有个约束必须满足∫ min(1, x²) ν(dx) ∞。这条把跳跃测度在无穷远处的增长和在零点附近的增长同时管住了。现在解释那个看起来很别扭的减号 -iux·1_{|x|1}。道理在于幅度大于 1 的大跳跃其总强度 ν({|x|≥1}) 是有限的这部分可以直接处理但幅度小于 1 的小跳跃强度可能是无穷的如果直接写 ∫ x N(dt,dx) 这一项它的期望可能根本不存在或者发散。减掉 iux 相当于对每一层小跳跃先做中心化让积分收敛。生活化的类比是这样的你每天要交无数笔小额手续费单笔金额都很小但笔数无穷多直接把所有笔加起来是发散的必须先对每一笔减去它的期望值这个总偏差才是收敛的、有意义的东西。这个减号的副作用必须记住漂移参数 b 依赖于截断函数的选择。如果你把 1_{|x|1} 换成 1_{|x|c}b 的数值就会变而 σ² 和 ν 保持不变。所以写论文或者写接口文档的时候一定要标注清楚用的是哪种截断约定否则别人照着复现一定对不上。我自己就栽过这个坑照着某篇文献的参数往代码里填模拟出来的路径均值怎么都偏查了整整两天最后发现是文献用的截断和我用的不一致b 差了 ∫_{|x|1} x ν(dx) 那一项。顺带给一个能直接用的一致性关系。在 1_{|x|1} 的约定下拆成漂移 扩散 补偿小跳 大跳之后期望满足 E[X_t] t·(b ∫_{|x|≥1} x ν(dx))因为补偿泊松积分是鞅、均值为零。写完模型先拿这个式子做一次数值校验能挡掉不少低级错误。2.3 Lévy-Itô 分解把过程切成四块如果说特征函数是看 Lévy 过程的方式那 Lévy-Itô 分解就是造 Lévy 过程的方式。它把任意 Lévy 过程唯一地拆成四块X_t bt σW_t ∫_0^t∫_{|x|1} x·Ñ(ds,dx) ∫_0^t∫_{|x|≥1} x·N(ds,dx)四块分别是确定性的线性漂移、标准的连续布朗部分、无穷多次被补偿的小跳跃、有限次不做补偿的大跳跃其中 Ñ N - ν(dx)ds 是补偿泊松随机测度。这四块之间相互独立这个独立性在模拟和蒙特卡洛里非常关键——意味着你可以分开采样再叠加不会引入相关性误差。这个分解最直接的价值是给了一切实操的起点。想模拟分别处理这四项。想推导定价方程把生成元写出来Lf(x) b·f(x) (σ²/2)·f(x) ∫ [f(xy) - f(x) - y·f(x)·1_{|y|1}] ν(dy)有了生成元套用 Feynman-Kac 就能得到跳扩散框架下的定价方程。以权益期权为例价格函数 V(t,S) 满足的是一个偏积分微分方程PIDE比 Black-Scholes 方程多的就是那个积分项∂V/∂t (r - λκ)S·∂V/∂S (σ²S²/2)·∂²V/∂S² ∫ [V(t, S·e^y) - V(t,S)] ν(dy) - rV 0这个积分项是求解时的头号麻烦。有限差分本来能生成稀疏的三对角矩阵加上积分项之后矩阵立刻变稠密直接求逆的代价从 O(N) 涨到 O(N³)。工程上的标准做法是把积分项当成卷积用 FFT 以 O(N log N) 算出来然后配合隐式处理微分算子、显式处理积分算子的 IMEX 时间离散格式。这套组合拳我在做美式期权 PDE 求解的时候用得最多实测稳定性比全隐式好速度比全显式快。3. 从属过程与时间变换最实用的一个技巧3.1 用随机时钟换掉日历时间从属过程subordinator是一类特殊的 Lévy 过程它单调不减只有正跳跳跃测度完全集中在正半轴并且满足比一般情形更强的约束 ∫ min(1, x) ν(dx) ∞。它的拉普拉斯变换有简洁形式 E[e^{-λS_t}] e^{-tΨ(λ)}其中 Ψ(λ) γλ ∫_0^∞ (1 - e^{-λx}) ν(dx) 叫拉普拉斯指数。从属过程最有价值的用法是做时间变换。取一个标准的布朗运动 W再取一个独立的从属过程 S构造 X_t μ·S_t σ·W_{S_t}。这个构造的直觉非常贴切布朗运动本身描述的是每一步都是独立小扰动的扩散但真实的成交和信息到达并不匀速——开盘、数据发布、突发事件的时候市场活跃度暴涨午间流动性枯竭的时候几乎不动。如果用一只有随机指针的时钟去给布朗运动计时活跃时段指针走得快、平静时段走得慢那么最终观察到的日历时间里的价格过程自然就带上了尖峰厚尾。这个物理含义在工程上有个更接地气的版本拿交易笔数或者成交量当业务时钟去替换日历时间很多统计性质会立刻变好。这不是玄学而是因为从属过程本质上就是把事件到达速率随机化而成交量恰好是事件到达速率的一个可观测量。我在处理高频数据的时候用成交量时钟重采样之后的收益率正态性检验的 p 值经常能从 10 的负若干次方提升到 0.05 以上。3.2 方差 Gamma 与 NIG把布朗运动拖慢的两个经典做法最简单的时间变换是拿 Gamma 过程当 S。Gamma 过程 S_t 服从形状参数 t/ν、尺度参数 ν 的 Gamma 分布因此 E[S_t] t、Var[S_t] ν·t。把它代入得到 X_t θ·S_t σ·W_{S_t}这就是方差 Gamma 过程它的特征函数是φ(u) (1 - iθνu σ²νu²/2)^{-t/ν}三个参数的分工非常清楚θ 控制偏度正负都能做ν 控制尾部厚度和超额峰度σ 是整体的扩散尺度。特别要注意 ν 趋近 0 的极限此时 Gamma 时钟退化成确定性的钟方差 Gamma 过程收敛到布朗运动。这个极限正好解释了后面要讲的参数识别噩梦。换成逆高斯从属过程得到的就是 NIG 过程正态逆高斯它和方差 Gamma 是同一个正态方差混合家族的两个成员区别只在于时钟的分布。两者都是纯跳过程没有连续成分但小跳跃强度无穷大所以在有限的分辨率下看起来几乎是连续的。单位时间内的累积量可以拿来做参数含义的定量判断参数表达式说明κ₁θ均值直接给出漂移方向κ₂σ² θ²ν方差扩散与跳动的叠加κ₃3θσ²ν 2θ³ν²三阶累积量偏度的来源κ₄3σ⁴ν 12θ²σ²ν² 3θ⁴ν³四阶累积量峰度的来源超额峰度等于 κ₄ / κ₂²你能从这张表直接看出想让数据显得更尖要么加大 ν要么加大 |θ|。这一步的推演在实际调参时非常省时间比盲试参数高效得多。3.3 方差 Gamma 的模拟代码直接按时间步采样 Gamma 增量是最省事的做法。注意增量独立性让你可以逐段生成然后累加。import numpy as np def simulate_vg(n_steps, T, theta, sigma, nu, n_paths1, seedNone): 模拟方差 Gamma 过程 X_t theta * S_t sigma * W_{S_t} S_t 是 Gamma 从属过程S_t ~ Gamma(shapet/nu, scalenu) 参数: n_steps: 时间步数 T: 总期限 theta: 偏度参数 sigma: 扩散尺度 nu: 尾部厚度参数 (nu - 0 退化为布朗运动) 返回: shape (n_paths, n_steps 1) 的路径数组首列为 0 rng np.random.default_rng(seed) dt T / n_steps if dt / nu 0.05: raise ValueError(dt/nu 过小Gamma 采样会大量下溢请增大步长或调小 nu) shape dt / nu # Gamma 增量的形状参数 scale nu # Gamma 增量的尺度参数 dS rng.gamma(shape, scale, size(n_paths, n_steps)) # 条件在 Gamma 时钟上布朗部分增量的方差是 sigma^2 * dS dW rng.standard_normal((n_paths, n_steps)) * np.sqrt(dS) dX theta * dS sigma * dW X np.cumsum(dX, axis1) return np.concatenate([np.zeros((n_paths, 1)), X], axis1)有一个坑必须提醒当 dt/ν 很小比如 ν 0.1、日频数据 dt 1/250商只有 0.04时Gamma 采样会产生大量接近零的值浮点数下溢之后再开方、再乘正态会引入明显的数值偏差你会看到路径方差不匹配。我一般把商控制在 0.05 以上做不到就改用步长自适应或者用 Gamma 桥在已知终点值的条件下做桥式插值来生成路径。实测下来Gamma 桥在 ν 小于 0.05 的场景里几乎是唯一可靠的选择代价是实现复杂度明显上升。4. 参数标定从数据走到模型4.1 用经验特征函数做拟合既然特征函数是初等形式最顺理成章的估计法就是最小化理论特征函数和经验特征函数之间的距离。这个方法对纯跳过程特别友好因为它们的密度要么难算要么根本不存在。import numpy as np from scipy.optimize import minimize def vg_cf(u, theta, sigma, nu): 方差 Gamma 过程单位时间的特征函数 phi(u) E[exp(i u X_1)] return (1 - 1j * theta * nu * u 0.5 * sigma**2 * nu * u**2) ** (-1.0 / nu) def fit_vg_by_cf(log_returns, u_max50.0, n_u80): 用经验特征函数与理论特征函数的加权距离估计 VG 参数。 log_returns: 单位时间如日频的对数收益序列 r np.asarray(log_returns, dtypefloat) r r - r.mean() # 先去掉均值逐渐比较形状 u_grid np.linspace(0.2, u_max, n_u) def ecf(u): return np.mean(np.exp(1j * u * r)) ecf_vals np.array([ecf(u) for u in u_grid]) # 权重相位信息集中在 |ECF| 较大处用 n/|ECF|^2 做最优加权 n len(r) w n / np.maximum(np.abs(ecf_vals) ** 2, 1e-8) def loss(p): theta, sigma, nu p if sigma 1e-6 or nu 1e-4: return 1e10 diff np.array([vg_cf(u, theta, sigma, nu) - e for u, e in zip(u_grid, ecf_vals)]) return float(np.sum(w * np.abs(diff) ** 2) / n_u) p0 np.array([0.0, r.std(), 0.2]) best None for nu0 in (0.05, 0.2, 0.5, 1.0): # 多初值避开局部极小 res minimize(loss, [p0[0], p0[1], nu0], methodNelder-Mead, options{xatol: 1e-8, fatol: 1e-10, maxiter: 4000}) if best is None or res.fun best.fun: best res theta, sigma, nu best.x # 还原均值到 theta 上单位时间的漂移 return theta r.mean(), sigma, nu, best.fun这段代码有三个地方值得展开讲。第一先减掉均值再拟合形状参数是因为均值和偏度参数 θ 存在直接耦合一起拟合会让目标函数出现长长的浅谷收敛极慢先固定均值、只拟合尺度与尾部最后再把均值加回 θ稳定性好很多。第二权重用 n/|φ̂(u)|² 而不是均匀权重这是有理论依据的渐近最优权重原因是高频部分的经验特征函数被噪声主导其估计偏差量级是 O(1/√n)不降权就会把噪声当信号拟合进去。第三多初值是必须的因为 CGMY 和 VG 这类模型的目标函数经常是多峰的单初值的 Nelder-Mead 在 ν 方向容易卡住不动。u 网格的上限选择也有讲究。日频对数收益的尺度在 0.01 左右主要相位信息集中在 u ≤ 50 这个区间u 再往上理论特征函数的模已经衰减到千分之一以下经验估计的噪声完全盖过信号。我一般的做法是画一条理论 ECF 与经验 ECF 的模随 u 变化的曲线找一个交叉点把 u_max 定在交叉点之前。4.2 累积量匹配快手但脆弱的做法如果你手上只有样本的前几阶矩用累积量匹配能很快得到一组粗略参数特别适合作为上面数值优化的初值。以方差 Gamma 为例记单位时间的样本累积量为 c₁、c₂、c₃、c₄则 θ c₁ν 满足一元二次方程θ³ν² - 3θc₂ν c₃ 0取哪个根要配合约束条件判断必须同时满足 σ² c₂ - θ²ν 0以及判别式 9c₂² ≥ 4θc₃。选好根之后 σ² 直接从 c₂ 里减出来c₄ 留着做一致性检验——如果代入 c₄ 的残差大得离谱说明数据里还有 VG 描述不了的结构比如波动率聚集这时候硬拟合就是自欺欺人。这套方法的短板在于它对高阶矩极其敏感。c₃ 和 c₄ 的样本估计标准误很大几个异常值就能把参数带飞。所以我从来不把它当最终答案只当数值优化的起跳点。真正靠谱的还是 4.1 节那种用完整特征函数曲线做最小化的方式。4.3 从数据里把跳跃抠出来模型参数估计完之后还有一个反向验证的问题数据里到底有没有跳跃这里的经典工具是双幂变差。它的出发点很聪明——已实现方差 RV Σr_i² 在存在跳跃时会包含跳跃变差而双幂变差 BV (π/2)·Σ|r_i||r_{i-1}| 对跳跃稳健只捕捉连续部分的积分方差。两者之差就是跳跃变差。import numpy as np def jump_variation(r, k1): 用双幂变差把已实现方差拆成连续变差和跳跃变差。 r: 高频对数收益已去均值 k: 滞后阶数通常取 1 r np.asarray(r, dtypefloat) absr np.abs(r) RV float(np.sum(r ** 2)) BV float((np.pi / 2.0) * np.sum(absr[k:] * absr[:-k])) JV max(RV - BV, 0.0) ratio JV / RV if RV 0 else 0.0 return {RV: RV, BV: BV, JV: JV, jump_share: ratio}我在用这个指标时有两条血泪经验。第一微观结构噪声会同时抬高 RV 和 BV在采样频率高于几秒的时候噪声项会彻底压过跳跃信号所以必须先做子采样或者噪声校正否则你会检测出一大堆根本不存在的跳跃。第二jump_share 在不同交易日之间波动极大单日的结果没有意义必须看滚动窗口的均值和分布再配合跨资产比较才能判断跳跃是不是一个稳定的模型特征。如果某个资产的 jump_share 长年在 5% 以下那给它配一个纯跳模型多半是过度设计老老实实上跳扩散甚至纯扩散更划算。5. 落到业务里定价、风控与排队5.1 特征函数定价一条公式吃遍所有模型Lévy 过程在衍生品定价里最爽的一点在于不管你选哪个具体模型定价流程是完全一样的因为大家都有特征函数。欧式看涨期权价格可以写成C S₀·Π₁ - K·e^{-rT}·Π₂其中 Π₂ Q(S_T K) 是风险中性测度下的行权概率Π₁ Q^S(S_T K) 是份额测度下的行权概率两者都能用 Gil-Pelaez 反演公式从特征函数直接算出来Π_j 1/2 (1/π)∫₀^∞ Re[e^{-iu·log K}·φ_j(u)/(iu)] du。份额测度的特征函数通过测度变换得到φ^S(u) φ(u - i)/φ(-i)。import numpy as np def call_price_by_cf(S0, K, r, T, cf, u_max200.0, n_u2**15): 用 Gil-Pelaez 反演计算欧式看涨期权价格。 cf: 特征函数 phi(u) E[exp(i u X_T)]X_T log(S_T / S0) 必须已经做过鞅修正满足 cf(-1j) exp(r*T) k np.log(K / S0) u np.linspace(1e-8, u_max, n_u) phi_u cf(u) # 风险中性测度 phi_ui cf(u - 1j) # 份额测度 phi_mi cf(-1j) # 应等于 exp(r*T) integ2 np.real(np.exp(-1j * u * k) * phi_u / (1j * u)) integ1 np.real(np.exp(-1j * u * k) * phi_ui / (1j * u * phi_mi)) Pi1 0.5 np.trapezoid(integ1, u) / np.pi Pi2 0.5 np.trapezoid(integ2, u) / np.pi return S0 * Pi1 - K * np.exp(-r * T) * Pi2鞅修正是这一段最关键也最容易漏的一步。以方差 Gamma 为例模型原始形式给出的对数收益特征函数对应的期望不是无风险利率你必须加一个漂移补偿量 ω满足ω r (1/ν)·ln(1 - θν σ²ν/2)这里的约束 1 - θν σ²ν/2 0 必须成立否则 E[S_T] 发散整个定价框架失效。我第一次做这块的时候忘了加 ω算出来的看涨期权价格比看跌还便宜对着屏幕愣了半天才反应过来。所以接口里最好加一句断言直接检查 cf(-1j) 是不是等于 exp(r*T)事后再检查价格合理性就晚了。u_max 和 n_u 的取值决定精度。深虚值期权对 u_max 极其敏感u_max 太小会把高频信息截掉价格错得离谱。判断方法很简单把 u_max 翻倍再算一遍看价格变化量是否落在你能容忍的范围内比如 0.01%。如果振荡明显就说明被积函数衰减太慢这时候需要上阻尼技巧——把看涨价格乘上 e^{αk}α 是正阻尼参数让被积函数在复平面上衰减得更快再配合 FFT 一次性算出一整条行权价曲线。做市系统的报价表基本都是这么生成的单次 FFT 算出几百个行权价比逐个反演快很多。5.2 保险精算里的破产概率Lévy 过程在保险里的用法甚至比金融领域更直接。经典的 Cramér-Lundberg 模型把保险公司的盈余写成 U_t u ct - Σ_{i1}^{N_t} Y_i其中 u 是初始资本c 是保费率N 是强度 λ 的泊松过程Y_i 是独立同分布的索赔额。这个 U 本身就是个 Lévy 过程——线性漂移加上一个复合泊松跳跃没有扩散项。在这个框架下安全负荷率 ρ 决定了保费率 c (1 ρ)λE[Y]调节系数 R 由方程 λ(M_Y(R) - 1) cR 唯一确定其中 M_Y 是索赔额的矩母函数。Lundberg 不等式给出破产概率的上界 ψ(u) ≤ e^{-Ru}指数衰减的速度直接由 R 控制。推广方向有两个都很自然。一是把复合泊松换成一般的从属过程得到更灵活的 Lévy 保险风险模型能描述索赔频率本身也随机变化的情形二是往盈余过程里加一个布朗扰动项得到扩散扰动风险模型用来描述保费收入的连续波动和投资收益。这两条推广在文献里统称Lévy 风险过程好处是所有分析工具都从特征函数和拉普拉斯指数直接继承过来不需要重新推一遍。有意思的是这里的 Pollaczek-Khinchine 公式和排队论里的稳态等待时间公式本质上是同一个式子。盈余过程的破产概率、M/G/1 队列的稳态工作量、反射 Lévy 过程在零点的局部时三者在数学上是一套东西的三个投影。我个人觉得理清这层关系比背十个公式有用得多。5.3 反射过程与排队反射 Lévy 过程的标准构造是 Q_t X_t L_t其中 L 是 Q 在零点处的局部时负责把过程推回非负半轴。它的离散版本就是 Lindley 递归 W_{n1} max(W_n X_n, 0)几乎每个写过排队模拟的人都用过。连续版本的价值在于稳态工作量的拉普拉斯变换可以用 Lévy 过程的谱理论写成显式形式从而避开长时间的瞬态模拟。这一块我印象最深的是数值稳定性问题。反射过程在接近零点时局部时会迅速累积如果模拟时步长太粗你会看到路径在零附近反复穿越局部时被严重高估稳态分布的右尾被压扁。我的做法是做步长收敛测试把步长减半两次比较稳态分布的 95% 分位数变化超过百分之一就继续细化。这个测试花不了几分钟但能省掉后面大量的返工。6. 模拟误差控制真正花钱的地方6.1 小跳跃截断的偏差从哪来对于跳跃测度在零点附近无界的模型α-稳定、CGMY 的 Y 0、方差 Gamma你不可能真的把每一次微小跳跃都模拟出来。标准做法是设一个截断阈值 ε只模拟 |x| ε 的跳跃把 |x| ≤ ε 的部分近似成扩散加漂移。具体地补偿漂移用小跳跃的一阶矩 b_ε ∫_{|x|≤ε} x ν(dx)扩散系数用二阶矩 σ_ε² ∫_{|x|≤ε} x² ν(dx)。这个替换带来两层误差。一层是矩的不匹配被替换掉的真实小跳跃有非零的三阶及以上累积量而替代的高斯扩散三阶以上累积量为零所以分布的形状偏了。另一层是随机性结构的差异一阶和二阶矩匹配了但跳跃在时间上的聚集特性和扩散完全不同。α-稳定过程下截断偏差随 ε 的收敛速度比较慢通常是 O(ε^{2-α}) 量级这意味着 α 越接近 2收敛越慢需要越小的 ε。工程上的做法很简单把 ε 依次取 1e-3、5e-4、2.5e-4算三次价格看这个序列的模式。如果差值按预期比例缩小说明已经在渐近区间内如果两次之间跳动很大而且方向不定说明还没进渐近区需要更小的 ε 或者检查代码是不是有 bug。我一般把这个测试写成一个 pytest 用例每次改动模拟器都跑一遍防止某次重构悄悄把截断逻辑改坏。6.2 方差缩减三个我常用的招第一个招是条件蒙特卡洛。以方差 Gamma 为例条件在 Gamma 时钟的整条路径 S 上资产价格的对数收益是正态的欧式期权的条件期望可以直接用 Black-Scholes 公式解析算出蒙特卡洛只需要处理从属过程这一低维随机源。这一步能把方差降低一到两个数量级是所有技巧里性价比最高的。对路径依赖的期权只能做到部分条件化但即使只条件化终点附近的时段方差下降也很明显。第二个招是重要性抽样。给跳跃测度做指数倾斜Esscher 变换把稀有的大跳跃拉到分布中心附近再抽样同时给每条路径乘上对应的似然比权重。这个技巧在计算深虚值期权和极端分位数风险度量比如 99.9% 违约阈值的时候几乎是必需的因为裸蒙特卡洛在这种场景下的有效样本量可能小到个位数。有意思的是Esscher 变换既是定价时的测度变换工具又是模拟时的方差缩减工具同一套数学两个用途。第三个招是控制变量。用连续部分把跳跃全部关掉的贴现收益作为控制变量因为它的解析期望已知而且和全模型的收益高度相关。实测下来对短期期权能把方差降低 50% 到 70%实现成本还很低属于白捡的收益。7. 常见问题与排查速查7.1 参数识别为什么会拟合出假的σ这是我在 Lévy 过程建模里遇到频率最高的问题。症状很典型拟合出来的 ν 一路往零跑σ 稳定在样本标准差附近θ 稳定在样本均值附近最终结果退化成布朗运动。原因不是优化器坏了而是数据本身没法区分无穷多极小跳跃和连续扩散。有限样本下两者对特征函数在这段区间的贡献几乎不可分目标函数在 ν 方向是一条长长的平底谷。诊断方法很简单固定 θ 和 σ把 ν 从 0 到 1 扫一遍画轮廓目标函数曲线。如果曲线在 ν 小于某个值之后基本是平的就说明这个参数在这个数据上不可识别。处理手段有三条按我实际用过的优先级排列换成更高频的数据小跳跃在细时间尺度上更突出识别力提升明显把 ν 固定到文献里的典型量级只估计剩下的参数或者上贝叶斯方法用弱信息先验把参数约束在一个合理范围让后验分布而不是点估计来回答问题。最后这条在样本量小的时候特别好用因为后验分布会把这个参数其实不太确定这件事诚实地表现出来而不是给你一个看起来很精确的假数字。7.2 数值反演的三个典型故障第一个故障是深虚值期权价格明显失真原因一般是 u_max 设得太小。判断方法是把 u_max 翻倍重算如果价格变化超过阈值就说明之前截断了有效信息。千万不要凭感觉设 u_max一定要做收敛测试。第二个故障是价格曲线出现上下振荡。它来源于被积函数随 u 衰减太慢离散化时的截断误差在反演积分里形成了类似 Gibbs 现象的振荡。解法是引入 Carr-Madan 的阻尼因子把被积函数乘上一个在虚轴方向增长、在实轴方向快速衰减的函数让尾部行为变好同时把积分网格加密配合 Simpson 权重而不是简单的梯形权重。第三个故障是模型隐含的远期价格对不上市场远期。这几乎总是鞅修正漏了或者做错了。记住那个检查点cf(-1j) 必须等于 e^{rT}。如果分红率、外汇双币种利率这些都在里面就得把对应项补全否则整个定价链路会在最基础的一致性上崩掉。7.3 常见问题速查表现象可能原因我的处理方式ν 一路趋近 0σ ≈ 样本标准差参数不识别小跳跃与扩散不可分换高频数据、固定 ν、或改用贝叶斯方法深虚值价格失真反演积分上限 u_max 太小翻倍 u_max 做收敛测试价格曲线上下振荡被积函数衰减慢、网格太粗Carr-Madan 阻尼因子 加密网格模拟路径方差与理论不符时间步长过粗、Gamma 采样下溢增大步长或改用 Gamma 桥E[S_T] 与 e^{rT} 对不上漏了漂移补偿 ω按 ω r ln(1-θνσ²ν/2)/ν 修正跳跃检测结果异常多微观结构噪声被当成跳跃子采样、噪声校正、提高阈值CGMY 拟合不收敛C 与 Y 强相关固定 Y 做网格搜索再优化其余参数截断阈值 ε 缩小时价格不收敛还没进渐近区或代码有 bug依次减半 ε 观察收敛模式这张表里的每一条我都真实踩过尤其是C 与 Y 强相关那条让一个模型调了两天才发现不是优化器的问题而是参数本身在数据长度不够的时候根本分不开。经验是CGMY 这种四参数模型日频数据至少要有三年以上、并且包含至少一次明显的跳跃事件否则不要指望能把四个参数都稳定估出来。最后分享一个我自己的习惯。每次上线一个新模型之前我一定做三件事拿布朗运动做退化工况测试把跳跃参数推向极限看代码能否平滑退化到 Black-Scholes 解析解拿一条人工构造的已知参数路径做回环测试生成数据、再拟合、看能不能把参数找回来以及拿市场上真实报价做一次反算隐含参数的一致性检查。这三步花掉的时间和事后排查一个生产事故相比便宜得不像话。踩过几次坑之后我现在看到任何拟合结果特别漂亮的模型第一反应都是先怀疑参数识别和数据质量而不是先庆祝。