Python中多维时间序列与单变量时间序列的因果推断方法问询

1. 扩展Granger因果检验（基于VAR模型） #

Granger因果的核心逻辑是“利用原因变量的历史数据能提升结果变量的预测精度”，statsmodels自带的grangercausalitytests仅支持一维对一维检验，但可以通过向量自回归（VAR）模型扩展到多维输入场景：

实现步骤： #

• 将多维输入信号（如多通道振动传感器数据）与单变量信号合并为同一时间序列数据集
• 拟合VAR模型后，调用模型的test_causality方法，直接检验多维输入对单变量的Granger因果性

代码示例： #

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

# 构造模拟数据：2维输入X（X1,X2） + 单变量Y
np.random.seed(42)
time_steps = 100
X1 = np.cumsum(np.random.randn(time_steps))
X2 = np.cumsum(np.random.randn(time_steps))
Y = 0.3 * X1[1:] + 0.2 * X2[1:] + np.random.randn(time_steps-1)  # Y依赖X的滞后项

# 合并数据并删除缺失值
df = pd.DataFrame({'X1': X1[:-1], 'X2': X2[:-1], 'Y': Y})

# 拟合VAR模型，选择合适的滞后阶数
model = VAR(df)
results = model.fit(maxlags=2)

# 检验X1、X2共同对Y的Granger因果性
causality_result = results.test_causality(caused='Y', causing=['X1', 'X2'], kind='f')
print(causality_result.summary())

• 若输出的p值小于预设显著性水平（如0.05），则拒绝“X不Granger引起Y”的原假设，说明二者存在因果关联。

2. 基于因果森林的非线性因果检验 #

如果数据存在非线性关系，线性Granger检验可能失效，此时可采用因果森林挖掘多维输入对单变量的非线性因果效应，适合非平稳或非线性时间序列：

实现步骤： #

• 构造特征集：包含多维输入的滞后项、单变量自身的滞后项
• 借助因果推断库构建模型，通过比较“包含/不包含X特征”的预测误差，验证X对Y的因果性

代码示例： #

import numpy as np
import pandas as pd
from causalml.inference.tree import CausalTreeRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# 构造滞后特征：取1阶、2阶滞后项
df['X1_lag1'] = df['X1'].shift(1)
df['X1_lag2'] = df['X1'].shift(2)
df['X2_lag1'] = df['X2'].shift(1)
df['X2_lag2'] = df['X2'].shift(2)
df['Y_lag1'] = df['Y'].shift(1)
df['Y_lag2'] = df['Y'].shift(2)
df = df.dropna()

# 拆分特征与目标：X的滞后+Y的滞后作为特征，当前Y作为目标
X_full = df[['X1_lag1', 'X1_lag2', 'X2_lag1', 'X2_lag2', 'Y_lag1', 'Y_lag2']]
X_onlyY = df[['Y_lag1', 'Y_lag2']]
y = df['Y']

X_train_full, X_test_full, y_train, y_test = train_test_split(X_full, y, test_size=0.2, random_state=42)
X_train_onlyY, X_test_onlyY = train_test_split(X_onlyY, test_size=0.2, random_state=42)

# 训练仅用Y滞后的线性模型
lr = LinearRegression()
lr.fit(X_train_onlyY, y_train)
mse_onlyY = np.mean((lr.predict(X_test_onlyY) - y_test)**2)

# 训练包含X特征的因果树
ct = CausalTreeRegressor(max_depth=3, random_state=42)
ct.fit(X_train_full, y_train, treatment=np.ones(X_train_full.shape[0]))
mse_withX = np.mean((ct.predict(X_test_full) - y_test)**2)

print(f"仅用Y滞后的预测MSE: {mse_onlyY:.4f}")
print(f"加入X特征的预测MSE: {mse_withX:.4f}")

• 若加入X特征后预测误差显著降低，且因果树估计的平均处理效应（ATE）显著不为0，说明X对Y存在非线性因果关联。

3. 传递熵（Transfer Entropy）：基于信息论的因果检验 #

传递熵是信息论视角下的因果指标，衡量“X的历史信息能减少Y的不确定性”的程度，天然支持多维输入，适合非线性、非高斯分布的数据：

实现步骤： #

• 使用信息论库计算从多维X到单变量Y的传递熵
• 通过置换检验（打乱X的时间顺序）判断结果的统计显著性

代码示例： #

import numpy as np
from pyinform import transfer_entropy

# 提取多维X和单变量Y的时间序列（需与之前的数据对齐）
X_multi = df[['X1', 'X2']].values
Y_single = df['Y'].values

# 计算X到Y的传递熵，滞后阶数设为1
te_value = transfer_entropy(X_multi, Y_single, k=1)
print(f"X到Y的传递熵值: {te_value:.4f}")

# 置换检验验证显著性
np.random.seed(42)
perm_te = []
for _ in range(1000):
    perm_X = np.random.permutation(X_multi)
    perm_te.append(transfer_entropy(perm_X, Y_single, k=1))

# 计算p值：原始传递熵大于置换样本的比例
p_value = np.mean(np.array(perm_te) > te_value)
print(f"置换检验p值: {p_value:.4f}")

• 若p值小于显著性水平，说明X对Y的因果关联具有统计显著性。

内容的提问来源于stack exchange，提问作者Α Πι