PCMCI: Causal Discovery in Multivariate Crypto Time Series

PCMCI由Jakob Runge（2018、2020）提出，是一种用于多变量时间序列的因果发现算法。它将图模型学习的PC算法与瞬时条件独立（MCI）检验结合起来；其条件集合被设计为同时消除三类影响：共同驱动因素、间接路径，以及源变量的自相关。最后这一点尤其值得仔细阅读——它正是PCMCI不会继承朴素lead-lag检验在动量很强的序列上出现高假阳性率的原因。

本文从第一性原理拆解这两阶段构造，使用tigramite库实现，并列出该方法在成为任何交易输入之前仍需要完成的真实数据研究。

它为现有内容补充了什么

博客已经完成了相关性方面的讨论。加密货币的共同波动很大程度上是共同因子造成的——交易对之间的信号相关性显示，在十个交易对组成的book中，单独PC1就吸收了65%的方差。高pairwise correlation并不等于可交易关系，因此统计套利与配对交易转而使用协整。而且这种依赖并不稳定：DCC-GARCH动态相关性列出了静态样本相关性所隐含的假设：依赖关系平稳、边际波动率恒定、正反方向对称，而这些都不成立。

这些方法都不能产生有方向的图。它们告诉你BTC、ETH和SOL一起波动，且共同因子解释了大部分波动；但不能告诉你BTC是否直接驱动SOL，还是这种关联完全通过ETH传递。Granger因果性通常用于寻找方向，但标准Granger是二元方法，无法将未观测的共同驱动因素与真正的因果关系分开。

这里还有第三种、而且比表面看起来更重要的失效模式。每个加密货币收益序列都带有正自相关，而自相关会放大在两条此类序列之间进行的依赖检验统计量：两个相互独立但具有自相关的过程，在有限样本中也会看起来显著相关。相关性分析把它当作样本量修正（signal-correlation文章中的effective_N讨论）。但对因果发现而言，这不是修正项，而是核心设计问题——Stage 2的条件集合正是为了解决它。

PCMCI算法

PCMCI分两个阶段运行。第一阶段为每个变量识别候选因果父节点集合。第二阶段使用精心构造的条件独立检验测试每条候选连接。

Stage 1：PC-Stable条件选择（PC1）

Stage 1的目标是为每个变量$X_t^j$找到其真实因果父节点集合的超集$\widehat{\mathcal{B}}(X_t^j)$。该集合稍后会作为Stage 2的条件集合。

算法从所有滞后变量的完整集合开始，把它们都作为潜在父节点:

$\widehat{\mathcal{B}}^0(X_t^j) = \{ X_{t-\tau}^i : i = 1, \ldots, N, \quad \tau = 1, \ldots, \tau_{\max} \}$

随后，算法迭代删除与$X_t^j$条件独立的变量。在第$p$次迭代中，算法测试每个候选父节点$X_{t-\tau}^i$与$X_t^j$的独立性，条件是剩余最强的$p$个父节点（不包括$X_{t-\tau}^i$本身）:

$X_{t-\tau}^i \perp\!\!\!\perp X_t^j \mid \mathbf{S} \quad \Rightarrow \quad \text{remove } X_{t-\tau}^i \text{ from } \widehat{\mathcal{B}}(X_t^j)$

其中，$\mathbf{S}$是包含$p$个变量的集合，来自$\widehat{\mathcal{B}}(X_t^j) \setminus \{X_{t-\tau}^i\}$中与$X_t^j$关联最强的变量。

"Stable"的含义是：同一迭代$p$内的删除决策都基于上一轮迭代的父节点集合，从而避免结果依赖处理顺序。这对可复现性至关重要。

此阶段的显著性水平$\alpha_{\text{PC}}$会有意设置得宽松（通常为0.2或更高，也可以通过tigramite中的pc_alpha=None选项自动确定）。这里的目标不是控制假阳性，而是在尽可能剔除无关变量的同时保留所有真实父节点。$\widehat{\mathcal{B}}$中有少量假阳性可以接受，但不能出现假阴性。

Stage 2：瞬时条件独立（MCI）检验

得到所有变量的估计父节点集合$\widehat{\mathcal{B}}(X_t^j)$后，Stage 2使用MCI统计量测试每条潜在因果连接$X_{t-\tau}^i \to X_t^j$:

$X_{t-\tau}^i \perp\!\!\!\perp X_t^j \mid \widehat{\mathcal{B}}(X_t^j) \setminus \{X_{t-\tau}^i\}, \quad \widehat{\mathcal{B}}(X_{t-\tau}^i)$

关键创新在于条件集合。我们同时进行以下条件化:

1. 目标变量的父节点 $\widehat{\mathcal{B}}(X_t^j) \setminus \{X_{t-\tau}^i\}$：消除共同驱动因素以及进入$X_t^j$的其他间接路径的影响。
2. 源变量的父节点 $\widehat{\mathcal{B}}(X_{t-\tau}^i)$：消除源变量的自相关，避免检验统计量膨胀。

同时对两组集合进行条件化后，MCI实际上隔离了$X_{t-\tau}^i$对$X_t^j$在特定滞后$\tau$下的直接、瞬时因果影响，并扣除了所有混杂因素、中介和自相关影响。

检验统计量可以使用任何有效的条件独立检验。对于线性关系，偏相关是标准选择。对于非线性依赖，可以使用通过$k$近邻方法估计的条件互信息（CMI），或GPDC（Gaussian Process Distance Correlation）检验。

偏相关作为默认检验

对于大多数金融应用，偏相关检验（tigramite中的ParCorr）是主力工具。给定变量$X$、$Y$和条件集合$\mathbf{Z}$，偏相关为:

$\rho_{XY \cdot \mathbf{Z}} = \text{corr}(r_X, r_Y)$

其中$r_X$和$r_Y$分别是将$X$和$Y$对$\mathbf{Z}$回归后得到的残差:

$r_X = X - \hat{X}(\mathbf{Z}), \quad r_Y = Y - \hat{Y}(\mathbf{Z})$

换句话说，偏相关衡量的是去除$X$与$Y$之间由条件变量$\mathbf{Z}$产生的线性影响后，两者的线性关联。在高斯数据的条件独立原假设下，检验统计量服从已知分布，因此可以精确计算$p$值。

偏相关取值在$[-1, 1]$之间，可以自然地对因果连接强度排序。例如，若在两个父节点集合条件下，BTC滞后1期与ETH滞后0期的偏相关为0.3，含义是："去除所有已识别的混杂因素和自相关后，BTC收益的一次标准差冲击会预测下一期ETH收益0.3个标准差的变动。"

从检验结果到因果图（DAG）

PCMCI的输出是一对矩阵:

• val_matrix：形状为$(N, N, \tau_{\max}+1)$，包含每条有向连接$X_{t-\tau}^i \to X_t^j$的MCI检验统计量（例如偏相关）。
• p_matrix：形状相同，包含对应的$p$值。

要构造因果图，我们在经过修正的显著性水平上对$p$矩阵设阈值:

$\text{Link } X_{t-\tau}^i \to X_t^j \text{ exists if } p_{ij\tau} < \alpha_{\text{corrected}}$

这里的修正绝不是可以一笔带过的细节。PCMCI检验$N^2 \times \tau_{\max}$条按构造相互依赖的连接——每条进入$X_t^j$的连接都与其他连接共享条件集合，而源序列本身也彼此相关。正是在这种情形下，普通Benjamini-Hochberg会失效，Benjamini-Yekutieli的调和惩罚$c(M) = \sum_{j=1}^{M} 1/j$才可以在检验之间任意依赖时保持有效；deflated Sharpe与多重检验具体衡量了使这一点落地的原假阳性发现率。发布任何图时，都应同时报告检验的连接数量和修正后的阈值。

得到的结构是一个时间序列图（TSG），即一种有向图，其中:

• 节点表示特定滞后下的变量。
• 有向边表示带有对应滞后的因果连接。
• 自环表示自回归效应。

这个TSG还可以折叠为汇总图，只显示变量之间因果连接的存在与方向（跨滞后聚合），通常更便于解释。

使用Tigramite实现

Tigramite是PCMCI的参考实现，由德国航空航天中心（DLR）Jakob Runge团队开发和维护。它提供干净的数据处理API、多种条件独立检验、多个PCMCI变体以及内置可视化。

安装

pip install tigramite

Recovery sanity check（不是结果）

在市场数据上运行之前，先确认管道能够找到你植入的结构是值得的。这是实现的单元测试，仅此而已——下面的ground truth被手工写入数据生成循环，因此恢复它只能证明代码可用，完全不能说明加密货币。请按这个范围理解。

import numpy as np
import tigramite
from tigramite import data_processing as pp
from tigramite.pcmci import PCMCI
from tigramite.independence_tests.parcorr import ParCorr
from tigramite import plotting as tp


np.random.seed(42)
T, N = 2000, 5
var_names = ["BTC", "ETH", "SOL", "BNB", "AVAX"]
data = np.zeros((T, N))
noise = np.random.randn(T, N) * 0.5

for t in range(2, T):
    data[t, 0] = 0.5 * data[t-1, 0] + noise[t, 0]
    data[t, 1] = 0.5 * data[t-1, 1] + 0.4 * data[t-1, 0] + noise[t, 1]
    data[t, 2] = 0.5 * data[t-1, 2] + 0.3 * data[t-1, 1] + noise[t, 2]
    data[t, 3] = 0.5 * data[t-1, 3] + 0.25 * data[t-1, 0] + noise[t, 3]
    data[t, 4] = 0.5 * data[t-1, 4] + 0.2 * data[t-2, 3] + noise[t, 4]

dataframe = pp.DataFrame(
    data,
    datatime=np.arange(T),
    var_names=var_names,
)

parcorr = ParCorr(significance="analytic")

pcmci = PCMCI(
    dataframe=dataframe,
    cond_ind_test=parcorr,
    verbosity=1,
)

results = pcmci.run_pcmci(
    tau_max=4,        # test lags up to 4 hours
    tau_min=1,        # only lagged (not contemporaneous) links
    pc_alpha=None,    # auto-select alpha for condition selection
    alpha_level=0.01, # significance threshold for final MCI test
)

print("\n--- Significant causal links ---")
pcmci.print_significant_links(
    p_matrix=results["p_matrix"],
    val_matrix=results["val_matrix"],
    alpha_level=0.01,
)

tp.plot_graph(
    val_matrix=results["val_matrix"],
    p_matrix=results["p_matrix"],
    var_names=var_names,
    link_colorbar_label="MCI (partial corr.)",
    node_colorbar_label="Auto-MCI",
    alpha_level=0.01,
    figsize=(10, 6),
)

tp.plot_time_series_graph(
    val_matrix=results["val_matrix"],
    p_matrix=results["p_matrix"],
    var_names=var_names,
    link_colorbar_label="MCI (partial corr.)",
    alpha_level=0.01,
    figsize=(14, 6),
)

阅读输出

在植入的VAR上，print_significant_links返回:

Variable BTC has 0 causal parent(s):

Variable ETH has 1 causal parent(s):
    BTC (lag -1): val = 0.38, p = 0.000

Variable SOL has 1 causal parent(s):
    ETH (lag -1): val = 0.28, p = 0.000

Variable BNB has 1 causal parent(s):
    BTC (lag -1): val = 0.24, p = 0.000

Variable AVAX has 1 causal parent(s):
    BNB (lag -2): val = 0.19, p = 0.000

植入的结构被恢复了，其中包括一个重要的negative result：没有直接的BTC -> SOL边，尽管BTC和SOL在该样本中高度相关。MCI将这种关联完全归因于ETH这个中介变量。这正是双重条件集合应产生的行为，验证它是这段代码的全部意义。这是一个通过的unit test，而不是关于crypto的finding。

使用真实市场数据

PCMCI需要弱平稳输入，因此应提供standardized log-return而不是price——ADF machinery以及原始price为何会失败的原因见统计套利与配对交易。得到standardized return matrix后，tigramite-specific部分只需两行:

data = log_returns.values          # (T, N) standardized log-returns
var_names = list(log_returns.columns)

dataframe = pp.DataFrame(data, var_names=var_names)

以下两点tigramite细节不同于常规return pipeline:

1. 缺失数据：Tigramite支持masked array来处理缺失observation。使用dataframe.mask标记gap——不要forward-fill，因为伪造的observation会传播到包含它的每个条件集合中。
2. 频率对齐：所有序列必须位于同一时间网格。对crypto来说通常很直接，因为交易所提供同步的OHLCV data。

该方法需要的研究

以上都是machinery。在该machinery指向真实数据并报告发现（包括答案是"没有稳定关系"）之前，这篇文章不适合在博客发布。需要执行的run如下:

• Setup：对真实basket的小时return运行PCMCI+，明确日期范围；根据domain reasoning固定tau_max，对全部$N^2 \times \tau_{\max}$条连接使用alpha_level = 0.01和Benjamini-Yekutieli修正。按照deflated Sharpe标准，报告检验连接数和修正阈值。
• 需要报告的结果：哪些edge通过修正、它们的偏相关幅度，以及有多少在out-of-sample中存活。
• 至少一项稳定性测量：没有它，这篇文章就只是带crypto ticker的tigramite README。可以报告rolling window之间edge持续的比例；或者恢复图对tau_max和pc_alpha的敏感性；或者在同一data上ParCorr与CMIknn的分歧。

如果图不稳定，这是可以发布的结果，而不是失败。"五个主要币种上的PCMCI：因果图无法在rolling window中保持稳定"符合博客在诚实的negative result中坚持的路线，也比一个能运行的tutorial更有用。

选择正确的参数

最大滞后$\tau_{\max}$

该参数限制因果发现的时间范围。设置过低可能漏掉传播缓慢的影响；设置过高则会增加计算成本和多重检验负担。

对于小时级crypto data，$\tau_{\max} = 6$到$24$是合理范围。对于日线data，$\tau_{\max} = 5$到$10$可以捕捉大多数lead-lag关系。这里需要domain knowledge：如果你知道funding rate effect需要8小时传播，就应设置$\tau_{\max} \geq 8$。

条件选择显著性$\alpha_{\text{PC}}$

设置pc_alpha=None后，tigramite会使用Akaike Information Criterion自动选择该参数，这是推荐默认值。如果希望手动设置，0.1到0.4通常效果良好。更低的值会让Stage 1更激进（保留更少父节点），降低计算成本，但有删除真实父节点的风险。

最终显著性水平$\alpha$

这是标准的假设检验阈值。探索性分析使用$\alpha = 0.05$即可。任何会进入交易决策的内容都应使用$\alpha = 0.01$或更严格的值，并且要在此前描述的dependence-aware修正之后应用——nominal level不是你实际检验的level。

PCMCI+：加入同时发生的连接

标准PCMCI只发现滞后因果连接（$\tau \geq 1$）。但在信息以秒级跨资产传播的crypto市场中，小时采样意味着许多因果效应会表现为同时发生（$\tau = 0$）。

PCMCI+（Runge，2020）扩展PCMCI，同时发现滞后和同时发生的因果连接。默认情况下，同时发生的连接是无向的（因为同一时间步内的时间顺序无法区分原因与结果），但其中一些可以依据标准PC算法的定向规则（碰撞点检测、无环约束）进行定向。

results_plus = pcmci.run_pcmciplus(
    tau_max=4,
    tau_min=0,         # include contemporaneous links
    pc_alpha=None,
)

对于小时或更低频率的crypto应用，PCMCI+通常比标准PCMCI更合适，因为许多跨资产效应发生得比采样间隔更快。

非线性扩展

金融时间序列经常呈现非线性依赖（例如波动率聚集、依赖状态的lead-lag effect）。Tigramite提供多种非线性条件独立检验:

Gaussian Process Distance Correlation（GPDC）

from tigramite.independence_tests.gpdc import GPDC

gpdc = GPDC(significance="analytic", gp_params=None)
pcmci = PCMCI(dataframe=dataframe, cond_ind_test=gpdc)

GPDC使用Gaussian process regression去除条件集合的影响，然后对残差应用distance correlation。它检测非线性效应的能力强于ParCorr，但速度显著更慢。

条件互信息（CMIknn）

from tigramite.independence_tests.cmiknn import CMIknn

cmiknn = CMIknn(significance="shuffle_test", knn=0.1, shuffle_neighbors=5)
pcmci = PCMCI(dataframe=dataframe, cond_ind_test=cmiknn)

CMIknn使用$k$近邻方法估计条件互信息。它完全是nonparametric的，可以检测任意函数依赖。代价是计算成本较高，并且需要更多data才能获得统计功效。

对于大多数crypto trading application，应从ParCorr开始。只有在有具体证据表明存在非线性因果机制且data足够时（通常$T > 5000$），才切换到非线性检验。

因果图可以接入哪里

因果图能提供而相关矩阵无法提供的一点是：出度。在20-50个asset上运行PCMCI，拥有许多出边的节点就是information leader——它们的波动对book其余部分具有预测内容，而这是对称相关矩阵在结构上无法给出的方向性结论。因子分解告诉你BTC主导variance；out-degree ranking告诉你BTC主导是因为它领先，还是因为包括BTC在内的所有asset都加载在同一个macro shock上。

如果最终能恢复稳定图，它接入博客已有工作的地方也很清楚。图拓扑在rolling window中的变化是structural-break signal，主题见使用HMM进行regime detection——open question是因果图能否比state model更早检测transition。将PCMCI应用于一个asset跨venue比较，可以研究cross-venue price leadership，这已经在smart order routing中以毫秒分辨率测量；小时bar上的图必须胜过那个结果，而不是重复它。funding mechanism每8小时产生的perp-spot link见funding rate arbitrage，这是一个已知答案的test case：PCMCI应该恢复一个可以独立验证的因果结构，因此它是validation target而不是discovery。

这些都不是finding，而是等待上述study验证的hypothesis。

局限与注意事项

PCMCI是强大的工具，但它有实践者必须理解的重要局限:

1. 因果充分性假设：PCMCI假设所有相关变量都被观测。如果隐藏的共同驱动因素（例如whale的交易活动、尚未发布的新闻）同时影响两个被观测asset，PCMCI可能错误地报告它们之间存在直接因果连接。LPCMCI变体允许潜在混杂因素，可部分解决这一问题，但代价是返回更少的有向edge。

2. 平稳性假设：假设因果结构在分析window内保持不变。实践中crypto market dynamics变化很快。使用rolling-window analysis检测structural break。

3. 线性与非线性：使用ParCorr只能检测线性因果效应。非线性机制（例如"只有当BTC跌幅超过5%时，BTC才会导致ETH下跌"）在线性检验中不可见。

4. 采样频率很重要：快于采样频率发生的因果效应，在标准PCMCI中会表现为同时发生（$\tau = 0$）连接，其方向可能不明确。使用PCMCI+并考虑更高频data。

5. 多重检验：$N^2 \times \tau_{\max}$个相互依赖的检验按上述方式修正——关于依赖结构为何排除显而易见的选择，参见deflated Sharpe与multiple testing。

6. 样本量要求：可靠的因果发现需要足够data。粗略指南是：使用ParCorr时，$T > 500$适用于$N < 10$个变量；非线性检验或更大的变量集合则目标为$T > 2000$。

PCMCI与其他方法的比较

PCMCI相对于Granger因果性和transfer entropy的主要优势是MCI检验：它在多变量环境中正确处理了自相关和共同驱动因素。这正是crypto market遇到的场景：几十个彼此相关且具有自相关的asset同时互动。

结论

值得带走的贡献是两阶段构造：先用稀疏条件选择限制父节点集合，再使用一个条件集合同时包含源和目标的父节点的MCI检验。后半部分正是PCMCI区别于Granger和transfer entropy的地方，也正是该方法能够用于自相关程度如crypto return一样高的序列的原因。

该方法没有自带的是证据。恢复出的图是关于information flow的hypothesis，而不是signal。在本博客中，除非它在真实data上运行，并明确basket、date range、dependence-aware correction以及至少一项测量图在window移动时是否仍然成立，否则不算任何结果。在那次run完成前，请把这里的一切视为tooling。

如果图最终不稳定，那就是文章本身。这不会是本博客第一个比可运行tutorial更有价值的negative result。

参考文献

• Runge, J., Nowack, P., Kretschmer, M., Flaxman, S., and Sejdinovic, D. (2019). Detecting and quantifying causal associations in large nonlinear time series datasets. Science Advances, 5(11), eaau4996.
• Runge, J. (2020). Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), PMLR 124:1388-1397.
• Tigramite documentation: https://jakobrunge.github.io/tigramite/
• Tigramite GitHub repository: https://github.com/jakobrunge/tigramite