#!/usr/bin/env python
# coding: utf-8

# # 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


# ## matlab实现
# ```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);
# ```