CEDAR：自回归过程的因果边发现

# CEDAR：自回归过程的因果边发现 来源：https://arxiv.org/html/2607.20696 ###### 摘要 我们提出 CEDAR（自回归过程的因果边发现），一种基于约束的方法，用于稀疏自回归时间序列中的滞后因果边发现。CEDAR 使用 AR(1) 残差化、\(\mathcal{U}\)-中心距离相关性筛选候选跨变量滞后，然后对每个显著的跨变量滞后候选执行两个有针对性的条件独立性检验，并每个有序对最多接受一个滞后。稳定的 MCI 剪枝步骤移除间接边，可选的确定性 C 节点则针对指定的趋势型非平稳性进行调整。在只有少数滞后通过筛选的稀疏场景中，CEDAR 在筛选后需要 \(O(d^2)\) 次 CI 检验，同时保留边级别的可解释性。当数据稀缺且变量表现出滞后-1 自动态时，CEDAR 最为有效；随着 \(T\) 增长或当高阶自回归或同时多滞后效应普遍存在时，具有更丰富条件集的方法更为可取。¹¹脚注文本：源代码和实验脚本可在 Causal-TS 包中获得；所有结果均可完全复现。

## 1 引言

从多元时间序列中进行因果发现对于气候科学、金融、神经科学和工业过程监测等领域至关重要 (Peters 等，2017 (https://arxiv.org/html/2607.20696#bib.bib1)；Runge 等，2019 (https://arxiv.org/html/2607.20696#bib.bib2)；Spirtes 等，2000 (https://arxiv.org/html/2607.20696#bib.bib3))。时间序列因果发现比独立同分布场景更具挑战性：观测值存在时间依赖性，变量可能跨多个滞后相互影响，强自相关或共同驱动者可能造成虚假关联。尽管时间顺序提供了有用信息——原因必须先于结果——但它也引入了更大且结构化的滞后变量搜索空间。这些挑战在现实系统中进一步加剧，其中因果机制可能是非平稳的，变量可能依赖于高阶自回归结构。因此，从时间序列中进行可靠的因果发现需要能够区分直接因果效应与自相关、滞后依赖性、间接链接以及随时间变化的机制的方法。

#### 相关工作。 时间序列因果发现算法大致分为两种范式。
*基于评分的方法*——DYNOTEARS (Pamfil 等，2020 (https://arxiv.org/html/2607.20696#bib.bib4))、VARLiNGAM (Hyvärinen 等，2010 (https://arxiv.org/html/2607.20696#bib.bib5))、神经格兰杰因果 (Tank 等，2021 (https://arxiv.org/html/2607.20696#bib.bib6))、CUTS+ (Cheng 等，2024 (https://arxiv.org/html/2607.20696#bib.bib7))——将发现问题表述为优化问题，搜索最能根据评分标准解释数据的图。它们能很好地扩展维度，并可以自然地纳入稀疏性惩罚，但需要针对二值图输出进行事后阈值调整，并且由于边缘源自全局目标而非每条链接的显式证据，单个边缘更难解释。

*基于约束的方法*——PCMCI/PCMCI+ (Runge 等，2019 (https://arxiv.org/html/2607.20696#bib.bib2)；Runge，2020 (https://arxiv.org/html/2607.20696#bib.bib8))、LPCMCI (Gerhardus and Runge，2020 (https://arxiv.org/html/2607.20696#bib.bib9))、CDNOTS (Sadeghi 等，2025 (https://arxiv.org/html/2607.20696#bib.bib10))——提供了一个更直接且可解释的统计过程，基于条件独立性 (CI) 检验添加或移除边。它们的主要代价是统计和计算上的：在稠密或高维图中，条件集可能迅速增大，当样本有限时，检验功效会降低。

#### 我们的方法。 为了解决这种低样本场景，我们引入 CEDAR，它从时间结构而非在大父节点子集上搜索来构建有针对性的条件集。该方法保留了边级别的可解释性：每条报告的链接都得到显式 CI 检验的支持，然后通过稳定的 MCI 剪枝步骤重新检查。CEDAR 建立在 SyPI (Mastakouri 等，2021 (https://arxiv.org/html/2607.20696#bib.bib11)) 的基础上，后者是一种针对具有指定汇目标的自回归时间序列的因果特征选择方法。我们通过自动化滞后选择、利用时间无环性放宽汇目标限制，并剪枝间接链接，将所有有序对推广。由此产生的过程是一个适用于稀疏类 AR(1) 系统的实用有限样本方法；其 oracle 保证需要附加的结构条件，这些条件在附录 A (https://arxiv.org/html/2607.20696#A1) 中给出。

#### 贡献。 我们的主要贡献是：
1. 1. 具有解析筛选能力的非线性滞后选择：原始 SyPI 使用 Lasso 回归进行滞后选择，这仅限于线性依赖。我们将其替换为 \(\mathcal{U}\)-中心距离相关性 (Székely and Rizzo，2014 (https://arxiv.org/html/2607.20696#bib.bib12))——一种非参数度量，能够捕捉任意非线性关联，且在独立性下期望为零——并结合解析的 t 检验 (Székely and Rizzo，2013 (https://arxiv.org/html/2607.20696#bib.bib13)) 进行快速筛选，使 CEDAR 能够在线性和非线性系统中筛选候选滞后。
2. 2. 稳定的 MCI 剪枝：在第一阶段的二条件 CI 测试之后，某些检测到的边可能是间接的——通过其他变量中介。我们引入一个稳定的、与顺序无关的瞬时条件独立性 (MCI) 剪枝步骤，该步骤以源和目标两者的已发现父节点为条件，在单次遍历中移除所有此类边，而不受迭代方法中的顺序依赖影响。
3. 3. 通过 C 节点处理非平稳性：现实世界的时间序列往往随时间发生分布变化，从而在共享共同趋势的变量之间产生虚假关联。我们引入一个编码多项式时间趋势的确定性 C 节点，将其视为纯原因，在产生虚假边之前针对指定的确定性趋势型非平稳性进行调整——这是原始 SyPI 框架所不具备的能力。

## 2 背景与符号

令 \(\mathbf{X} = (X^1, \ldots, X^d)\) 为一个 \(d\) 变量时间序列，有 \(T\) 次观测值和最大滞后 \(L\)。我们关注跨变量滞后图：\(G_X[c, e, \ell] = 1\)，对于 \(c \neq e\) 且 \(\ell \geq 1\)，表示直接边 \(X^c(t-\ell) \to X^e(t)\)。滞后-1 自环由 AR(1) 模型假定，且不被学习；可选的 C 节点边在滞后 0 处单独存储为 \(G_C[r, e]\)。

CEDAR 改编了 SyPI (Mastakouri 等，2021 (https://arxiv.org/html/2607.20696#bib.bib11)) 的二条件思想，该思想测试候选 \(X^c(t-w) \to X^e(t)\) 是否是指定目标的直接原因。对于每个候选源 \(c\)，SyPI 以来自每个其他候选时间序列的一个移位节点 \(\mathcal{S}^c = \{X^k(t-w_k-1): k \neq c\}\) 为条件，并测试
(C1) \(X^c(t-w) \not\perp\perp X^e(t) \mid \{\mathcal{S}^c, X^e(t-1)\}\)，
(C2) \(X^c(t-w-1) \perp\perp X^e(t) \mid \{\mathcal{S}^c, X^c(t-w), X^e(t-1)\}\)。
条件 1 建立剩余依赖性，而条件 2 验证所选滞后阻断了源更深层过去的影响。该逻辑依赖于严格的 AR(1) 自动态以及不存在滞后不匹配的旁路路径；附录 A (https://arxiv.org/html/2607.20696#A1) 陈述了确切的充分条件及一个反例。

对于滞后筛选，CEDAR 使用 Székely and Rizzo (2014 (https://arxiv.org/html/2607.20696#bib.bib12)) 的有偏校正平方距离相关 \( \mathcal{R}^{*2}_n \)。由于其底层的 \(\mathcal{U}\)-中心距离协方差在独立性下期望为零，它提供了一种快速的非线性边际评分；定义见附录 B (https://arxiv.org/html/2607.20696#A2)。

## 3 方法

CEDAR 分三个主要阶段进行：(1) 滞后选择，(2) 通过条件 1 和 2 进行 CI 测试，以及 (3) MCI 剪枝。算法 1 (https://arxiv.org/html/2607.20696#alg1) 给出了完整过程；我们下面详细描述每个组件。

算法 1 CEDAR 总结。完整实现细节在第 3.1 (https://arxiv.org/html/2607.20696#S3.SS1)–3.6 (https://arxiv.org/html/2607.20696#S3.SS6) 节中阐述。

**要求**：数据 \(\mathbf{X} \in \mathbb{R}^{T \times d}\)，最大滞后 \(L\)，CI 检验 \(\mathcal{T}\)，水平 \(\alpha_1, \alpha_2, \alpha_{\rm lag}\)；可选的 C 节点基函数 \(g_1, \ldots, g_R\)。
1: 可选地附加确定性纯原因列 \(C_r(t) = g_r(t/T)\)。
2: 对每个非 C 目标进行残差化：\(\tilde{X}^e(t) = X^e(t) - \hat{\beta}_e X^e(t-1)\)。
3: 使用 \(\mathcal{R}_n^{*2}(X^c(t-\ell), \tilde{X}^e(t))\) 对所有跨变量滞后进行评分；对 C 节点仅在滞后 0 处评分。
4: 在 \(\alpha_{\rm lag}\) 水平下应用 BH 筛选。
5: **对于** 每个目标 \(e\) 和具有筛选后滞后(们)的候选源 \(c\)：
6:   从其他筛选后的非 C 候选者在 \(t-\hat{w}_{k,e}-1\) 处以及筛选后的 C 节点在滞后 0 处构建 \(\mathcal{S}^c\)。
7:   **如果** \(c\) 是 C 节点：
8:     仅测试条件 1，如果显著则设置 \(G_C[r, e] = 1\)。
9:   **否则**：
10:    按评分降序测试筛选后的滞后；接受第一个通过条件 1–2 的滞后并设置 \(G_X[c, e, w] = 1\)。
11:  **结束如果**
12: **结束循环**
13: 在相同的初始化 \(G_X\) 上运行所有 MCI 剪枝测试，并同时移除接受的独立性。

### 3.1 滞后选择

对于每对变量 \((X^c, X^e)\) 和候选滞后 \(\ell \in \{1, \ldots, L\}\)，我们使用对 AR(1) 残差化目标进行的 \(\mathcal{U}\)-中心距离相关性 (Székely and Rizzo，2014 (https://arxiv.org/html/2607.20696#bib.bib12)) 对滞后重要性进行评分。具体来说，我们首先移除目标的线性自回归分量：
\[\tilde{X}^e(t) = X^e(t) - \hat{\alpha}_e - \hat{\beta}_e X^e(t-1), \quad \hat{\beta}_e = \frac{\sum_t (X^e(t) - \bar{X}^e)(X^e(t-1) - \bar{X}^e_-)}{\sum_t (X^e(t-1) - \bar{X}^e_-)^2}, \quad \hat{\alpha}_e = \bar{X}^e - \hat{\beta}_e \bar{X}^e_-,\]
其中 \(\bar{X}^e\) 和 \(\bar{X}^e_-\) 分别表示 \(X^e(t)\) 和 \(X^e(t-1)\) 的样本均值（即带截距的普通最小二乘）。然后计算滞后重要性评分：
\[M[c, e, \ell] = \mathcal{R}_n^{*2}(X^c(t-\ell), \tilde{X}^e(t)),\]
其中 \(n = T - L\) 是有效样本量。残差化防止目标的 AR 动态在无关滞后上夸大重要性评分。注意，这是简单的线性残差化，并非 Székely–Rizzo 意义上的偏距离相关性 (Székely and Rizzo，2014 (https://arxiv.org/html/2607.20696#bib.bib12))；如果目标的自回归动态是非线性的，则移除不完整。

#### 解析 t 检验。 由于在独立性下 \(\mathbb{E}[\operatorname{dCov}_n^{*}] = 0\) ((Székely and Rizzo，2014 (https://arxiv.org/html/2607.20696#bib.bib12)) 中的命题 1)，我们可以利用 Székely and Rizzo (2013 (https://arxiv.org/html/2607.20696#bib.bib13)) 推导出的 \(\mathcal{R}_n^{*2}\) 的渐近零分布。定义 \(v = n(n-3)/2\)；在 \(H_0: X \perp\perp Y\) 下，统计量
\[T_n = \frac{\sqrt{v-1} \cdot \mathcal{R}_n^{*2}}{\sqrt{1 - (\mathcal{R}_n^{*2})^2}}\]
近似服从自由度为 \(v-1\) 的 Student-t 分布，产生 p 值 \(p = 1 - F_{t(v-1)}(T_n)\)。这消除了对置换检验的需要。t 检验假设样本独立同分布；对于强自相关数据，有效样本量小于 \(n\)，可能导致 p 值过于宽松。本文报告的所有实验均使用解析检验。附录 C (https://arxiv.org/html/2607.20696#A3) 报告了 AR(1) 零假设下的校准：I 类错误略有膨胀（名义 5% 时为 6.4–9.0%），在强自相关和小 \(T\) 时最严重。

### 3.2 多滞后筛选与 BH 校正

现实因果系统通常表现出多个滞后的效应。原始 SyPI 框架测试每对的最小滞后 (Mastakouri 等，2021 (https://arxiv.org/html/2607.20696#bib.bib11))，识别单一因果滞后。CEDAR 将所有 \(d(d-1)L\) 个跨变量滞后显著性假设应用 Benjamini–Hochberg (BH) 校正 (Benjamini and Hochberg，1995 (https://arxiv.org/html/2607.20696#bib.bib14)) 作为多重检验筛选启发式。设 \(p_{(1)} \leq p_{(2)} \leq \ldots \leq p_{(m)}\) 为来自滞后 t 检验的有序 p 值，其中 \(m = d(d-1)L\) 个总跨变量假设。如果满足下式，BH 过程拒绝假设 \(i\)：
\[p_{(i)} \leq \frac{i}{m} \cdot \alpha_{\text{lag}},\]
其中 \(\alpha_{\text{lag}}\) 是筛选阈值（默认 0.05）。如果使用 C 节点，其滞后 0 的筛选 p 值也包含在同一筛选族中，除非另有说明。BH 校正后，每对 \((X^c, X^e)\) 可能有多个显著滞后。我们按重要性降序（\(\mathcal{R}_n^{*2}\) 最高者优先）测试这些滞后，并仅接受第一个同时通过条件 1 和条件 2 的滞后（keep = “first”）。也就是说，CEDAR 每个有序对选择一个滞后；多滞后筛选识别最佳候选滞后，而非恢复一对之间的所有同时滞后。

### 3.3 稳定的 MCI 剪枝

在初始基于 CI 的发现（阶段 1–2）之后，某些检测到的边可能是间接的——通过已发现图中的其他变量中介。CEDAR 应用一个受 PCMCI (Runge 等，2019 (https://arxiv.org/html/2607.20696#bib.bib2)) 启发但适用于成对发现设置的瞬时条件独立性 (MCI) 剪枝步骤。对于每条已发现的边 \(X^c(t-\tau) \to X^e(t)\)，我们测试：
\[X^c(t-\tau) \perp\perp X^e(t) \mid \mathcal{S}_{\text{MCI}},\]
其中 MCI 条件集结合了目标和源两者的父节点，各自在其对应时刻：
\[\mathcal{S}_{\text{MCI}} = \bigl(\text{Pa}(X^e(t)) \setminus \{X^c(t-\tau)\}\bigr) \cup \{X^e(t-1)\} \cup \text{Pa}(X^c(t-\tau)) \cup \{X^c(t-\tau-1)\} \cup \mathcal{C}_e.\]
这里 \(\text{Pa}(X^e(t))\) 表示已发现的 \(X^e\) 的非 C 父节点，相对于 \(t\) 的滞后时间（不包括正在测试的边），\(X^e(t-1)\) 是目标的自回归项，\(\text{Pa}(X^c(t-\tau))\) 表示已发现的 \(X^c\) 的父节点，移位到