Notebook

Granger因果检验¶

模型介绍¶

经济学中常常需要确定因果关系究竟是从 x 到 y，还是从 y 到 x，亦或双向因果关系。格兰杰提出的检验方法基于以下思想。如果 x 是 y 的因，但 y 不是 x 的因，则 x 的过去值可以帮助预测 y 的未来值，但 y 的过去值却不能帮助预测 x 的未来值。考虑以下时间序列模型：

$y_{t}=\gamma+\sum_{m=1}^{p} \alpha_{m} y_{t-m}+\sum_{m=1}^{p} \beta_{m} x_{t-m}+\varepsilon_{t}\tag{1}$

其中，滞后阶数 $p$ 可根据“信息准则”或“由大到小的序贯 $t$ 规则”来确定。检验原假设“ $H_{0}: \beta_{1}=$ $=\beta_{p}=0 ",$ 即 $x$ 的过去值对预测 $y$ 的未来值没有帮助。如果拒绝 $H_{0},$ 则称 $x$ 是 $y$ 的“格兰杰因"（ Granger cause）。将以上回归模型中 x 与 y 的位置互换，则可以检验 $y$ 是否为 $x$ 的格兰杰因。在实际操作中，常将（ $x, y$ ）构成一个二元 VAR 系统，然后在 VAR 的框架进行格兰杰因果检验。

需要指出的是，格兰杰因果关系并非真正意义上的因果关系。它充其量只是一种动态相关关系，表明的是一个变量是否对另一变量有"预测能力"（ predictability）。从某种意义上来说，它顶多是因果关系的必要条件（如果不考虑非线性的因果关系）。另外，格兰杰因果关系也可能由第三个变量所引起。在经济学的实证研究中，由于通常不可能进行"控制实验"，能够最有说服力地说明因果关系的当属随机实验与自然实验。

另外，格兰杰因果检验仅适用于平稳序列，或者有协整关系的单位根过程。对于不存在协整关系的单位根变量，则只能先差分，得到平稳序列后再进行格兰杰因果检验。

原理讲解¶

wald 检验介绍¶

自回归分布滞后模型做 wald 检验

对于线性同归模型，检验原假设 $H_{0}: \boldsymbol{\beta}=\boldsymbol{\beta}_{0}$ 其中 $\boldsymbol{\beta}_{K \times 1}$ 为未知参数 $, \boldsymbol{\beta}_{0}$ 已知，共有 $K$ 个约束。

沃尔德检验：通过研究 $\boldsymbol{\beta}$ 的无约束估计量 $\hat{\boldsymbol{\beta}}_{v}$ 与 $\boldsymbol{\beta}_{0}$ 的距离来进行检验。其基本思想是，如果 $H_{0}$ 正确，则 $\left(\hat{\boldsymbol{\beta}}_{v}-\boldsymbol{\beta}_{0}\right)$ 的绝对值不应该很大。沃尔德统计量为：

$W \equiv\left(\hat{\boldsymbol{\beta}}_{v}-\boldsymbol{\beta}_{0}\right)^{\prime}\left[\operatorname{Var}\left(\hat{\boldsymbol{\beta}}_{v}\right)\right]^{-1}\left(\hat{\boldsymbol{\beta}}_{v}-\boldsymbol{\beta}_{0}\right) \stackrel{d}{\longrightarrow} \chi^{2}(K)\tag{2}$

其中， $K$ 为约束条件的个数（即解释变量的个数）。

Granger 因果检验¶

$H_{0}: \quad \beta_{1}=\beta_{2}=0$

$H_{0}:\left(\begin{array}{c}\beta_{1} \\ \beta_{2}\end{array}\right)=\underbrace{\left(\begin{array}{cccc}0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{array}\right)}_{\mathrm{R}}\left(\begin{array}{c}\gamma \\ \alpha_{1} \\ \beta_{1} \\ \beta_{2}\end{array}\right)=\left(\begin{array}{c}0 \\ 0\end{array}\right)$

$H_0\text{：}\underset{m\times K}{\underbrace{R}}\underset{K\times 1}{\underbrace{\beta }}=\underset{m\times 1}{\underbrace{r}}$

根据沃尔德检验原理，由于 $\hat{\boldsymbol{\beta}}$ 是 $\beta$ 的估计量，故如果 $H_{0}$ 成立, 则 $(\boldsymbol{R} \hat{\boldsymbol{\beta}}-\boldsymbol{r})$ 应比较接近 0 向量

这种接近程度可用其二次型来衡量,比如

$(R \hat{\beta}-r)^{\prime}[\operatorname{Var}(R \hat{\beta}-r)]^{-1}(R \hat{\beta}-r)\tag{3}$

其中, $(\boldsymbol{R} \hat{\boldsymbol{\beta}}-\boldsymbol{r})$ 的协方差矩阵可写为

$\begin{aligned} \text{Var}\left( \boldsymbol{R\hat{\beta}}-\boldsymbol{r} \right) &=\text{Var}\left( \boldsymbol{R\hat{\beta}} \right) \,\, \left( \,\,\text{去掉常数，方差不变} \right)\\ &=\boldsymbol{R}\text{Var}\left( \boldsymbol{\hat{\beta}} \right) \boldsymbol{R}^{\prime}\,\, \left( \,\,\text{夹心估计量的公式} \right)\\ &=\sigma ^2\boldsymbol{R}\left( \boldsymbol{X}^{\prime}\boldsymbol{X} \right) ^{-1}\boldsymbol{R}^{\prime}\,\,\left( \,\,\text{因}为 \text{Var}\left( \boldsymbol{\hat{\beta}} \right) =\sigma ^2\left( \boldsymbol{X}^{\prime}\boldsymbol{X} \right) ^{-1} \right)\\ \end{aligned}$

其中， $\sigma^{2}$ 可由 $s^{2}$ 来估计。

In [14]:

import re
import pandas as pd
import statsmodels.api as sm

data = pd.read_excel('../数据/上证指数与沪深300.xlsx')
Y = data['hs300']
X = data['sz']

def lag_list(Y, X, p=1, q=1):
    '''
    待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q)
    获取自回归分布滞后模型的估计向量

    Parameters
    ----------
    Y : 被估计变量
    X : 估计变量
    p : ADL 模型 Y 的滞后阶数，标量默认为1
    q : ADL 模型 X 的滞后阶数，标量默认为1

    Returns
    -------
    ADLy : ADL 模型被解释变量
    ADLx : ADL 模型解释变量

    '''
    ADLx = pd.DataFrame()
    T = len(Y)
    ADLy = list(Y[max(p, q):T])
    for i in range(1, p+1):
        name = f'y_{i}'
        ADLx[name] = list(Y[max(p, q)-i:T-i])
    for i in range(1, q+1):
        name = f'x_{i}'
        ADLx[name] = list(X[max(p, q)-i:T-i])
    return ADLy, ADLx

p = 2
q = 2
ADLy, ADLx = lag_list(Y, X ,p, q)
ADLx = sm.add_constant(ADLx)
mod = sm.OLS(ADLy, ADLx)
res = mod.fit()
# res.summary()

wald = ''
for i, value in enumerate(ADLx):
    if i > p:
        wald += f'{value}='
wald = wald + '0'
wald = str(res.f_test(wald))
walds = float(re.findall('array\(\[\[(.*?)\]\]', wald)[0])
wald

Out[14]:

'<F test: F=array([[7.30924536]]), p=0.000751239087419008, df_denom=453, df_num=2>'

matlab实现¶

function [beta,wald,wald_P,AIC,p,q,cov_mat] = Granger(Y,X)

%待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q)
% 获取Granger因果检验wald值与p值
%  X   - 解释变量，列向量
%  Y   - 被解释变量，列向量
%  p   - ADL模型Y的滞后阶数,标量
%  q   - ADL模型X的滞后阶数,标量
%  wald - wald统计量
%  wald_P - wald统计量的P值
%  AIC - 最小AIC值

% 1.计算各滞后阶的AIC值
AIC = zeros(5,5);
for p = 1:5
    for q = 1:5
        [ADLy,ADLx] = ADLxx(Y,X,p,q);
        [~,~,resid] = regress(ADLy,ADLx);
        T = length(ADLy);
        AIC(p,q) = log(resid'*resid/T)+2*(p+q+1)/T;
    end
end

% 2.找到最小的AIC值所对应的阶数
[p,q] = find(AIC == min(min(AIC)));
[ADLy,ADLx] = ADLxx(Y,X,p,q);
[beta,~,resid] = regress(ADLy,ADLx);

% 3.计算wald统计量
R = [zeros(q,1+p),eye(q)];
% 3.1 计算协方差矩阵
T = length(ADLy);
nu = T-(1+p+q);
siga2 = sum(resid.^2)/nu;
cov_mat = siga2*inv(ADLx'*ADLx);
Rcov_mat = siga2*R*inv(ADLx'*ADLx)*R';
% 3.2 计算wald
wald = (R*beta)'*inv(Rcov_mat)*(R*beta);
wald_P = 1 - chi2cdf(wald,q);