时变参数 TVP-VAR 模型：把「系数的关系也会变」做成可估计的量 • halo 的技术博客

因子的 beta 上周还说 0.6，这周突然塌到 0.2；美股和原油的相关，平静时是弱正、危机时是强正。这些「变量之间的关系本身在随时间变」的现象，是量化里最容易被常系数模型忽略、却最致命的一类——你以为模型拟合得好，其实它把一段缓慢漂移、中间还可能跳一下的系数，强行压成了一个平均数。

最朴素的应对，是切一段窗口做滚动 OLS：每 40 天重跑一次回归，看系数怎么走。它的问题很硬：窗口内的系数是常数假设，跳跃被平滑成一段斜坡；而且它对近期变化天然滞后约半个窗长（≈20 期）——等你「看到」系数变了，市场已经走完半程。

结论先放这：时变参数 VAR（TVP-VAR）把每个系数当成一条随时间演化的状态，用卡尔曼滤波逐期更新，既能紧跟漂移、又能定位跳跃，还天然给出不确定度。 在我们合成的二维 VAR(1)（系数缓慢漂移 + t=200 处一次结构性跳跃）上，卡尔曼滤波估计相对滚动 OLS，把系数追踪的均方误差从 1.62e-01 降到 2.03e-02，降幅 87%；脉冲响应显示，同一单位冲击在 t=100（跳跃前）与 t=300（跳跃后）的传导形状明显不同。诚实地说：当系数确实恒定（纯随机游走式关系）时，TVP 不会凭空造出结构，它只是把真实存在的时变与跳跃精确地描出来。

二维 VAR(1) 合成序列与真实时变系数

1. 为什么常系数 VAR 会误判关系#

VAR(p) 把一串变量联合建模为各自对「全体变量的滞后项」做回归：

$y_t = A_1 y_{t-1} + A_2 y_{t-2} + \cdots + \varepsilon_t$

常系数版本里，所有 $A_i$ 是常数。它的隐含假设是「经济结构稳定」。但现实里：

缓慢漂移：货币政策反应函数随通胀预期微调，因子对风格的相关性随拥挤度渐变。
结构性跳跃：一次危机、一个监管新规、一个流动性事件，能让传导系数在几期内跳变（比如 Copula 相关系数在 2008 年从 0.2 跳到 0.8）。

常系数 VAR 把这两类都压成单一常数，要么滞后（滚动窗口），要么抹平（全样本）。TVP 框架把 $A_t$ 写成状态变量：

$y_t = A_t\,y_{t-1} + \varepsilon_t, \qquad \text{vec}(A_t) = \text{vec}(A_{t-1}) + \eta_t$

系数在「接近上期值」的先验下，用新观测逐期更新——这就是标准状态空间 + 卡尔曼滤波。

2. 卡尔曼滤波：逐期更新时变系数#

把第 $i$ 个方程单独看，回归量 $X_t=[y_{t-1,1},y_{t-1,2}]$ ，待估系数 $\beta_t$ 是状态：

状态方程： $\beta_t = \beta_{t-1} + \eta_t$ （随机游走过程先验）
观测方程： $y_{t,i} = X_t\,\beta_t + \varepsilon_t$

卡尔曼滤波在每一期做「预测→更新」两步，输出系数的滤波估计 $\hat\beta_t$ 及其协方差（即不确定度）。

import numpy as np

def kalman_tvp(y, X, Q=2e-4, R=0.09):
    """对单个方程做 TVP 回归，返回滤波后的状态序列 (len(y), n)。"""
    n = X.shape[1]
    beta_f = np.zeros((len(y), n))
    P = np.eye(n) * 1.0          # 状态协方差
    beta = np.zeros(n)           # 状态均值
    I = np.eye(n)
    for t in range(len(y)):
        x = X[t]
        # 预测步（恒等转移，加过程噪声 Q）
        P = P + Q * I
        pred = x @ beta
        # 更新步
        S = x @ P @ x + R
        K = P @ x / S
        beta = beta + K * (y[t] - pred)
        P = (I - np.outer(K, x)) @ P
        beta_f[t] = beta
    return beta_f

python

$Q$ 控制系数「允许变多快」（越大越灵敏、越抖）， $R$ 是观测噪声方差（由数据波动标定）。二者是模型最关键的超参数。

3. 合成数据 + 系数路径对比#

我们生成 $T=400$ 的二维 VAR(1)：系数随 $t$ 缓慢线性漂移，并在 $t=200$ 处给 $A_{11},A_{22}$ 各加一次跳跃（模拟一次 regime 切换）。

rng = np.random.default_rng(20260721)
T = 400
A_true = np.zeros((T, 2, 2))
for t in range(T):
    j = 0.35 if t >= 200 else 0.0           # t=200 处跳跃
    a00 = 0.55 + 0.0008 * t + j + rng.normal(0, 0.004)
    a01 = -0.20 + 0.0002 * t
    a10 = 0.15 - 0.0003 * t
    a11 = 0.45 - 0.0006 * t + 0.25 * (1 if t >= 200 else 0) + rng.normal(0, 0.004)
    A_true[t] = [[a00, a01], [a10, a11]]

y = np.zeros((T, 2)); y[0] = [0, 0]
for t in range(1, T):
    eps = rng.multivariate_normal([0, 0], [[0.09, 0.02], [0.02, 0.09]])
    y[t] = A_true[t] @ y[t - 1] + eps

X = np.column_stack([y[:-1, 0], y[:-1, 1]])
Ahat1 = kalman_tvp(y[1:, 0], X, Q=2e-4, R=0.09)   # 方程1 时变系数
# 滚动 OLS 基线（窗口 40）
def rolling_ols(y, X, W=40):
    out = np.zeros_like(Ahat1)
    for t in range(len(y)):
        if t < W:
            out[t] = Ahat1[0]; continue
        coef, *_ = np.linalg.lstsq(X[t-W:t], y[t-W:t], rcond=None)
        out[t] = coef
    return out
Roll1 = rolling_ols(y[1:, 0], X, W=40)

python

时变系数：真实 vs 卡尔曼滤波 vs 滚动OLS

上图是第一个方程的系数路径。关键对比一目了然：

卡尔曼滤波（蓝）：紧贴真实路径（金），漂移跟得上、t=200 的跳跃也跳得出来。
滚动 OLS（红）：明显滞后——t=200 的跳跃它被平滑成一段缓慢斜坡，且整体比真实值晚约半个窗长（≈20 期）才到位。

把两套估计对齐真实 $A_t$ ，逐元素算均方误差：

$\text{MSE} = \frac1{T}\sum_t \big\|\hat A_t - A_t^{\text{true}}\big\|_F^2$

结果：TVP-KF 的 MSE = 2.03e-02，滚动 OLS = 1.62e-01，降幅 87%。 这不是因为 TVP「更会拟合噪声」——当系数真恒定时两者都会收敛；差距完全来自 TVP 处理时变与跳跃的能力。

4. 脉冲响应：同一冲击，不同时点形状不同#

VAR 的 $h$ 步脉冲响应就是 $A^h$ 乘冲击向量。系数时变时，响应函数也时变：

def irf(A_t, shock, H=12):
    resp = np.zeros((H, 2)); cur = shock.astype(float)
    resp[0] = cur
    for h in range(1, H):
        cur = A_t @ cur; resp[h] = cur
    return resp

shock = np.array([1.0, 0.0])          # 给变量1一个单位冲击
irf_early = irf(A_true[100], shock)   # 跳跃前
irf_late  = irf(A_true[300], shock)   # 跳跃后

python

脉冲响应：跳跃前后传导形状不同

给变量 1 一个单位冲击，在 t=100（跳跃前）冲击快速衰减、对变量 2 的传导弱；在 t=300（跳跃后）冲击更持久、对变量 2 的溢出更强。这是用常系数 VAR 做 IRF 会直接漏掉的结构——常系数只能给你「一条平均的」响应曲线，而真实经济里，同样的加息冲击在 2006 年和 2020 年的传导路径本就不同。TVP 的 IRF 把「时间维度」还给了脉冲响应。

5. 五类真实陷阱（中高阶）#

过程噪声 Q 是双刃剑：Q 太小，滤波系数僵化、跟不上漂移（退化为近常数）；Q 太大，系数被单期噪声牵着抖、过拟合。Q 没有银弹，必须按「你预期系数最小变化速度」标定，并做敏感性扫描（本文用 Q=2e-4 是合成数据的示范值）。
观测噪声 R 需与数据波动匹配：R 定错会把截距/水平误当成系数变化。实务中 R 用回归残差的截面/时序方差标定，而非拍脑袋。
状态空间可识别性：VAR 同时估所有系数时，状态维度 = $k^2$ （二维就是 4 个），维度随变量数平方爆炸。变量 > 5 时建议配收缩先验（如明尼苏达先验）或降维（因子 TVP-VAR），否则协方差矩阵病态、估计发散。
跳跃 vs 漂移的混淆：本文把跳跃当成「Q 允许下的快速漂移」捕获。但若跳跃幅度极大、Q 又小，滤波会滞后定位。极端结构断点应配专门的变点检测（见本专栏贝叶斯结构断点），TVP 负责「连续时变」、变点检测负责「离散跳」。
前视偏差：所有系数估计在 $t$ 时刻只能用 $t$ 及之前的数据。任何把 $t$ 之后样本喂进卡尔曼滤波或滚动窗口的写法，都会让「完美追踪」变成回测幻象。

6. 实战落地点：从系数路径到交易信号#

TVP-VAR 不是论文玩具，它直接对应三条实战监控线，而且和本专栏其他工具天然互补：

时变 beta 暴露监控：把组合收益对风格因子做 TVP 回归，beta 不再是「一个历史平均数」，而是一条实时更新的路径。当某因子 beta 在滤波不确定度之外跌破阈值（比如从 0.6 掉到 0.1），往往比「看净值回撤」更早预警「这个因子可能失效了」——因为 beta 坍塌通常发生在净值崩坏之前。
跨资产相关实时图：把相关性矩阵当成一组 TVP 系数逐期更新，危机时相关飙升会比「用 60 天滚动窗口算相关」早约半个窗长被捕捉到，给风险平价、对冲比率调整留出反应时间。
IRF 情景分析：用当前时点的 $A_t$ 算脉冲响应，回答「如果美联储现在加息 25bp，我的组合 3 个月后大概承受多大传导」——这是常系数 VAR 给不出、却又是风控最该问的问题。
与结构断点工具分工：TVP 负责「连续、平滑的时变」，本专栏的贝叶斯结构断点检测负责「离散、剧烈的跳」。生产环境里两者并联：TVP 给出平滑 beta 路径，断点检测器在 TVP 残差上找突变，互不重叠又互相印证。

需要提醒的是：TVP 估出的系数自带不确定度（卡尔曼协方差 $P_t$ ）。实战里一定要把不确定度画出来——当 $P_t$ 很大（数据稀薄、刚经历跳变）时，所谓「beta 信号」可能只是估计噪声，不应触发交易。把不确定度当滤波器，是 TVP 从「能算」到「能用」的关键一步。

7. 完整 Python 代码#

下面是与本文全部数字、配图一一对应的端到端复现脚本（自洽合成数据，仅演示方法；真实落地需替换为你自己的收益率/宏观序列与标定过的 Q、R）：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams["font.sans-serif"] = ["Arial Unicode MS", "PingFang SC", "SimHei", "DejaVu Sans"]
rcParams["axes.unicode_minus"] = False

def kalman_tvp(y, X, Q=2e-4, R=0.09):
    n = X.shape[1]; beta_f = np.zeros((len(y), n)); P = np.eye(n) * 1.0
    beta = np.zeros(n); I = np.eye(n)
    for t in range(len(y)):
        x = X[t]; P = P + Q * I; pred = x @ beta
        S = x @ P @ x + R; K = P @ x / S
        beta = beta + K * (y[t] - pred); P = (I - np.outer(K, x)) @ P
        beta_f[t] = beta
    return beta_f

def rolling_ols(y, X, W=40):
    out = np.zeros_like(Ahat1 if False else np.zeros((len(y), X.shape[1])))
    for t in range(len(y)):
        if t < W: out[t] = 0; continue
        coef, *_ = np.linalg.lstsq(X[t-W:t], y[t-W:t], rcond=None); out[t] = coef
    return out

def irf(A_t, shock, H=12):
    resp = np.zeros((H, 2)); cur = shock.astype(float); resp[0] = cur
    for h in range(1, H): cur = A_t @ cur; resp[h] = cur
    return resp

# 合成
rng = np.random.default_rng(20260721); T = 400
A_true = np.zeros((T, 2, 2))
for t in range(T):
    j = 0.35 if t >= 200 else 0.0
    a00 = 0.55 + 0.0008*t + j + rng.normal(0,0.004)
    a01 = -0.20 + 0.0002*t
    a10 = 0.15 - 0.0003*t
    a11 = 0.45 - 0.0006*t + 0.25*(1 if t>=200 else 0) + rng.normal(0,0.004)
    A_true[t] = [[a00,a01],[a10,a11]]
y = np.zeros((T,2)); y[0]=[0,0]
for t in range(1,T):
    eps = rng.multivariate_normal([0,0], [[0.09,0.02],[0.02,0.09]])
    y[t] = A_true[t] @ y[t-1] + eps

X = np.column_stack([y[:-1,0], y[:-1,1]])
Ahat1 = kalman_tvp(y[1:,0], X); Roll1 = rolling_ols(y[1:,0], X, W=40)
Ahat2 = kalman_tvp(y[1:,1], X); Roll2 = rolling_ols(y[1:,1], X, W=40)

# 系数追踪 MSE 对比（对齐真实 A_t）
At = A_true[1:]
Ahat_stack = np.column_stack([Ahat1, Ahat2])
Roll_stack = np.column_stack([Roll1, Roll2])
Atrue_stack = np.column_stack([At[:,0,:], At[:,1,:]])
mse_tvp = np.mean((Ahat_stack - Atrue_stack)**2)
mse_roll = np.mean((Roll_stack - Atrue_stack)**2)
print(f"MSE TVP={mse_tvp:.4e}  Roll={mse_roll:.4e}  降幅={100*(1-mse_tvp/mse_roll):.0f}%")

# 脉冲响应
shock = np.array([1.0, 0.0])
irf_early = irf(A_true[100], shock); irf_late = irf(A_true[300], shock)
print("IRF@100 末步:", irf_early[-1].round(3), " IRF@300 末步:", irf_late[-1].round(3))

python

真实落地提示：把 $y_t$ 换成你的资产收益率/宏观因子、把 $Q,R$ 用样本残差标定、变量多时加明尼苏达收缩先验；输出可直接接「时变 beta 暴露监控」「跨资产相关实时图」「regime 切换预警」等实战模块。