4、协整关系建模：Johansen协整检验、协整向量估计、误差修正模型（ECM）、半衰期计算

协整关系，说白了就是两个或多个非平稳时间序列，它们的线性组合却是平稳的。这在统计套利里是个核心概念——你想想看，两只股票各自随机游走，但它们的价差却像被一根橡皮筋牵着，偏离太远总会弹回来。这就是我们赚钱的基础。

我个人习惯把协整建模拆成四个步骤：先检验是否存在协整关系，再估计出具体的协整向量，然后建立误差修正模型来刻画短期动态，最后算一下半衰期，看看这个均值回复的速度有多快。下面我们一个一个来过。

4.1 Johansen协整检验

Engle-Granger两步法虽然简单，但有个硬伤——它只能处理两个变量，而且你得事先指定哪个是因变量。Johansen检验就灵活多了，它能同时检验多个变量之间是否存在多个协整关系。

Johansen检验的核心思想，是基于向量自回归（VAR）模型的。假设我们有一个k维的时间序列向量 y_t，我们可以把它写成：

Δy_t = Π y_{t-1} + Σ Γ_i Δy_{t-i} + ε_t

这里的Π矩阵是关键。它的秩r就决定了协整关系的个数。如果r=0，说明没有协整关系；如果r=k，说明所有变量都是平稳的；如果0 < r < k，那就存在r个协整关系。

Johansen给出了两种检验统计量：

• 迹检验（Trace Test）：H0: 协整关系个数 ≤ r，H1: 协整关系个数 > r
• 最大特征值检验（Max Eigenvalue Test）：H0: 协整关系个数 = r，H1: 协整关系个数 = r+1

我在项目中遇到过一个问题：两个检验结果不一致怎么办？嗯，我的经验是优先看迹检验的结果，因为它对备择假设更稳健。当然，你也可以结合经济意义来判断。

下面是一个Python实现的例子：

import statsmodels.tsa.vector_ar.vecm as vecm
import pandas as pd

# 假设我们有三只股票的价格序列
data = pd.DataFrame({
    'stock_a': prices_a,
    'stock_b': prices_b,
    'stock_c': prices_c
})

# 确定最优滞后阶数
# 这里用AIC准则，你也可以试试BIC
model_select = vecm.select_order(data, maxlags=10)
print(model_select.summary())

# 进行Johansen协整检验
# 注意：det_order=0表示没有确定性趋势，1表示有常数项
johansen_result = vecm.coint_johansen(data, det_order=1, k_ar_diff=2)

# 输出迹统计量和临界值
print("迹统计量:", johansen_result.lr1)
print("95%临界值:", johansen_result.cvt)
print("最大特征值统计量:", johansen_result.lr2)
print("95%临界值:", johansen_result.cvm)

4.2 协整向量估计

一旦Johansen检验确认存在协整关系，下一步就是估计协整向量。说白了，就是找到那个线性组合的系数，让组合后的序列变得平稳。

Johansen方法通过求解一个广义特征值问题来估计协整向量。具体来说，它会给出所有可能的协整向量，并按特征值从大到小排序。特征值越大，说明该协整关系越显著。

在Python中，我们可以直接从Johansen检验的结果中提取协整向量：

# 提取协整向量（特征向量）
# 注意：statsmodels返回的是转置后的矩阵
eigenvectors = johansen_result.evec
print("协整向量:\n", eigenvectors)

# 通常我们取第一个协整向量（对应最大特征值）
coint_vector = eigenvectors[:, 0]
print("第一个协整向量:", coint_vector)

# 计算价差序列
spread = data.dot(coint_vector)
print("价差序列的前10个值:\n", spread.head(10))

你可能会问：协整向量需要归一化吗？我的建议是，归一化到第一个变量上，这样解释起来更直观。比如，如果第一个变量是股票A，归一化后协整向量就变成了(1, -β₂, -β₃)，这样价差就是 A - β₂B - β₃C。

# 归一化到第一个变量
normalized_vector = coint_vector / coint_vector[0]
print("归一化后的协整向量:", normalized_vector)

4.3 误差修正模型（ECM）

协整关系告诉我们长期均衡是什么，但短期偏离后怎么调整回来？这就是误差修正模型（ECM）要解决的问题。

ECM的基本形式是这样的：

Δy_t = α * (y_{t-1} - β'x_{t-1}) + γ * Δx_t + ε_t

这里的 (y_{t-1} - β'x_{t-1}) 就是上一期的误差修正项，也就是价差。α是调整速度系数，它告诉我们当价差偏离均衡时，系统会以多快的速度往回拉。

我个人习惯把ECM理解为「一个带纠错机制的动态模型」。你想想看，如果价差昨天偏离了均衡，今天就会有一部分修正回来。α的绝对值越大，修正速度越快。

下面是用Python估计ECM的代码：

from statsmodels.tsa.vector_ar.vecm import VECM

# 构建VECM模型
# k_ar_diff是差分项的滞后阶数，这里设为1
vecm_model = VECM(data, k_ar_diff=1, coint_rank=1, deterministic='ci')
vecm_result = vecm_model.fit()

# 查看结果摘要
print(vecm_result.summary())

# 提取误差修正系数（alpha）
alpha = vecm_result.alpha
print("误差修正系数:\n", alpha)

# 提取短期动态系数（gamma）
gamma = vecm_result.gamma
print("短期动态系数:\n", gamma)

4.4 半衰期计算

半衰期，就是价差偏离均衡后，恢复到一半所需的时间。这个指标在统计套利里特别实用——它直接告诉我们，这个策略的持仓周期大概有多长。

计算半衰期的方法很简单。假设价差序列服从一个均值回复过程：

Δspread_t = θ * (μ - spread_{t-1}) + ε_t

其中θ是回复速度。半衰期就是：

half_life = ln(2) / θ

注意，这里的θ是负的（因为均值回复），所以半衰期是正数。

具体计算时，我们可以对价差做一阶自回归：

import numpy as np
from sklearn.linear_model import LinearRegression

# 计算价差
spread = data.dot(normalized_vector)

# 构建回归：Δspread_t = α + β * spread_{t-1} + ε_t
spread_lag = spread.shift(1).dropna()
spread_diff = spread.diff().dropna()

# 对齐数据
y = spread_diff.values
X = spread_lag.values.reshape(-1, 1)

# 线性回归
model = LinearRegression()
model.fit(X, y)
beta = model.coef_[0]

# 计算半衰期
half_life = np.log(2) / (-beta)
print(f"半衰期: {half_life:.2f} 个交易日")

我曾经遇到过一个案例，半衰期算出来是200多天。一开始我还以为代码写错了，后来仔细一看，原来是两只股票的基本面已经发生了变化——一家公司被收购了。所以，半衰期异常长的时候，一定要回头检查一下协整关系是否还成立。

知识体系总览

下面这张图总结了协整关系建模的完整流程：

整个流程环环相扣。先通过Johansen检验确认协整关系存在，然后估计出具体的协整向量，接着用ECM刻画短期调整机制，最后算半衰期来指导实际交易。每一步都有坑，但每一步也都有对应的解决办法。

嗯，协整关系建模就讲到这里。记住，统计套利的核心不是找到完美的协整关系，而是找到那些协整关系稳定、半衰期适中的配对。剩下的，就是严格执行交易纪律了。