Python实战Merton跳跃扩散模型：量化分析商品期货隐藏的24%尾部风险

如果你还在依赖传统的几何布朗运动（GBM）模型来评估商品期货的风险敞口，那么你可能正严重低估市场的真实破坏力。

原油、天然气这类大宗商品的价格波动从来都不“温和”。它们时刻受到供应冲击、地缘政治突发事件和流动性危机的洗礼，由此产生的剧烈价格“跳跃”是标准正态分布根本无法捕捉的“肥尾”风险。

本文将带你用Python完整实现一个Merton跳跃扩散（MJD）模型的校准与模拟引擎，并通过对比揭示传统模型的致命缺陷。整个流程涵盖：

1. 数据清洗：精准识别并处理期货合约换月带来的“假跳跃”噪声。
2. 参数校准：使用最大似然估计（MLE）拟合跳跃扩散过程的五个核心参数。
3. 蒙特卡洛模拟：利用向量化操作，高效生成上千条未来价格路径。
4. 风险对比：量化比较GBM与MJD模型在衡量极端尾部风险时的巨大差距。

最终的结论值得所有风险管理者警惕：标准GBM模型可能低估了约24%的尾部风险。

一、为何传统模型在商品市场失效？

1.1 GBM模型的根本局限

几何布朗运动（GBM）假设资产价格的变化是连续且平滑的，其收益率服从正态分布。然而，真实的商品市场充满了不连续性：

• 突发的供应中断（如OPEC意外减产）
• 地缘政治事件导致的瞬间暴涨暴跌
• 收益率分布呈现出明显的“尖峰厚尾”特征

这些由离散事件驱动的价格跳跃，是GBM模型完全无法描述的。

1.2 实践中的两个常见误区

误区一：忽视肥尾，低估风险。使用GBM计算的风险价值（VaR）或预期损失（CVaR）会严重偏乐观，导致风险准备金不足。
误区二：数据噪声误判为跳跃。直接使用拼接的期货连续合约数据时，会将合约换月产生的确定性价差误认为是随机跳跃，从而污染模型校准。

二、数据清洗：屏蔽合约换月噪声

当你从yfinance获取CL=F（WTI原油期货）这类数据时，得到的是由多个不同到期月份合约拼接而成的连续序列。合约换月时产生的价差是确定性的展期收益或成本，而非我们模型要捕捉的随机市场跳跃。

2.1 5-Sigma窗口检测法

我们的清洗策略是：在每个月度的最后3个交易日内，将收益率超过5倍标准差的数据点识别为换月噪声，并将其收益率置零（而非删除，以保持时间序列的连续性）。

import numpy as np
import pandas as pd
import yfinance as yf

class Config:
    """全局配置参数"""
    PRIMARY_TICKER = 'CL=F'  # 原油期货
    START_DATE = "2020-01-01"
    END_DATE = "2025-12-18"
    ROLL_WINDOW_DAYS = 3      # 换月窗口天数
    OUTLIER_SIGMA = 5.0       # 异常值阈值
    DT = 1 / 252.0            # 日时间步长

def fetch_and_clean_data(ticker: str) -> pd.DataFrame:
    """
    获取期货数据并清洗换月噪声

    参数:
        ticker: 期货代码
    返回:
        清洗后的 DataFrame
    """
    # 下载原始数据
    df = yf.download(
        ticker,
        start=Config.START_DATE,
        end=Config.END_DATE,
        auto_adjust=False,
        progress=False
    )

    # 计算对数收益率
    df['Log_Returns'] = np.log(df['Close'] / df['Close'].shift(1))
    std_returns = df['Log_Returns'].std()

    # 标记换月窗口：每月最后 3 个交易日
    df['Is_Roll_Window'] = False
    month_ends = df.index.to_series().groupby(pd.Grouper(freq='ME')).max()

    for me in month_ends:
        # 获取该月末前 3 天的索引
        window_indices = df.index[df.index <= me][-Config.ROLL_WINDOW_DAYS:]
        df.loc[window_indices, 'Is_Roll_Window'] = True

    # 在换月窗口内，屏蔽超过 5 倍标准差的收益率
    mask_condition = (
        (df['Is_Roll_Window']) & 
        (np.abs(df['Log_Returns']) > Config.OUTLIER_SIGMA * std_returns)
    )

    outlier_count = mask_condition.sum()
    print(f"检测到 {outlier_count} 个换月噪声点，已屏蔽")

    # 将噪声收益率置零（而非删除，保持时间连续性）
    df.loc[mask_condition, 'Log_Returns'] = 0.0

    return df.dropna(subset=['Log_Returns'])

# 执行数据清洗
clean_df = fetch_and_clean_data(Config.PRIMARY_TICKER)
print(f"清洗后数据点数: {len(clean_df)}")

运行结果：

检测到 2 个换月噪声点，已屏蔽
清洗后数据点数: 1500

三、Merton跳跃扩散模型核心原理

3.1 模型思想

Merton跳跃扩散模型（MJD）巧妙地将价格运动分解为两个独立部分：

• 连续扩散部分：描述资产价格的常规随机波动，与GBM模型中的布朗运动相同。
• 离散跳跃部分：由一个泊松过程驱动，用来刻画突然发生的、幅度随机的价格冲击。

模型涉及五个关键参数：

参数	含义
$\mu$	资产收益的漂移率
$\sigma$	扩散部分的波动率
$\lambda$	跳跃强度（单位时间内跳跃发生的平均次数）
$\mu_J$	跳跃幅度的均值
$\sigma_J$	跳跃幅度的标准差

3.2 不可或缺的补偿项

为了确保模型的期望收益率与漂移率参数$\mu$一致，必须引入一个补偿项$k$：
$k = e^{\mu_J + \frac{1}{2}\sigma_J^2} - 1$
这个项会调整漂移部分，以抵消跳跃带来的平均影响。

四、最大似然估计（MLE）校准

MJD模型的对数收益率概率密度函数是一个无限混合分布：由无数个具有不同均值和方差的正态分布，按其泊松概率加权求和。在实际计算中，我们需要对跳跃次数进行截断（例如，求和到k=10次），并使用logsumexp技巧来保证数值稳定性。

import math
from scipy.stats import norm
from scipy.special import logsumexp
from statsmodels.base.model import GenericLikelihoodModel

class MertonMLE(GenericLikelihoodModel):
    """
    Merton 跳跃扩散模型的最大似然估计器

    参数向量: [mu, sigma, lambda, mu_j, sigma_j]
    """

    def __init__(self, endog, dt=1/252, max_jumps=10):
        super().__init__(endog)
        self.dt = dt              # 时间步长
        self.max_jumps = max_jumps  # 泊松求和截断项数

    def nloglikeobs(self, params):
        """计算每个观测点的负对数似然"""
        mu, sigma, lamb, mu_j, sigma_j = params
        y = self.endog  # 对数收益率序列
        dt = self.dt

        # 参数有效性检查
        if sigma <= 1e-6 or lamb < 0 or sigma_j <= 1e-6:
            return np.ones_like(y) * 1e10

        log_probs = []

        # 对跳跃次数 k 从 0 到 max_jumps 求和
        for k in range(self.max_jumps):
            # 泊松权重（对数形式）
            log_weight = (
                -lamb * dt + 
                k * np.log(np.maximum(lamb * dt, 1e-12)) - 
                np.log(float(math.factorial(k)))
            )

            # 给定 k 次跳跃时的正态分布参数
            mean_k = (mu - 0.5 * sigma**2) * dt + k * mu_j
            var_k = sigma**2 * dt + k * sigma_j**2

            # 正态分布的对数概率密度
            lp = log_weight + norm.logpdf(y, loc=mean_k, scale=np.sqrt(var_k))
            log_probs.append(lp)

        # 使用 logsumexp 合并（保证数值稳定）
        total_log_ll = logsumexp(np.vstack(log_probs), axis=0)

        return -total_log_ll  # 返回负对数似然

def calibrate_mjd(returns: np.ndarray) -> np.ndarray:
    """
    使用 MLE 校准 Merton 跳跃扩散模型

    参数:
        returns: 对数收益率数组
    返回:
        校准后的参数数组 [mu, sigma, lambda, mu_j, sigma_j]
    """
    dt = Config.DT

    # 初始参数猜测
    start_params = [
        np.mean(returns) / dt,      # mu: 年化漂移
        np.std(returns) / np.sqrt(dt),  # sigma: 年化波动率
        5.0,                        # lambda: 初始跳跃强度
        0.0,                        # mu_j: 跳跃均值
        np.std(returns)             # sigma_j: 跳跃波动率
    ]

    # 创建模型并拟合
    model = MertonMLE(returns, dt=dt)
    result = model.fit(
        start_params=start_params,
        method='nm',       # Nelder-Mead 算法，对非凸优化更稳健
        maxiter=1000,
        disp=0
    )

    return result.params

# 校准模型
returns = clean_df['Log_Returns'].values
mle_params = calibrate_mjd(returns)

# 打印校准结果
param_names = ['漂移率 (μ)', '扩散波动率 (σ)', '跳跃强度 (λ)',
               '跳跃均值 (μ_j)', '跳跃波动率 (σ_j)']
print("\n========== 校准结果 ==========")
for name, value in zip(param_names, mle_params):
    print(f"{name}: {value:.4f}")

校准结果示例：

========== 校准结果 ==========
漂移率 (μ): 0.3275
扩散波动率 (σ): 0.3212
跳跃强度 (λ): 17.5882
跳跃均值 (μ_j): -0.0148
跳跃波动率 (σ_j): 0.0859

结果解读：模型从原油期货历史数据中识别出，平均每年会发生约17.6次价格跳跃。跳跃均值$\mu_J$为负，表明下行冲击（暴跌）的频率或幅度略高于上行冲击（暴涨）。

五、向量化蒙特卡洛模拟

为了高效地生成大量（例如1000条）未来价格路径以进行风险分析，我们采用NumPy的向量化操作，彻底避免低效的Python循环。

def simulate_mjd_paths(
    S0: float,
    params: np.ndarray,
    steps: int = 252,
    paths: int = 1000
) -> np.ndarray:
    """
    向量化生成 MJD 价格路径

    参数:
        S0: 初始价格
        params: 模型参数 [mu, sigma, lambda, mu_j, sigma_j]
        steps: 模拟步数（天数）
        paths: 路径数量
    返回:
        价格矩阵，形状为 (steps+1, paths)
    """
    mu, sigma, lamb, mu_j, sigma_j = params
    dt = Config.DT

    # 计算补偿项 k
    k = np.exp(mu_j + 0.5 * sigma_j**2) - 1

    # 补偿后的漂移
    drift = (mu - 0.5 * sigma**2 - lamb * k) * dt

    # 生成扩散部分（标准正态）
    Z = np.random.standard_normal((steps, paths))
    diffusion = sigma * np.sqrt(dt) * Z

    # 生成跳跃部分
    # 泊松分布决定每个时间步的跳跃次数
    N = np.random.poisson(lamb * dt, (steps, paths))
    # 跳跃幅度 = 次数 × 均值 + 根号(次数) × 标准差 × 正态随机数
    J = N * mu_j + np.sqrt(N) * sigma_j * np.random.standard_normal((steps, paths))

    # 合并对数增量
    log_increments = drift + diffusion + J

    # 累积求和并转换为价格
    path_matrix = S0 * np.exp(np.cumsum(log_increments, axis=0))

    # 在首行添加初始价格
    path_matrix = np.vstack([np.ones(paths) * S0, path_matrix])

    return path_matrix

# 获取最新价格
current_price = clean_df['Close'].iloc[-1]

# 模拟 MJD 路径
mjd_paths = simulate_mjd_paths(current_price, mle_params, steps=252, paths=1000)
print(f"生成 {mjd_paths.shape[1]} 条路径，每条 {mjd_paths.shape[0]} 个时间点")

六、风险度量与模型对比

6.1 计算VaR与CVaR

我们采用条件风险价值（CVaR），也称预期损失，来衡量尾部风险。CVaR衡量的是在损失超过风险价值（VaR）的条件下的平均损失，比VaR更能反映极端损失的严重程度。

def calculate_risk_metrics(paths: np.ndarray, confidence: float = 0.99):
    """
    计算 VaR 和 CVaR（预期损失）

    参数:
        paths: 价格路径矩阵
        confidence: 置信水平
    返回:
        (VaR, CVaR) 元组
    """
    # 计算终端收益率
    final_returns = (paths[-1] / paths[0]) - 1

    # VaR: 给定置信水平下的最大损失分位点
    alpha = 1 - confidence
    var = np.percentile(final_returns, alpha * 100)

    # CVaR: 超过 VaR 损失的平均值
    cvar = final_returns[final_returns <= var].mean()

    return var, cvar

# 计算 MJD 模型的风险
mjd_var, mjd_cvar = calculate_risk_metrics(mjd_paths)

# 生成 GBM 基准路径（跳跃参数设为 0）
gbm_params = np.array([mle_params[0], mle_params[1], 0.0, 0.0, 0.0])
gbm_paths = simulate_mjd_paths(current_price, gbm_params, steps=252, paths=1000)
gbm_var, gbm_cvar = calculate_risk_metrics(gbm_paths)

# 结果对比
print("\n========== 风险度量对比 (99% 置信水平) ==========")
print(f"{'模型':<20} {'VaR':<15} {'CVaR (预期损失)':<15}")
print("-" * 50)
print(f"{'MJD (跳跃扩散)':<20} {mjd_var:.2%}       {mjd_cvar:.2%}")
print(f"{'GBM (标准模型)':<20} {gbm_var:.2%}       {gbm_cvar:.2%}")
print("-" * 50)
print(f"风险低估差距: {abs(mjd_cvar - gbm_cvar):.2%}")

输出结果：

========== 风险度量对比 (99% 置信水平) ==========
模型                  VaR            CVaR (预期损失)
--------------------------------------------------
MJD (跳跃扩散)        -56.23%        -66.75%
GBM (标准模型)        -35.12%        -42.41%
--------------------------------------------------
风险低估差距: 24.34%

6.2 结果解读

指标	GBM 模型	MJD 模型	差距
CVaR (99%)	-42.41%	-66.75%	24.34%

这个差距是触目惊心的。它意味着，如果你基于传统的GBM模型来设置风险限额、计算保证金或评估投资组合风险，那么当真实世界中发生一次剧烈的价格跳跃事件时，你所遭受的实际损失可能会比你预期的最坏情况还要高出24个百分点。对于杠杆交易者而言，这种程度的低估足以导致爆仓。

七、完整代码示例

以下是将所有步骤整合在一起的完整脚本，方便你直接运行和实验。

"""
Merton 跳跃扩散模型完整实现
用于商品期货的风险分析
"""

import numpy as np
import pandas as pd
import yfinance as yf
import math
from scipy.stats import norm
from scipy.special import logsumexp
from statsmodels.base.model import GenericLikelihoodModel

# ==================== 配置 ====================
class Config:
    PRIMARY_TICKER = 'CL=F'
    START_DATE = "2020-01-01"
    END_DATE = "2025-12-18"
    ROLL_WINDOW_DAYS = 3
    OUTLIER_SIGMA = 5.0
    DT = 1 / 252.0
    SIM_PATHS = 1000
    SIM_STEPS = 252
    CONFIDENCE = 0.99

# ==================== 数据处理 ====================
def fetch_and_clean_data(ticker: str) -> pd.DataFrame:
    """获取并清洗期货数据"""
    df = yf.download(ticker, start=Config.START_DATE,
                     end=Config.END_DATE, auto_adjust=False, progress=False)
    df['Log_Returns'] = np.log(df['Close'] / df['Close'].shift(1))
    std_returns = df['Log_Returns'].std()

    # 标记换月窗口
    df['Is_Roll_Window'] = False
    month_ends = df.index.to_series().groupby(pd.Grouper(freq='ME')).max()
    for me in month_ends:
        window_indices = df.index[df.index <= me][-Config.ROLL_WINDOW_DAYS:]
        df.loc[window_indices, 'Is_Roll_Window'] = True

    # 屏蔽换月噪声
    mask = (df['Is_Roll_Window']) & (np.abs(df['Log_Returns']) > Config.OUTLIER_SIGMA * std_returns)
    df.loc[mask, 'Log_Returns'] = 0.0

    return df.dropna(subset=['Log_Returns'])

# ==================== MLE 校准 ====================
class MertonMLE(GenericLikelihoodModel):
    def __init__(self, endog, dt=1/252, max_jumps=10):
        super().__init__(endog)
        self.dt = dt
        self.max_jumps = max_jumps

    def nloglikeobs(self, params):
        mu, sigma, lamb, mu_j, sigma_j = params
        y = self.endog
        dt = self.dt

        if sigma <= 1e-6 or lamb < 0 or sigma_j <= 1e-6:
            return np.ones_like(y) * 1e10

        log_probs = []
        for k in range(self.max_jumps):
            log_weight = -lamb * dt + k * np.log(np.maximum(lamb * dt, 1e-12)) - np.log(float(math.factorial(k)))
            mean_k = (mu - 0.5 * sigma**2) * dt + k * mu_j
            var_k = sigma**2 * dt + k * sigma_j**2
            lp = log_weight + norm.logpdf(y, loc=mean_k, scale=np.sqrt(var_k))
            log_probs.append(lp)

        return -logsumexp(np.vstack(log_probs), axis=0)

def calibrate_mjd(returns):
    start_params = [np.mean(returns)/Config.DT, np.std(returns)/np.sqrt(Config.DT), 5.0, 0.0, np.std(returns)]
    model = MertonMLE(returns, dt=Config.DT)
    result = model.fit(start_params=start_params, method='nm', maxiter=1000, disp=0)
    return result.params

# ==================== 模拟与风险计算 ====================
def simulate_paths(S0, params, steps, paths):
    mu, sigma, lamb, mu_j, sigma_j = params
    dt = Config.DT
    k = np.exp(mu_j + 0.5 * sigma_j**2) - 1
    drift = (mu - 0.5 * sigma**2 - lamb * k) * dt

    Z = np.random.standard_normal((steps, paths))
    diffusion = sigma * np.sqrt(dt) * Z
    N = np.random.poisson(lamb * dt, (steps, paths))
    J = N * mu_j + np.sqrt(N) * sigma_j * np.random.standard_normal((steps, paths))

    log_increments = drift + diffusion + J
    path_matrix = S0 * np.exp(np.cumsum(log_increments, axis=0))
    return np.vstack([np.ones(paths) * S0, path_matrix])

def calculate_cvar(paths):
    final_returns = (paths[-1] / paths[0]) - 1
    var = np.percentile(final_returns, (1 - Config.CONFIDENCE) * 100)
    return final_returns[final_returns <= var].mean()

# ==================== 主程序 ====================
if __name__ == "__main__":
    # 1. 数据准备
    df = fetch_and_clean_data(Config.PRIMARY_TICKER)
    returns = df['Log_Returns'].values
    S0 = df['Close'].iloc[-1]

    # 2. 模型校准
    params = calibrate_mjd(returns)
    print(f"跳跃强度: {params[2]:.2f} 次/年")

    # 3. 风险对比
    mjd_cvar = calculate_cvar(simulate_paths(S0, params, Config.SIM_STEPS, Config.SIM_PATHS))
    gbm_cvar = calculate_cvar(simulate_paths(S0, np.array([params[0], params[1], 0, 0, 0]), Config.SIM_STEPS, Config.SIM_PATHS))

    print(f"\nMJD CVaR: {mjd_cvar:.2%}")
    print(f"GBM CVaR: {gbm_cvar:.2%}")
    print(f"风险低估: {abs(mjd_cvar - gbm_cvar):.2%}")

总结

本文构建了一个从数据清洗到风险量化的完整Merton跳跃扩散模型分析流程，其核心价值点在于：

数据清洗：提出了基于时间窗口和标准差阈值的过滤器，有效剔除了期货合约换月带来的结构性噪声，为模型校准提供了“干净”的数据基础。

模型校准：实现了基于最大似然估计（MLE）的稳健参数拟合方法，其中跳跃强度λ是识别市场是否存在显著跳跃风险的关键指标。

风险度量：采用向量化蒙特卡洛模拟进行高效路径生成，并强调使用CVaR而非VaR来更准确地衡量极端尾部风险的潜在损失。

核心发现：在原油期货的实证分析中，传统的GBM模型严重低估了风险，其预期损失（CVaR）比MJD模型的结果低了约24%。这一差距在风险管理和资金配置中是不可忽视的。

对于希望将此模型应用于生产环境的从业者，建议进一步采取滚动窗口校准策略以动态捕捉市场状态的变迁，并可以考虑使用期货交易所公布的精确展期价差来进行更高级的数据调整。

本文探讨的数据科学方法在量化金融与风险管理中具有广泛的应用前景。如果你想与其他开发者交流此类技术的实现细节或挑战，云栈社区的相关板块是一个不错的去处。