跳到正文

VAR 模型:向量自回归

Sims 用'无理论'的 VAR 掀翻了联立方程宏观模型。但你对 VAR 提出的每一个有意思的问题——脉冲响应、方差分解、政策分析——都会把识别假设从后门请回来。推导 VAR、Granger 因果、IRF、FEVD 和结构识别。

第 12 篇 · 共 12 篇 约 13 分钟

Sims (1980) 之前,宏观经济学家建巨型联立方程模型,用”理论”给识别约束找理由。Sims 指出这些理由基本站不住:约束不可检验、模型拟合差、政策结论依赖没人能辩护的假设。他的解法是无约束向量自回归——让每个变量对所有变量的滞后回归,不施加任何”不可信的约束”。这个提议简单到迫使整个领域反思建模方式。

但这篇文章想让你带走的不只是 VAR 怎么推,还有一条贯穿始终的线:“无理论”只是把理论从估计端挪到了解读端。估计一个 VAR 确实不需要经济学;可你一旦问它有意思的问题——冲击怎么传导、谁驱动谁、政策改变会怎样——识别假设就会从后门回来,而且每一步都得你自己选。看清这一点,就看清了从 Cholesky 排序到符号约束的整个谱系在争什么。参考 Hamilton (1994, 第 11 章)、Lütkepohl (2005)、Tsay (2010, 第 4 章)。

基本 VAR

KK 变量 pp 阶 VAR:

yt=c+A1yt−1+A2yt−2+⋯+Apyt−p+ut\mathbf{y}_t = \mathbf{c} + \mathbf{A}_1 \mathbf{y}_{t-1} + \mathbf{A}_2 \mathbf{y}_{t-2} + \cdots + \mathbf{A}_p \mathbf{y}_{t-p} + \mathbf{u}_t

其中 yt=(y1t,y2t,…,yKt)′\mathbf{y}_t = (y_{1t}, y_{2t}, \ldots, y_{Kt})' 是 KK 维向量,c\mathbf{c} 是截距向量,Ai\mathbf{A}_i 是 K×KK \times K 系数矩阵,ut=(u1t,…,uKt)′\mathbf{u}_t = (u_{1t}, \ldots, u_{Kt})' 是创新向量,满足:

E[ut]=0,E[utut′]=Σ,E[utus′]=0 for t≠sE[\mathbf{u}_t] = \mathbf{0}, \quad E[\mathbf{u}_t \mathbf{u}_t'] = \boldsymbol{\Sigma}, \quad E[\mathbf{u}_t \mathbf{u}_s'] = \mathbf{0} \text{ for } t \neq s

注意 Σ\boldsymbol{\Sigma} 不要求对角——同期的创新可以相关,这个”可以相关”是后面一半故事的起点。

平稳性条件:det⁡(IK−A1z−⋯−Apzp)=0\det(\mathbf{I}_K - \mathbf{A}_1 z - \cdots - \mathbf{A}_p z^p) = 0 的所有根在单位圆外。伴随形式让这一点变透明。定义 Kp×1Kp \times 1 状态向量:

ξt=(yt′,yt−1′,…,yt−p+1′)′\boldsymbol{\xi}_t = (\mathbf{y}_t', \mathbf{y}_{t-1}', \ldots, \mathbf{y}_{t-p+1}')'

那么:

ξt=Fξt−1+Gut\boldsymbol{\xi}_t = \mathbf{F} \boldsymbol{\xi}_{t-1} + \mathbf{G} \mathbf{u}_t

其中

F=(A1A2⋯Ap−1ApIK0⋯000IK⋯00⋮⋱⋮00⋯IK0),G=(IK0⋮0)\mathbf{F} = \begin{pmatrix} \mathbf{A}_1 & \mathbf{A}_2 & \cdots & \mathbf{A}_{p-1} & \mathbf{A}_p \\ \mathbf{I}_K & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_K & \cdots & \mathbf{0} & \mathbf{0} \\ \vdots & & \ddots & & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{I}_K & \mathbf{0} \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \mathbf{I}_K \\ \mathbf{0} \\ \vdots \\ \mathbf{0} \end{pmatrix}

任何 VAR(pp) 都被压成一个大 VAR(1),平稳当且仅当 F\mathbf{F} 所有特征值的模小于 1。

平稳的 VAR 可以求逆成 VMA(∞\infty) 表示:

yt=μ+∑j=0∞Ψjut−j\mathbf{y}_t = \boldsymbol{\mu} + \sum_{j=0}^\infty \boldsymbol{\Psi}_j \mathbf{u}_{t-j}

其中 μ=(I−A1−⋯−Ap)−1c\boldsymbol{\mu} = (\mathbf{I} - \mathbf{A}_1 - \cdots - \mathbf{A}_p)^{-1} \mathbf{c},Ψj\boldsymbol{\Psi}_j 满足递推 Ψj=A1Ψj−1+⋯+ApΨj−p\boldsymbol{\Psi}_j = \mathbf{A}_1 \boldsymbol{\Psi}_{j-1} + \cdots + \mathbf{A}_p \boldsymbol{\Psi}_{j-p},初始条件 Ψ0=IK\boldsymbol{\Psi}_0 = \mathbf{I}_K、Ψj=0\boldsymbol{\Psi}_j = \mathbf{0} 对 j<0j < 0。这些移动平均矩阵记录 jj 期前的冲击怎么传到当下——先记住它们,脉冲响应整个建在上面。

估计:便宜得可疑

估计是 VAR 最省心的部分。把 KK 个方程摆开,逐方程 OLS 就够了,而且逐方程 OLS 等价于系统 GLS——因为每个方程的回归子完全相同(都是滞后的 y\mathbf{y} 加截距),跨方程的误差相关不带来效率损失。对每个方程 ii,把 yity_{it} 对常数和 yt−1,…,yt−p\mathbf{y}_{t-1}, \ldots, \mathbf{y}_{t-p} 回归,得到 A^j\hat{\mathbf{A}}_j 和残差 u^t\hat{\mathbf{u}}_t;协方差矩阵用残差的样本协方差,带自由度修正(每方程 Kp+1Kp + 1 个回归子):

Σ^=1T−Kp−1∑t=p+1Tu^tu^t′\hat{\boldsymbol{\Sigma}} = \frac{1}{T - Kp - 1} \sum_{t=p+1}^T \hat{\mathbf{u}}_t \hat{\mathbf{u}}_t'

标准正则条件下估计量一致且渐近正态(Hamilton 1994, 定理 11.4 附近的结果)。堆叠系数向量 β^=vec([c^,A^1,…,A^p])\hat{\boldsymbol{\beta}} = \text{vec}([\hat{\mathbf{c}}, \hat{\mathbf{A}}_1, \ldots, \hat{\mathbf{A}}_p]) 的渐近协方差:

Var(β^)=1TΓ−1⊗Σ\text{Var}(\hat{\boldsymbol{\beta}}) = \frac{1}{T} \boldsymbol{\Gamma}^{-1} \otimes \boldsymbol{\Sigma}

其中 Γ=E[xtxt′]\boldsymbol{\Gamma} = E[\mathbf{x}_t \mathbf{x}_t'] 是回归子(截距与各期滞后 y\mathbf{y} 堆叠成的 xt\mathbf{x}_t)的二阶矩矩阵,由 y\mathbf{y} 各滞后的自协方差拼成。实际上这些软件都替你做了——Python 里 statsmodels.tsa.api.VAR,R 里 vars 包。

便宜的地方在参数量上收费。每个方程 Kp+1Kp + 1 个参数,系统一共 K2p+KK^2 p + K 个,外加协方差矩阵的 K(K+1)/2K(K+1)/2 个。K=5K = 5、p=4p = 4 就是一百多个系数,自由度立刻紧张。所以滞后阶选择不是走过场:AIC、BIC、HQ 都能用,BIC 倾向于比 AIC 更小的 pp;Sims (1980) 的建议是对滞后慷慨一点以避免动态设定偏差,渐近理论 Lütkepohl (2005, 第 4 章) 有综述。参数灾难更彻底的解法是贝叶斯收缩,后面讲。

Granger 因果:VAR 里最容易的检验,最容易被误读的词

变量 y1ty_{1t} Granger 引起 y2ty_{2t},是说在 y2ty_{2t} 自身滞后之外,y1ty_{1t} 的滞后对预测 y2ty_{2t} 有增量帮助。在 VAR 里这是个直接的系数检验。写出 y2ty_{2t} 的方程:

y2t=c2+∑j=1p(β2,1jy1,t−j+β2,2jy2,t−j)+u2ty_{2t} = c_2 + \sum_{j=1}^p (\beta_{2,1j} y_{1,t-j} + \beta_{2,2j} y_{2,t-j}) + u_{2t}

零假设是 y1ty_{1t} 的滞后全不进入:H0:β2,1,1=⋯=β2,1,p=0H_0: \beta_{2,1,1} = \cdots = \beta_{2,1,p} = 0,用 Wald 或 F 检验。等价地比较无约束和有约束模型的似然:

LR=2(ℓunrestricted−ℓrestricted)→dχp2LR = 2(\ell_{\text{unrestricted}} - \ell_{\text{restricted}}) \xrightarrow{d} \chi^2_p

LRLR 大就拒绝。对每个有序对跑一遍,得到完整的 Granger 因果矩阵。另有一个单独的问题叫瞬时因果:控制动态之后,u1tu_{1t} 和 u2tu_{2t} 的同期相关是否为零——看 Σ^\hat{\boldsymbol{\Sigma}} 的非对角元,零假设下样本协方差渐近正态、方差量级 1/T1/T。

现在说那个词。“Granger 因果”在任何哲学意义上都不是因果,它的全部含义是”对预测有增量信息”。一个变量可以 Granger 引起另一个而不是它的原因(两者被第三个变量共同驱动);一个变量可以是原因却不 Granger 引起结果(效应是瞬时的,被同期相关吃掉了)。这个区分不是学究气——它直接决定检验结果能拿来说什么。实证里的常见问法:货币供给是否 Granger 引起通胀(宏观数据里通常是);收益率曲线是否 Granger 引起衰退(期限利差的预测能力是文献里的稳定发现);油价对股票收益(有争议);成交量对波动率(常见的肯定结果)。每一条都是预测性陈述,不是机制陈述。

检验本身也有几个会咬人的地方。漏掉共同原因,两个变量会显得互相 Granger 引起。变量 I(1) 且不协整时,检验的分布非标准、容易拒绝过多——回到平稳性那块地基。关系还会随样本期漂移:1960-1990 成立的 1990-2020 未必成立。最后,标准检验是线性的,非线性的可预测性它抓不到。

脉冲响应:识别假设从这里进场

IRF 回答的问题:今天给系统一个冲击,所有变量未来 hh 期怎么走?答案就藏在 VMA 矩阵里。在时刻 0 给创新一个冲击 δ\boldsymbol{\delta}、其余期为零:

E[yh∣u0=δ,u−1=u−2=⋯=0]−E[yh∣u0=0,…]=ΨhδE[\mathbf{y}_h \mid \mathbf{u}_0 = \boldsymbol{\delta}, \mathbf{u}_{-1} = \mathbf{u}_{-2} = \cdots = \mathbf{0}] - E[\mathbf{y}_h \mid \mathbf{u}_0 = \mathbf{0}, \ldots] = \boldsymbol{\Psi}_h \boldsymbol{\delta}

对第 ii 个创新的 hh 期响应就是 Ψh\boldsymbol{\Psi}_h 的第 ii 列。到这里还不需要任何经济学。

问题马上来了:Σ\boldsymbol{\Sigma} 有非对角元,创新在方程之间同期相关。“只冲击 u1tu_{1t}、其他保持为零”这个思想实验在数据里不存在——它们通常一起动。要把冲击掰开,得先让它们正交。标准做法是 Cholesky 分解:

Σ=PP′\boldsymbol{\Sigma} = \mathbf{P} \mathbf{P}'

P\mathbf{P} 下三角、对角元为正。定义正交化冲击 et=P−1ut\mathbf{e}_t = \mathbf{P}^{-1} \mathbf{u}_t,它们互不相关、方差为单位:

E[etet′]=P−1Σ(P−1)′=P−1PP′(P′)−1=IKE[\mathbf{e}_t \mathbf{e}_t'] = \mathbf{P}^{-1} \boldsymbol{\Sigma} (\mathbf{P}^{-1})' = \mathbf{P}^{-1} \mathbf{P} \mathbf{P}' (\mathbf{P}')^{-1} = \mathbf{I}_K

对第 jj 个正交冲击的 IRF 是 ΨhP\boldsymbol{\Psi}_h \mathbf{P} 的第 jj 列。

注意刚刚发生了什么:Cholesky 不是中性的数学变换,它强加了递归结构——变量 1 的冲击同期影响所有变量,变量 2 的冲击同期影响变量 2 到 KK 但不影响 1,依此类推。这是一个因果排序:变量 1 同期最”外生”,变量 KK 最内生。这就是开头说的那件事——“无理论”的 VAR 在你画第一张正交化 IRF 图的瞬间就重新需要理论了,理论就藏在排序里。惯例是按反应速度排:慢变量(实际 GDP、资本存量)在前,快变量(利率、货币总量)在后。但惯例不等于真理:把货币政策放在产出前面和后面,会得到不同的答案。所以拿到 IRF 先做一件事——换几个排序重跑。结果定性不变,识别还算结实;定性翻转,你报告的是排序不是数据。

预测误差方差分解

FEVD 问的是长期账目:变量 ii 的 hh 步预测误差方差里,每个正交冲击各占多少?hh 步预测误差:

yT+h−y^T+h∣T=∑j=0h−1ΨjuT+h−j\mathbf{y}_{T+h} - \hat{\mathbf{y}}_{T+h|T} = \sum_{j=0}^{h-1} \boldsymbol{\Psi}_j \mathbf{u}_{T+h-j}

MSE 矩阵:

MSE(h)=∑j=0h−1ΨjΣΨj′\text{MSE}(h) = \sum_{j=0}^{h-1} \boldsymbol{\Psi}_j \boldsymbol{\Sigma} \boldsymbol{\Psi}_j'

代入 Σ=PP′\boldsymbol{\Sigma} = \mathbf{P}\mathbf{P}',用正交冲击改写:

yT+h−y^T+h∣T=∑j=0h−1ΨjPeT+h−j\mathbf{y}_{T+h} - \hat{\mathbf{y}}_{T+h|T} = \sum_{j=0}^{h-1} \boldsymbol{\Psi}_j \mathbf{P} \mathbf{e}_{T+h-j}

冲击 jj 对变量 ii 在 hh 期 MSE 的占比:

ωij(h)=∑k=0h−1[(ΨkP)ij]2MSEii(h)\omega_{ij}(h) = \frac{\sum_{k=0}^{h-1} \left[ (\boldsymbol{\Psi}_k \mathbf{P})_{ij} \right]^2}{\text{MSE}_{ii}(h)}

对 jj 求和等于 1。解读:ω1,2(24)=0.6\omega_{1,2}(24) = 0.6 意味着变量 1 的 24 步预测误差方差有 60% 归到变量 2 的冲击——这就是”变量 2 长期驱动变量 1”这类说法的技术底座。当然,它继承了 IRF 的一切识别依赖:换个排序,分解就换一副面孔。

结构 VAR:把识别摆上台面

Cholesky 只是识别谱系的起点。把问题正式化:结构 VAR 假设存在经济上有含义的结构冲击 εt\boldsymbol{\varepsilon}_t,结构形式:

A0yt=c+A1∗yt−1+⋯+Ap∗yt−p+Bεt\mathbf{A}_0 \mathbf{y}_t = \mathbf{c} + \mathbf{A}_1^* \mathbf{y}_{t-1} + \cdots + \mathbf{A}_p^* \mathbf{y}_{t-p} + \mathbf{B} \boldsymbol{\varepsilon}_t

其中 E[εtεt′]=IKE[\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_t'] = \mathbf{I}_K(正交、单位方差),A0\mathbf{A}_0 捕捉同期交互,B\mathbf{B} 把结构冲击映射到简化型误差。两边乘 A0−1\mathbf{A}_0^{-1} 回到简化型:

yt=c+A1yt−1+⋯+ut\mathbf{y}_t = \mathbf{c} + \mathbf{A}_1 \mathbf{y}_{t-1} + \cdots + \mathbf{u}_t

其中 Ai=A0−1Ai∗\mathbf{A}_i = \mathbf{A}_0^{-1} \mathbf{A}_i^*,ut=A0−1Bεt\mathbf{u}_t = \mathbf{A}_0^{-1} \mathbf{B} \boldsymbol{\varepsilon}_t,于是 Σ=A0−1BB′(A0−1)′\boldsymbol{\Sigma} = \mathbf{A}_0^{-1} \mathbf{B} \mathbf{B}' (\mathbf{A}_0^{-1})'。识别的算术:待定的映射有 K2K^2 个自由元,而 Σ\boldsymbol{\Sigma} 只提供 K(K+1)/2K(K+1)/2 个可识别的量,缺口 K2−K(K+1)/2=K(K−1)/2K^2 - K(K+1)/2 = K(K-1)/2 个约束——这些约束数据给不了,只能理论给。往哪里要,分出了三大流派。

递归识别(Sims, 1980):A0\mathbf{A}_0 下三角、对角线为 1,B\mathbf{B} 对角。这正是 Cholesky——下三角结构天然提供 K(K−1)/2K(K-1)/2 个约束,前一节的排序问题就是在选这个三角形。

长期约束(Blanchard & Quah, 1989):不约束同期,约束长期。识别假设是需求冲击对产出没有长期效应,而供给冲击可以有永久效应——约束落在长期乘子 ∑j=0∞ΨjA0−1B\sum_{j=0}^\infty \boldsymbol{\Psi}_j \mathbf{A}_0^{-1} \mathbf{B} 上。适合理论对”什么是暂时的、什么是永久的”有明确话可说的场景。

符号约束(Uhlig, 2005):最弱的一档。只要求某些冲击对某些变量在某些期取某些符号——比如紧缩性货币冲击不该让价格上涨。代价是从点识别退到集合识别:与约束相容的是一个冲击区域,贝叶斯方法从中抽样。当前宏观实证里这已是标配。

三派没有免费的:递归识别的排序、长期约束的理论、符号约束的先验,都是你替数据做的选择。VAR 把这些选择从联立方程时代的暗处搬到了明处——这大概才是 Sims 革命真正的遗产。

水平里有单位根:去 VECM

以上全部假设平稳。如果水平 yt\mathbf{y}_t 非平稳、但某个线性组合 β′yt\boldsymbol{\beta}' \mathbf{y}_t 平稳,那么 yt\mathbf{y}_t 是协整的,对水平估 VAR 是设定错误,有限样本里标准误也会失真。全差分能救平稳性,但把长期信息扔了;正确的框架是向量误差修正模型:

Δyt=α(β′yt−1−μ)+∑j=1p−1ΓjΔyt−j+ut\Delta \mathbf{y}_t = \boldsymbol{\alpha}(\boldsymbol{\beta}' \mathbf{y}_{t-1} - \boldsymbol{\mu}) + \sum_{j=1}^{p-1} \boldsymbol{\Gamma}_j \Delta \mathbf{y}_{t-j} + \mathbf{u}_t

β′yt−1\boldsymbol{\beta}' \mathbf{y}_{t-1} 是协整关系,α\boldsymbol{\alpha} 是调整速度。协整值得两篇文章的篇幅,紧接这篇之后。

贝叶斯 VAR:参数灾难的解法

回到估计一节埋的雷:K2p+KK^2 p + K 个参数。K=10K = 10、p=4p = 4 就是 410 个,而宏观数据 T≈200T \approx 200——无约束估计根本立不住。Litterman (1986) 和 Sims & Zha (1998) 的贝叶斯 VAR 用先验把系数往合理的默认值收缩。主力是 Minnesota 先验,它的默认信念是”每个变量近似自己的随机游走”:自身一阶滞后的系数先验均值为 1,其余系数先验均值为 0;先验方差随滞后阶增大而收缩(远处的滞后先验上更接近零),交叉变量的滞后再按两个变量残差方差之比缩放,让先验的松紧在不同量纲的变量之间可比。后验均值于是成为 OLS 估计和先验均值的加权平均,权重由先验相对数据的精度决定。效果实打实:BVAR 的样本外预测好,尤其在 TT 相对参数不富裕时——它绕开了无约束 VAR 的过拟合。

预测与落地

VAR 的 hh 期预测递推地做:

y^T+h=c^+∑j=1pA^jy^T+h−j\hat{\mathbf{y}}_{T+h} = \hat{\mathbf{c}} + \sum_{j=1}^p \hat{\mathbf{A}}_j \hat{\mathbf{y}}_{T+h-j}

h≤ph \leq p 时右边混着观测值和已算出的预测,h>ph > p 之后全是预测。预测误差:

yT+h−y^T+h=∑j=0h−1ΨjuT+h−j\mathbf{y}_{T+h} - \hat{\mathbf{y}}_{T+h} = \sum_{j=0}^{h-1} \boldsymbol{\Psi}_j \mathbf{u}_{T+h-j}

MSE 矩阵还是 MSE(h)=∑j=0h−1ΨjΣΨj′\text{MSE}(h) = \sum_{j=0}^{h-1} \boldsymbol{\Psi}_j \boldsymbol{\Sigma} \boldsymbol{\Psi}_j',高斯误差下预测区间:

y^T+h±1.96⋅diag(MSE(h))\hat{\mathbf{y}}_{T+h} \pm 1.96 \cdot \sqrt{\text{diag}(\text{MSE}(h))}

代码上一切都是现成的:

import pandas as pd
from statsmodels.tsa.api import VAR


model = VAR(y)
lag_order = model.select_order(maxlags=8).aic
print("Selected lag order:", lag_order)

results = model.fit(lag_order)
print(results.summary())

# Granger 因果
print(results.test_causality('gdp', caused=['inflation', 'interest_rate'], kind='f'))

# IRF 和 FEVD
irf = results.irf(periods=20)
irf.plot(orth=True)

fevd = results.fevd(periods=20)
fevd.plot()

用在哪。宏观预测:BVAR 是 GDP、通胀、失业 nowcasting 和短期预测的标配,各国央行的研究部门用它做政策分析。货币政策分析:对产出、价格、利率、货币估 VAR,识别一个紧缩性冲击(比如经符号约束),追踪产出和价格的反应——货币政策冲击是现代实证货币经济学的基石。财政政策同理,把政府支出和税收当冲击,乘数估计(政府多花 1 美元,产出反应多少)就出自财政 SVAR。金融这边:对收益率曲线因子(水平、斜率、曲率)跑 VAR 研究期限结构动态,对 Fama-French 因子跑 VAR 研究风险溢价,对收益率和宏观变量跑 VAR 研究收益可预测性。混合频率 VAR(MF-VAR)用周度、月度数据实时更新对当季度变量的估计,这就是 nowcasting 的机器。

不论用在哪,交付结论前把这篇的主线再过一遍:估计端你没用理论,解读端你用了哪些?滞后阶是谁选的、排序换过没有、水平里有没有单位根、“因果”两个字指的是预测还是机制。VAR 的诚实之处在于这些问题都有明确的位置可查——前提是你真的去查。

参考文献

  • Blanchard, O. J., & Quah, D. (1989). The dynamic effects of aggregate demand and aggregate supply disturbances. American Economic Review, 79(4), 655-673.
  • Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Ch. 11.
  • Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer.
  • Sims, C. A. (1980). Macroeconomics and reality. Econometrica, 48(1), 1-48.
  • Sims, C. A., & Zha, T. (1998). Bayesian methods for dynamic multivariate models. International Economic Review, 39(4), 949-968.
  • Stock, J. H., & Watson, M. W. (2001). Vector autoregressions. Journal of Economic Perspectives, 15(4), 101-115.
  • Tsay, R. S. (2010). Analysis of Financial Time Series (3rd ed.). Wiley. Ch. 4.
  • Uhlig, H. (2005). What are the effects of monetary policy on output? Results from an agnostic identification procedure. Journal of Monetary Economics, 52(2), 381-419.
  • Litterman, R. B. (1986). Forecasting with Bayesian vector autoregressions. Five years of experience. Journal of Business & Economic Statistics, 4(1), 25-38.