格兰杰因果关系检验代码

以下是 Python 中使用格兰杰因果关系检验的代码示例：

import numpy as np
from scipy.stats import pearsonr

def granger_causality(x, y, maxlag):
    """
    格兰杰因果关系检验
    :param x: 一维时间序列 x
    :param y: 一维时间序列 y
    :param maxlag: 最大滞后阶数
    :return: p 值
    """
    n = len(x)
    lag = range(1, maxlag+1)
    xlags = np.zeros((n - maxlag, maxlag))
    ylags = np.zeros((n - maxlag, maxlag))
    for i in range(maxlag):
        xlags[:, i] = x[maxlag-i-1:n-i-1]
        ylags[:, i] = y[maxlag-i-1:n-i-1]
    xy = np.hstack([ylags, xlags])
    yx = np.hstack([xlags, ylags])
    beta1 = np.linalg.inv(xy.T @ xy) @ xy.T @ y[maxlag:]
    beta2 = np.linalg.inv(yx.T @ yx) @ yx.T @ x[maxlag:]
    resid1 = y[maxlag:] - xy @ beta1
    resid2 = x[maxlag:] - yx @ beta2
    rss1 = resid1.T @ resid1
    rss2 = resid2.T @ resid2
    f = ((rss2-rss1)/maxlag) / (rss1/(n-2*maxlag-1))
    p = 1 - stats.f.cdf(f, maxlag, n-2*maxlag-1)
    return p

其中，`x` 和 `y` 分别表示两个时间序列，`maxlag` 表示最大滞后阶数，`p` 值表示格兰杰因果关系的显著性水平。如果 `p` 值小于某一显著性水平（例如 0.05），则可以拒绝零假设，即认为存在格兰杰因果关系。