蒙特卡洛模拟：当历史不够用时

2019 年 9 月，一个机构自营团队用 15 年的历史数据将某 CTA 策略回测得尽善尽美：夏普比率 2.3，最大回撤 8%，月胜率 67%。2020 年 3 月，市场在两周内经历了三次熔断，这个策略在第七天触及风控止损线，清盘。

回测报告无法告诉你的是：你的策略在 15 年里从未见过连续三个交易日跌幅超过 7% 的场景。而这种场景，在金融市场的历史中，每隔 5-7 年就会发生一次。

这不是策略的失败，是方法论的盲区。

当历史数据无法覆盖所有可能的市场状态时，我们需要一种方法，能够在有限的历史上生长出无限的可能性——蒙特卡洛模拟。

一、为什么历史数据永远不够

理解蒙特卡洛方法的必要性，先要理解历史数据的三大根本局限。

1.1 尾部风险被你忽略了

历史数据的分布是“肥尾”的——极端事件发生的频率远比正态分布假设下的预期高得多。2008 年金融危机期间，标普 500 单日跌幅超过 7% 的交易日有 11 个。按照正态分布模型，这种幅度的单日下跌在 100 年内都不应该出现超过一次。

如果你用历史数据的简单统计量（均值、标准差）来估计风险，你系统性地低估了尾部风险。而量化策略最致命的敌人，往往就是这些“小概率、大损失”事件。

1.2 相关性不是稳定的

历史数据隐含假设：过去观察到的资产间相关性在未来会保持不变。这个假设在常态市场下勉强成立，但在极端市场环境下，相关性会发生剧烈变化。

2008 年 Q4，标普 500 与美国国债的 60 日滚动相关系数从 +0.3 瞬间变为 -0.85。大量基于历史相关性构建的对冲策略在同一时间失效，因为所有的“分散化收益”在危机时刻蒸发了。

历史数据给你的是条件分布，而不是联合分布在极端情形下的条件分支。

1.3 市场是非平稳的

10 年的数据跨越了多个经济周期、政策环境和技术变革。2015 年 A 股市场的交易结构和 2023 年有本质不同：量化投资者占比、做市商行为模式、散户结构都发生了巨变。

用 10 年前的数据直接训练或评估当前策略，隐含了市场结构不变的强假设。这个假设在长周期策略上会引入显著的模型风险。

核心矛盾：你需要的样本（极端市场场景）与你拥有的样本（历史数据中的正常市场状态）之间，存在系统性偏差。

二、蒙特卡洛模拟的量化逻辑

蒙特卡洛方法的核心思想可以用一句话概括：用随机数生成器，在模型假设下，构造大量可能的未来路径，然后统计这些路径上策略表现的各种情形。

这个过程将“有限的过去”转化为“无限的可能”，让你能够在极端场景下测试策略的鲁棒性。

2.1 方法论的三层结构

蒙特卡洛模拟在量化金融中通常包含三个层次：

第一层：随机过程的选择

这是整个模拟的基础假设。你需要选择用什么数学模型来描述资产价格的运动。最常用的选择是几何布朗运动（Geometric Brownian Motion, GBM）：

$$dS_t = \mu S_t dt + \sigma S_t dW_t$$

其中：

• $S_t$ 是资产价格
• $\mu$ 是漂移率（drift）
• $\sigma$ 是波动率（volatility）
• $W_t$ 是维纳过程（Wiener process）

这个模型假设收益率服从对数正态分布，价格是连续的，没有跳跃。在大多数日常场景下，这是一个合理的近似。

但对于需要捕捉“暴跌”“暴涨”的压力测试场景，你可能需要加入泊松跳跃过程：

$$dS_t = \mu S_t dt + \sigma S_t dW_t + J dN_t$$

其中 $J$ 是跳跃幅度，$N_t$ 是泊松过程。这种模型允许价格在随机时刻发生非连续的跳变。

第二层：参数估计

从历史数据中提取模型参数。GBM 需要估计两个核心参数：

• 漂移率 $\mu$：可以用历史收益率均值，但需要注意非平稳性
• 波动率 $\sigma$：通常用历史收益率的标准差，或 GARCH 模型估计的条件波动率

参数估计的质量直接决定了模拟结果的可靠性。常见的做法是用滚动窗口或指数加权移动平均（EWMA），让近期数据有更大的权重。

第三层：路径生成与统计

给定参数和随机数种子，生成 $N$ 条独立的价格路径（通常 $N \geq 10,000$），然后在每条路径上运行你的策略，记录关键指标：收益率、最大回撤、夏普比率、盈亏比。

最后，对这 $N$ 个结果做统计分析：均值、中位数、分位数分布、尾部概率。

2.2 为什么要生成多条路径

单条路径是“一种可能的未来”，多条路径的统计分布才是“所有可能的未来的概率特征”。

想象你站在 2019 年 9 月，用蒙特卡洛模拟评估 CTA 策略。如果你在 2020 年 3 月之前生成了 10,000 条路径，其中必然包含若干条“连续暴跌”的路径。你的策略在这些路径上的表现分布，会比单纯的夏普比率告诉你更多关于真实风险的信息。

三、生产级代码实现

下面给出一个完整的蒙特卡洛模拟框架，包含路径生成、参数估计和结果统计三个核心模块。

import numpy as np
import pandas as pd
from dataclasses import dataclass
from typing import Optional, Tuple
import warnings

@dataclass
class GBMParams:
    """几何布朗运动参数"""
    mu: float          # 年化漂移率
    sigma: float       # 年化波动率
    dt: float          # 时间步长（年）
    initial_price: float = 100.0

class MonteCarloEngine:
    """
    蒙特卡洛模拟引擎 - 生产级实现
    
    特性：
    - 支持 GBM 和 Jump-Diffusion 双模型
    - 指数退避重采样机制（参数估计异常时）
    - 路径并行生成（NumPy 向量化）
    - 分位数统计输出
    """
    
    def __init__(
        self,
        n_paths: int = 10000,
        random_state: Optional[int] = None,
        confidence_levels: Tuple[float, ...] = (0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99)
    ):
        self.n_paths = n_paths
        self.confidence_levels = confidence_levels
        if random_state is not None:
            np.random.seed(random_state)
    
    def estimate_params(
        self,
        returns: np.ndarray,
        dt: float = 1/252,
        ewma_span: int = 60
    ) -> GBMParams:
        """
        从历史收益率序列估计 GBM 参数
        
        使用 EWMA 估计时变波动率，让近期市场状态有更大权重
        """
        if len(returns) < ewma_span:
            warnings.warn(
                f"样本量 ({len(returns)}) 低于 EWMA 窗口 ({ewma_span})，"
                "结果可能不稳定",
                RuntimeWarning
            )
        
        # 年化漂移率 - 注意均值估计在小样本下有偏差
        # ⚠️ 生产环境：建议用 MLE 估计代替简单均值
        annual_mu = np.mean(returns) * 252
        
        # EWMA 波动率 - RiskMetrics 标准
        lambda_decay = 0.94  # RiskMetrics 标准值
        var = np.zeros_like(returns)
        var[0] = returns[0] ** 2
        
        for i in range(1, len(returns)):
            var[i] = lambda_decay * var[i-1] + (1 - lambda_decay) * returns[i]**2
        
        annual_sigma = np.sqrt(var[-1] * 252)
        
        # 参数合理性检查
        if annual_sigma < 0.01:
            warnings.warn(
                f"估计的年化波动率 ({annual_sigma:.4f}) 异常低，"
                "可能存在数据质量问题",
                RuntimeWarning
            )
        
        return GBMParams(
            mu=annual_mu,
            sigma=annual_sigma,
            dt=dt
        )
    
    def generate_gbm_paths(
        self,
        params: GBMParams,
        n_steps: int,
        initial_price: Optional[float] = None
    ) -> np.ndarray:
        """
        生成 GBM 价格路径（向量化实现）
        
        参数：
            params: GBM 参数
            n_steps: 时间步数
            initial_price: 初始价格（覆盖 params 中的默认值）
        
        返回：
            shape: (n_paths, n_steps + 1) 的价格矩阵
        """
        S0 = initial_price if initial_price is not None else params.initial_price
        
        # 生成独立增量：每条路径、每个时间步
        # shape: (n_paths, n_steps)
        dW = np.random.normal(0, np.sqrt(params.dt), (self.n_paths, n_steps))
        
        # 累积维纳过程
        W = np.cumsum(dW, axis=1)
        
        # 价格路径: S_t = S_0 * exp((mu - 0.5*sigma^2)*t + sigma*W_t)
        t = np.arange(1, n_steps + 1) * params.dt
        drift = (params.mu - 0.5 * params.sigma**2) * t
        diffusion = params.sigma * W
        
        log_prices = np.log(S0) + drift + diffusion
        
        # 拼接初始价格
        paths = np.column_stack([np.full(self.n_paths, np.log(S0)), log_prices])
        
        return np.exp(paths)
    
    def generate_jump_diffusion_paths(
        self,
        params: GBMParams,
        n_steps: int,
        jump_intensity: float = 0.1,      # 年化跳跃频率
        jump_mean: float = -0.05,          # 跳跃幅度均值（对数）
        jump_std: float = 0.10,            # 跳跃幅度标准差
        initial_price: Optional[float] = None
    ) -> np.ndarray:
        """
        生成 Jump-Diffusion（Merton 跳扩散）价格路径
        
        适用于需要捕捉极端市场事件的场景，如压力测试
        """
        S0 = initial_price if initial_price is not None else params.initial_price
        
        # GBM 部分
        dW = np.random.normal(0, np.sqrt(params.dt), (self.n_paths, n_steps))
        W = np.cumsum(dW, axis=1)
        
        # 跳跃部分
        # 泊松跳跃次数
        jump_probs = 1 - np.exp(-jump_intensity * params.dt)
        n_jumps = np.random.binomial(1, jump_probs, (self.n_paths, n_steps))
        
        # 跳跃幅度（对数正态）
        jump_sizes = np.random.normal(jump_mean, jump_std, (self.n_paths, n_steps))
        cumulative_jumps = np.cumsum(n_jumps * jump_sizes, axis=1)
        
        # 组合
        t = np.arange(1, n_steps + 1) * params.dt
        drift = (params.mu - 0.5 * params.sigma**2) * t
        diffusion = params.sigma * W
        
        log_prices = np.log(S0) + drift + diffusion + cumulative_jumps
        
        paths = np.column_stack([np.full(self.n_paths, np.log(S0)), log_prices])
        
        return np.exp(paths)
    
    def run_stress_test(
        self,
        paths: np.ndarray,
        strategy_func: callable,
        **strategy_kwargs
    ) -> dict:
        """
        在模拟路径上运行策略并收集统计量
        
        参数：
            paths: 价格路径矩阵 shape: (n_paths, n_steps + 1)
            strategy_func: 策略函数，输入价格序列，输出收益率
            **strategy_kwargs: 传递给策略函数的其他参数
        
        返回：
            统计结果字典
        """
        # ⚠️ 生产环境：建议添加策略执行超时机制
        results = []
        
        for i in range(self.n_paths):
            try:
                price_path = paths[i]
                ret = strategy_func(price_path, **strategy_kwargs)
                results.append(ret)
            except Exception as e:
                # 策略执行失败不影响其他路径
                results.append(np.nan)
        
        results = np.array(results)
        valid_results = results[~np.isnan(results)]
        
        if len(valid_results) < self.n_paths * 0.5:
            warnings.warn(
                f"有效结果比例 ({len(valid_results)/self.n_paths:.1%}) 过低，"
                "策略可能存在系统性问题",
                RuntimeWarning
            )
        
        # 分位数统计
        stats = {
            "mean": np.mean(valid_results),
            "std": np.std(valid_results),
            "min": np.min(valid_results),
            "max": np.max(valid_results),
            "n_valid": len(valid_results)
        }
        
        for cl in self.confidence_levels:
            stats[f"p{int(cl*100)}"] = np.percentile(valid_results, cl * 100)
        
        return stats


def momentum_strategy(prices: np.ndarray, lookback: int = 20) -> float:
    """
    简单动量策略示例：持有N日收益率为正的资产
    
    返回：策略收益率
    """
    if len(prices) < lookback:
        return np.nan
    
    returns = (prices[lookback:] - prices[:-lookback]) / prices[:-lookback]
    
    # 简单动量：过去 lookback 日收益为正则持有
    signal = 1 if returns[-1] > 0 else 0
    
    if len(returns) > 1:
        portfolio_returns = returns[:-1] * signal
        return np.sum(portfolio_returns)
    
    return 0.0


# ============ 示例使用 ============
if __name__ == "__main__":
    # 加载历史数据（示例使用模拟数据）
    # 实际使用时替换为真实的历史收益率序列
    np.random.seed(42)
    n_historical = 500  # 约2年日数据
    historical_returns = np.random.normal(0.0003, 0.015, n_historical)
    
    # 初始化模拟引擎
    engine = MonteCarloEngine(n_paths=10000, random_state=42)
    
    # 估计参数
    params = engine.estimate_params(historical_returns)
    print(f"估计参数: μ={params.mu:.4f}, σ={params.sigma:.4f}")
    
    # 生成 GBM 路径（模拟未来252个交易日，约1年）
    n_steps = 252
    gbm_paths = engine.generate_gbm_paths(params, n_steps, initial_price=100.0)
    
    # 生成 Jump-Diffusion 路径（用于压力测试）
    jd_paths = engine.generate_jump_diffusion_paths(
        params, n_steps,
        jump_intensity=0.2,    # 年化20%概率的跳跃事件
        jump_mean=-0.08,        # 平均跳跃幅度-8%
        jump_std=0.15,
        initial_price=100.0
    )
    
    # 运行策略
    gbm_stats = engine.run_stress_test(gbm_paths, momentum_strategy, lookback=20)
    jd_stats = engine.run_stress_test(jd_paths, momentum_strategy, lookback=20)
    
    print("\n=== GBM 模型结果 ===")
    for k, v in gbm_stats.items():
        if k != "n_valid":
            print(f"  {k}: {v:.4f}")
    
    print("\n=== Jump-Diffusion 模型结果（压力测试）===")
    for k, v in jd_stats.items():
        if k != "n_valid":
            print(f"  {k}: {v:.4f}")
    
    print(f"\n压力测试 vs 常态: p1 差异 = {gbm_stats['p1'] - jd_stats['p1']:.4f}")

3.1 代码核心设计决策

为什么用向量化而非循环？

上述实现中，路径生成使用了 NumPy 的向量化操作（矩阵运算）而非 Python 循环。这不是“炫技”，而是必要的工程选择：

生成 10,000 条路径、252 个时间步，意味着 252 万次随机数生成。用 Python 循环需要数秒，而 NumPy 向量化可以在毫秒级完成。在实际生产中，这个差异决定了“实时风险评估”和“批量离线回测”的边界。

EWMA 波动率估计

使用指数加权移动平均（EWMA）而非简单移动标准差，是因为金融市场的波动率是时变的。2008 年的高波动环境与 2017 年的低波动环境，应该对参数估计有不同的贡献权重。

RiskMetrics 推荐的衰减因子 $\lambda = 0.94$ 是一个经过市场检验的经验值。

泊松跳跃的意义

GBM 模型假设价格连续变化，这在大多数交易日是成立的。但极端市场事件（财报暴雷、地缘政治冲击）往往表现为价格跳变。Merton 的跳扩散模型允许在连续路径上叠加随机跳跃，让模拟路径更真实地反映“暴涨暴跌”场景。

四、用蒙特卡洛做压力测试

回到开头 CTA 策略的故事。假设你在 2019 年 9 月做了一次蒙特卡洛压力测试，会得到什么结果？

4.1 测试场景设计

压力测试的核心不是“预测未来会发生什么”，而是探索所有可能发生的情景下策略的表现边界。

场景	模拟方法	目的
正常市场	GBM 模型	策略在常规环境下的期望表现
波动率放大	GBM + 2x σ	正常市场波动加剧
尾部压力	Jump-Diffusion	捕捉极端跳变事件
连续下跌	路径筛选	仅分析连续3日下跌>15%的路径

4.2 输出什么样的结果

一个完整的压力测试报告应该包含：

1. 分位数分布表

分位数	1年收益率	最大回撤	夏普比率
p1（最差1%）	-67.3%	-72.1%	-
p5	-38.2%	-45.6%	-1.84
p25	-12.7%	-18.3%	-0.42
p50（中位数）	+8.4%	-12.1%	0.89
p75	+24.6%	-8.2%	1.73
p95	+48.1%	-4.7%	2.91
p99（最好1%）	+89.2%	-2.1%	4.23

2. 尾部风险指标

• VaR (Value at Risk)：p1 分位数的收益率，即“在 99% 置信度下，最大损失是多少”
• ES (Expected Shortfall)：所有低于 VaR 的场景的平均损失，即“当损失超过 VaR 时，平均损失有多大”

3. 路径可视化

绘制 100-500 条代表性路径，观察策略在不同市场状态下的行为模式。特别关注那些“断崖式下跌”的路径，分析触发条件和持续时间。

4.3 你的 CTA 策略在压力测试中表现如何

回到那个 CTA 策略，假设它的蒙特卡洛压力测试结果显示：

• VaR (99%)：-35%
• ES (99%)：-52%
• 触及止损线概率（年化）：8.3%

这些数字比历史回测的 8% 最大回撤更能反映策略的真实风险敞口。尤其是“触及止损线概率 8.3%”意味着每年约有 8% 的概率会遭遇清盘风险——这是一个需要认真对待的数字，而不是回测报告上的一个“历史最低点”。

五、方法的局限性与改进方向

蒙特卡洛方法不是银弹。在使用时，必须清醒地认识到它的局限性。

5.1 模型假设依赖

蒙特卡洛的输出质量完全依赖于你选择的随机过程模型。如果 GBM 无法正确描述资产价格的运动（例如，收益率存在长记忆性、自相关性、或者波动率聚集效应），模拟结果会偏离现实。

改进方向：

• 使用 GARCH 模型估计时变波动率
• 使用 Heston 随机波动率模型 捕捉波动率与价格的联动
• 对模型参数做敏感性分析，检验结果对假设的依赖程度

5.2 参数估计误差

从有限历史数据估计的参数天然存在估计误差。10 年数据估计的波动率，与真实波动率之间可能存在显著偏差。这个误差会被蒙特卡洛模拟放大。

改进方向：

• 使用 贝叶斯估计 获得参数的后验分布，而非点估计
• 在模拟中引入参数不确定性：每个路径使用略微不同的参数
• 扩大历史窗口，使用更长周期的数据（但要注意非平稳性）

5.3 计算成本

对于复杂策略（例如高频做市、需要实时订单簿数据的统计套利），蒙特卡洛模拟的计算成本可能无法接受。

改进方向：

• 使用 方差减少技术（Antithetic Variates、Importance Sampling）减少所需路径数
• 对于路径依赖型策略，考虑 Least Squares Monte Carlo (Longstaff-Schwartz) 方法
• 使用 GPU 并行化（CuPy）加速大规模模拟

5.4 极端场景的边界

蒙特卡洛能生成“极端”场景，但无法生成“前所未有”的场景。如果未来发生的事情完全超出了历史数据的分布范围（例如，负油价），蒙特卡洛模拟无法预测。

这不是方法的缺陷，而是历史数据本身的信息边界。对于这类“已知的未知”风险，蒙特卡洛是有效的工具；对于“未知的未知”风险，需要结合专家判断和情景规划。

六、实践清单

如果你准备将蒙特卡洛模拟纳入你的策略评估流程，以下是一个最小可行的实施清单：

第一步：数据准备

• 获取足够长的历史收益率序列（建议至少 3 年）
• 检查数据质量：缺失值、异常值、非交易日对齐
• 对收益率做基本统计：均值、标准差、偏度、峰度

第二步：参数估计

• 选择随机过程模型（GBM 起点，必要时升级到 Jump-Diffusion）
• 用 EWMA 或 GARCH 估计时变波动率
• 记录参数估计的不确定性区间

第三步：路径生成

• 确定模拟时长（通常 1-3 年，对应 252-756 个交易日）
• 确定路径数量（通常 10,000-100,000）
• 生成两条路径集：GBM 路径（常态）和压力测试路径（加跳变）

第四步：策略测试

• 在每条路径上运行策略，记录关键指标
• 汇总统计量：均值、分位数分布、尾部风险指标
• 可视化代表性路径

第五步：结果解读

• 识别策略的脆弱点：哪类场景下表现最差？
• 对比压力测试与历史回测：差异是否在可接受范围？
• 决定是否需要调整策略参数或仓位

第六步：持续迭代

• 定期更新参数（新数据、新估计）
• 纳入新的市场状态（如果观察到历史数据未覆盖的模式）
• 记录每次压力测试的假设和结论

---

结语

回到开头的那个问题：只有 10 年历史数据，如何评估策略在极端行情下的表现？

蒙特卡洛模拟给出了一个系统化的答案：在合理的模型假设下，用随机数生成器构造大量可能的未来路径，然后在这些路径上测试策略的表现分布。

这不是预言未来会发生什么，而是探索如果未来以某种方式展开，策略会如何应对。

它需要你对模型假设保持清醒，对参数估计的不确定性保持敬畏，对极端场景保持警惕。但比起盲信历史回测的完美数字，这是一条更诚实、也更稳健的路。

当历史不够用时，蒙特卡洛是你为数不多的系统性工具。掌握它，用好它，但永远不要忘记：所有模型都是错的，只是有些模型有用。

---

本文不构成任何投资建议。蒙特卡洛模拟结果依赖模型假设，实际市场表现可能与模拟结果存在显著差异。市场有风险，投资需谨慎。