合成价格路径与相关收益

周五晚上,某私募量化研究员要对一个 20 只股票的行业轮动策略做半年回测,需要一个 (T=252, N=20) 的日收益矩阵。问题是平台的合规决策写得很清楚:不接行情数据牌照,所有训练样例只能跑合成数据。CSV 里没有,卖方接口也没有,只能自己生成。这一课给出最小可复现的配方:一颗确定的随机种子、对数欧拉离散化的 GBM 一步、用 Cholesky 分解构造相关矩阵、再留一根扣子接学生 t 分布与跳跃。产出是一个 pd.DataFrame,带 DatetimeIndex,可以直接喂给 3.2.2 L5 的业绩归因 pipeline 和 3.2.3 L4 的最小方差组合 pipeline。

单资产 GBM 的对数欧拉一步

连续时间的几何布朗运动 SDE 是 $dS_t = \mu S_t \, dt + \sigma S_t \, dW_t$,其中 $W_t$ 是标准布朗运动。本课只用它的离散化配方;存在唯一性定理、欧拉-丸山方案的强 1/2 阶 / 弱 1 阶收敛性证明全在 2.4.x 章节的随机微积分里。对一个交易日,取时间步 $dt = 1/252$ 年,抽取 $T$ 个标准正态扰动 $Z_t$,然后应用对数欧拉更新:

公式里那个 -0.5 * sigma**2 * dt 是伊藤修正项,推导见 2.4.1。少了它,离散化后的漂移每步会向上偏 0.5 * sigma**2 * dt,半年下来叠出来的偏差肉眼可见。算术欧拉形式 $S_{t+1} = S_t + \mu S_t \, dt + \sigma S_t \sqrt{dt}\, Z_t$ 在粗格子上可能给出负价格;对数形式不会。新代码一律默认对数形式。

import numpy as np
rng = np.random.default_rng(seed=42)
dt = 1.0 / 252
n_days = 252
mu, sigma, S0 = 0.08, 0.20, 100.0
Z = rng.standard_normal(n_days)
log_returns = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z
prices = S0 * np.exp(np.cumsum(log_returns))

把 rng = np.random.default_rng(seed=42) 当成本模块的全局纪律:每一次随机抽样都从这只显式 Generator 走,SciPy 的 .rvs(...) 用 random_state=rng 关键字把同一只 Generator 传进去。永远不要再写 np.random.seed / np.random.randn 这种全局状态 API,种子一旦从全局漏出去,失败的测试会变成幽灵 bug。下面这句话是不容妥协的可复现性纪律:两次调用 simulate_basket(..., seed=42) 必须产出逐字节相同的 pd.DataFrame,用 pd.testing.assert_frame_equal(a, b) 验证。

用 Cholesky 构造相关多资产收益

目标:一个形状 (T, N) 的标准正态扰动矩阵,其 N x N 的样本相关矩阵,在 T 增大时收敛到事先指定的相关矩阵 R。标准配方:抽 (T, N) 的独立同分布标准正态,左乘协方差矩阵的下三角 Cholesky 因子,再让漂移与波动率按列广播。R 必须对称且严格正定:R == R.T 且所有特征值大于零,验证手段是对应的 3.2.3 L4 对称特征分解一课。构造方法本质上是高斯协方差矩阵的方根分解,得到的协方差矩阵在构造上等于 L @ L.T = R。

import scipy.linalg as sla
N = 5
R = np.full((N, N), 0.30); np.fill_diagonal(R, 1.0)
L = sla.cholesky(R, lower=True)
Z = rng.standard_normal((n_days, N))
Z_corr = Z @ L.T

接下来把 Z_corr 翻译成每只资产的对数收益:log_returns = drift_vector * dt + vol_vector * sqrt(dt) * Z_corr,这里 drift_vector 与 vol_vector 都是形状 (N,) 的一维向量,沿轴 1 广播过去,严格依赖 3.2.1 L1 教过的广播规则。再做一次取指数再累乘:prices = S0 * np.exp(np.cumsum(log_returns, axis=0)),就拿到了价格矩阵。

一行代码切换到学生 t 分布的厚尾扰动

把高斯抽样换成学生 t 抽样,只动一行;但要补一句方差校正,否则厚尾的同时方差也膨胀了。学生 t 在自由度 $\nu$ 处的方差是 $\nu/(\nu - 2)$,只在 $\nu > 2$ 时有定义,所以做公平对比要乘 sqrt((df - 2) / df) 把方差扳回 1。

import scipy.stats as st
df = 4
Z = st.t.rvs(df=df, size=(n_days, N), random_state=rng) * np.sqrt((df - 2) / df)

3.2.3 L1 教过的高斯分布偏度峰度诊断在这里反过来用:用 df=4 模拟出来的收益,反向拟合应该回到 $\nu_t \approx 4$;超额峰度按 $6/(\nu - 4)$ 计算,在 $\nu = 4$ 处发散,因此面板回归需要的拟合峰度,建议起步用 df=5 或 df=6,先把数值稳住再调厚尾强度。

可选的跳跃扰动

跳跃扩散路径在对数收益里加一项复合泊松扰动:每一步抽 $N_{jumps} \sim \text{Poisson}(\lambda \cdot dt)$,典型 $\lambda \approx 1$ 次/年,所以 lam * dt ≈ 1/252,每步几乎不跳;真跳的时候,跳幅从 $N(\mu_J, \sigma_J)$ 抽出来再加到这一步的对数收益。这是简化的 Merton 跳跃扩散模型,完整模型(包含风险中性测度下的漂移修正、欧式期权闭式定价公式、隐含波动率曲面校准)归 1.4 衍生品 / 2.4 随机微积分。这里它只是一个开关,用来制造像样的隔夜跳空,服务于风险演示。

打包成 simulate_basket

def simulate_basket(n_days: int, tickers: list[str], mu: np.ndarray, sigma: np.ndarray, R: np.ndarray, seed: int, S0: float = 100.0, start: str = '2024-01-02') -> pd.DataFrame:
    import pandas as pd
    import scipy.linalg as sla
    rng = np.random.default_rng(seed)
    dates = pd.bdate_range(start=start, periods=n_days)
    L = sla.cholesky(R, lower=True)
    Z = rng.standard_normal((n_days, len(tickers)))
    Z_corr = Z @ L.T
    dt = 1.0 / 252
    log_returns = mu * dt - 0.5 * sigma**2 * dt + sigma * np.sqrt(dt) * Z_corr
    prices = S0 * np.exp(np.cumsum(log_returns, axis=0))
    return pd.DataFrame(prices, index=dates, columns=tickers)

参数校验:len(tickers) == len(mu) == len(sigma) == R.shape[0] == R.shape[1],并且 R 在 np.allclose 容差内对称。本课的工作篮子用 {600519.SH, 000001.SZ, 600036.SH, 510300.SH, 510500.SH},以沪深300 ETF (510300.SH) 作大盘指数代理、中证500 ETF (510500.SH) 作小盘指数代理,带宽中证 500 是与沪深 300 高度相关、与单票相关一般的双层结构。说明性参数显式写死:mu = np.array([0.08, 0.06, 0.06, 0.05, 0.07])(年化漂移),sigma = np.array([0.30, 0.25, 0.25, 0.18, 0.22])(年化波动率),R 对角线 1.0,非对角线默认 0.30,但 510300.SH ↔ 510500.SH 这对索引代理写 0.60。这些数字只为教学清晰,不对应任何具体市场的校准。

公式 -0.5 * sigma**2 * dt 在对数欧拉 GBM 步里是伊藤修正项;省了它会让离散化漂移每步向上偏 0.5 * sigma**2 * dt。底层连续过程是 布朗运动,其离散化产生 GBM;独立同分布的扰动来自高斯分布,其样本均值收敛到正态在 中心极限定理 的脚注里说明;Cholesky 分解的目标是单位方差的协方差矩阵(即相关矩阵 R)。完成后,做一次健全性检查:

S_{t+1} = S_t \cdot \exp((\mu - 0.5\sigma^2) dt + \sigma\sqrt{dt}\, Z_t)

国内 A 股有春节周与十一长假的非交易日缺口,合成数据这里为简化直接用 pd.bdate_range,严格的交易日历讨论留给 3.2.2 L4。simulated_returns.corr() 与目标 R 的偏差用 np.allclose(..., atol=0.1) 比较,容差给得宽,是因为 252 天样本的样本相关矩阵本身就有抽样误差;若要更细的容差,3.2.3 L2 的 bootstrap CI 工具是首选。

参考阅读:SciPy 官方用户指南 linalg + stats 章节;NumPy 官方用户指南 random.Generator 章节;《Python 金融大数据分析》(Yves Hilpisch 中译,第二版)第 12 章关于随机过程与蒙特卡洛模拟的实现章节;国内常用 pip 镜像(清华 / 阿里)安装 scipy / pyarrow 的实务提示。

下一课把 N 从固定的 5 只扩到 200 只的合成横截面:行业、市值、因子分都按真实分布抽,再把每只股票的 CAPM + 因子模型期望收益接回本课的 simulate_basket,顺手再做一份带买卖价差与市场冲击的合成 tick 流,服务于成本演示。