跳转至

Statistical Inference for Local Granger Causality

作者: Yan Liu, Masanobu Taniguchi, Hernando Ombao
来源: Statistica Sinica
主题: 因果推断
相关性: 8/10
链接: 期刊页 · arXiv


一、领域脉络与小综述

这个方向是什么

本子方向解决的核心问题是:如何从(可能非平稳的)多元时间序列数据中,定义并检验随时间变化的(Granger)因果关系——即在频域上刻画并推断一个时间序列的分量是否在某个时刻超前于另一个分量提供预测信息。当前对这一问题的统计处理主要依赖于将因果关系嵌入到时变谱密度矩阵中,通过对谱函数的频域分解来定义“局部”的意义。成熟度:局部平稳过程(locally stationary processes)的理论基础已经较完备(Dahlhaus 1997),但把 Granger 因果从平稳推广到这一框架并完成严格的渐近推断(估计量分布 + 检验统计量分布)的工作此前不完整。本论文是填补这个缺口的一个系统工作。

发展脉络

从引言引用的工作出发,可以梳理出以下时间轴(每个节点均标注作者-年份和一句话判断,引用来源为用户提供的论文原文中的陈述):

  1. 奠基工作 I:频域 Granger 因果的经典形式 (Geweke 1982, 1984; Hosoya 1991);Brillinger (2001) 的专著《Time Series: Data Analysis and Theory》。

    • 作者原话 (引用句):“Granger causality has been employed to investigate causality relations between components of stationary multiple time series.”
    • 贡献:在线性、平稳、频域框架下定义了 Granger 因果的谱表示,给出了线性反馈测度(measure of linear feedback)的分解。
    • 留下口子:只适用于全局平稳假设,不能处理实时演化。
  2. 奠基工作 II:局部平稳过程的系统理论 (Dahlhaus 1997, 2000, 2001)。

    • 作者原话 (引用句):“Our method is based on the locally stationary process framework developed by Dahlhaus (1997).”
    • 贡献:定义了局部平稳过程——利用一个平滑的“时间”索引将非平稳过程局部地近似于平稳过程,并推导了局部周期图的渐近性质。
    • 留下口子:这个框架虽能处理谱的平滑演化,但并没有将因果结构作为目标来专门定义和检验。
  3. 主要进展:时变谱估计 (Guo et al. 2003; Ombao et al. 2005)。

    • 作者原话 (引用):“Guo, Dai and Ombao (2003) … proposed a time-varying autoregressive modeling approach.” “Ombao, von Sachs and Guo (2005) … extended the SLEX (smooth localized complex exponentials) approach to nonstationary processes.”
    • 贡献:将时间序列谱估计推广到时变情形,采用分段平稳 + 自适应窗口,使谱估计能够捕捉断点或平滑转变。
    • 留下口子:这些方法主要是为估计谱(而非因果检验) 而设计的,其因果推断需额外假设。
  4. 当前 frontier → 本文的位置

    • 作者原话 (引用句中隐含的判断):“Our proposed local Granger causality approach captures time-evolving causality relationships in nonstationary processes.”
    • 本文相对于之前工作的增量是:在 Dahlhaus 的局部平稳框架下,将 Granger 因果的频域检验显式地参数化并通过局部 Whittle 似然进行估计,并给出了估计量渐近正态检验统计量二次型分布的严格证明。
    • 不同于SLEX等分段方法,本文的局部窗口是通过核函数平滑移动的(类似局部Whittle),因此因果的时变是平滑的

子线索聚类

这些被引文献大致落在两条线索上:

  • 线索 A:时间序列因果推断的谱表示方法与频率戳
    代表:Geweke (1982, 1984), Hosoya (1991), Brillinger (2001), Dhamala et al. (2008) (后被引于本文讨论部分)。
    核心:在频域上定义因果性时,因果方向由“谱密度矩阵的零模式”表征——频率ω下从函数j→k无Granger因果等价于该频率处的谱密度矩阵分块对角块中相应的交叉谱密度处处为0。这条线索的不足是普遍假设过程全局平稳

  • 线索 B:时变谱估计与局部平稳过程
    代表:Dahlhaus (1997, 2000, 2001), Guo et al. (2003), Ombao et al. (2005)。
    核心:建立了时变谱估计的统计基础,大多采用局部窗 → 加权似然的结构。
    不足:B中鲜有研究直接聚焦于Granger因果的定义与检验;因果推断往往需要事后在估计谱上施加检验(如permutation test/Bootstrap),缺乏大样本有效性理论。

注:作者在引言中未直接与非谱方法(如时域 VAR 参数向量自回归 + 结构变化检验)进行比较——这构成一个可能的竞争路线,后文中会提及。

这个方向在追问的核心问题

  1. 如何定义“局部”Granger 因果以保证其在频域下的可识别,且与传统的全局平稳定义相容?
  2. 在时变谱的估计中,如何克服局部窗口内数据减少与参数自由度膨胀之间的权衡,同时仍给出可论证的渐近分布(而非仅仅依赖Bootstrap)?
  3. 能否在保持“局部”形式下,将 Granger 因果检验从线性扩展到非线性
  4. 对于高维(p >> n)局部平稳时间序列,频域因果推断是否可行?—— 当前工作仅处理有限维(p固定,n→∞的情形)。

⚠️ 作者的 framing(必须明确标注成“这是作者的说法”)

  • 作者把缺口 frame 成:“现有 Granger 因果方法假设全局平稳,而现实数据(脑电、金融)常是非平稳的,因此需要发展局部版本。” 作者在引言中称:“The proposed local Granger causality is a natural extension of the classical Granger causality in the frequency domain.”
  • 哪些竞争路线被回避? 作者几乎没有讨论时域含潜在变量情形下的 Granger 因果识别(即 instrument-free 方法),也未涉及非线性因果模型(如信息理论测度)。作者也回避了分段突变(非平滑) 的应用场景——论文的局部平稳框架本质上假设谱是时变但光滑的(假设 2.1 的高阶导数有界)。
  • 什么明显该被引 / 该存在、却没出现在 intro 里?
  • 高维时间序列因果推断(比如 Lustig et al. 2017 或 Basu & Michailidis 2015 关于 Group Lasso Granger 因果的渐近理论)——因为本文的公式自然可以扩展到高维,而作者并未讨这一线文献。
  • 时域分段线性 VAR 的结构断点检验(如 Qu & Perron 2007)——这是一种在时域直接检验时变因果关系的不同路线。
  • 这些缺失值得去查[属于作者的 omission,研究者可自问其是否有意回避]。

张力

未见明显对立引用。被引文献彼此和谐:Geweke (1982) 给的频域因果与 Dahlhaus (1997) 的局部平稳过程之间没有内在矛盾,只是领域不同。但一个张力值得注意:局部窗口(Local Whittle)与分段平稳(SLEX)哪个更适于处理因果结构的问题——作者在本文倾向平滑,但并未在理论上证明平滑窗口优于分段窗口。


二、最核心、最简单的例子 / 数学问题

第一步:把符号、模型、可观测数据交代清楚

符号:
- \( X_t = (X_{1t}, X_{2t}, \dots, X_{pt})^\top \):p-维复随机向量(时间序列),下标 t = 1, …, T。
- \( T \):总样本长度。
- \( u = t/T \in (0,1) \):规范化时间,作为连续参数出现在局部平稳过程的均值和自协方差中。
- \( \boldsymbol{\theta}(u) \):在时刻 u 处的有限维参数向量(例如时变 AR 系数,ARCH 系数等)。它是本论文直接估计的目标。
- \( f(\omega, u) \):时间 u 处、频率 ω 处的时变谱密度矩阵,p×p——从局部平稳过程的时变自协方差函数通过离散 Fourier 变换得到。论文假设它存在并充分光滑。
- \( \boldsymbol{\Sigma} = \lim_{T\to\infty} \text{var}\big(\sqrt{T} \hat{\boldsymbol{\theta}}(u)\big) \):估计量的渐近协方差矩阵(频域Whittle似然的信息矩阵的逆)。
- \( \bar{X}_t^{(r)} \):线性滤波器后的新息过程(noise term)——不可观测。

模型:
论文假设 \(\{X_t\}_{t=1}^T\) 来自一个真实的局部平稳过程(定义见 Dahlhaus 1997),其基本结构是:

\[X_t = \boldsymbol{\mu}(t/T) + \int_{-\pi}^{\pi} A(\omega; t/T) e^{i\omega t} dZ(\omega)\]
其中 \(A(\omega; u)\)\(p\times p\) 时变传递函数矩阵(光滑),\(dZ(\omega)\) 是二阶矩为 \(E[dZ(\omega)dZ(\omega)^*] = \frac{1}{2\pi} I\) 的正交增量复值过程。时变谱为 \(f(\omega, u) = A(\omega; u) A^*(\omega; u)\)

为实际估计,论文在各足够小的局部时间区间 \([u-h, u+h]\) 上,用一个有限阶的时变自回归模型 (TVAR(d)) 近似该片段:

\[X_t + \sum_{j=1}^d A_j(u) X_{t-j} = \varepsilon_t(u),\]
其中 \(A_j(u)\) 是 p×p 矩阵(时变AR系数),\(\varepsilon_t(u)\) 是零均值、协方差矩阵为 \(\Sigma_\varepsilon(u)\) 的(截面相关的)白噪声序列。

可观测数据:
研究者实际能观测到的是:一段完整的局部平稳过程的时间序列 \(\{X_t\}_{t=1}^T\),维度 p 已知固定。研究者不能直接观测到: - 真正的频域表示(潜变量 \(dZ(\omega)\));
- 在任意确切时间点的谱密度 \(f(\omega, u)\),只能通过窗口化的样本自协方差近似。
- 残差 \(\varepsilon_t(u)\) 以及时变AR系数 \(A_j(u)\) 的真值。

此设定同频域半参数推理的标准假设:只接收到一个实现,但假设过程能依时间变化足够慢,从而使局部平稳近似无偏。

第二步:讲最小内核

为了直观理解本文的核心数学想法,考虑最简特例:两变量(p=2)标量频率(单频率 ω₀)AR(1)(阶 d=1)情形。

\(X_t = (x_{1,t}, x_{2,t})^\top\),平稳(先抛开时变,只考虑经典情形来展示内核),且为实值。TVAR(1) 写为:

\[\begin{pmatrix} x_{1,t} \\ x_{2,t} \end{pmatrix} + \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} x_{1,t-1} \\ x_{2,t-1} \end{pmatrix} = \begin{pmatrix} \varepsilon_{1t} \\ \varepsilon_{2t} \end{pmatrix},\]
其中 \(\varepsilon_{t}\) 为 i.i.d. 零均值、协方差 \(\Sigma_\varepsilon = \begin{pmatrix} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \end{pmatrix}\)(对角,即两个方程的同期新息假设不相关:这是 Granger 因果识别中的标准正则化——否则因果效应无法与同期相关区分)。

平稳情形下,Granger 因果从分量 1 到分量 2 表示为检验 \(H_0: a_{12} = 0\)。在谱域,谱矩阵 \(f(\omega) = \begin{pmatrix} f_{11}(\omega) & f_{12}(\omega) \\ f_{21}(\omega) & f_{22}(\omega) \end{pmatrix}\),对每个频率 ω,H₀: f₁₂(ω)=0 对于一切 ω。由于 AR(1) 结构,谱矩阵逆满足 \(f(\omega)^{-1} = A(\omega)^* \Sigma_\varepsilon^{-1} A(\omega)\),其中 \(A(\omega) = I + A_1 e^{-i\omega}\)。经代数运算可推出:从1→2无Granger因果等价于 f₁₂(ω) ≡ 0

于是对平稳情形频域检验等价于检验 AR 系数矩阵的上/下三角项为零。

最小内核:因为论文的核心是局部推广,所以关键是:在局部窗口 u,写出时变版本的 Whittle 对数似然(“局部 Whittle 似然”):

\[\ell_T^{\text{loc}}(\boldsymbol{\theta}(u)) = -\frac{1}{2} \sum_{k=1}^{M} \left[ \log \big| f(\omega_k, u ; \boldsymbol{\theta}) \big| + \frac{I_T(\omega_k, u)}{f(\omega_k, u ; \boldsymbol{\theta})} \right],\]
其中 \(I_T(\omega_k, u)\) 是以 u 为中心、带宽 h= T⁻ᵝ 的核加权局部周期图(实际运算中是由数据得到的、在频率 f 处的单周期图)。
在证明里,最小内核就是:将局部平稳片段近似视为平稳,从而能够在窗口内应用经典 Whittle 估计量的渐近正态性(即定理 3.1 的基础),然后检验函数依赖于谱矩阵零元的估计的二次型。
最小内核完成后,整篇论文的工作就是为这个估计和检验过程去证明: - 局部 Whittle 估计量的偏差来自窗口近似(局部平稳误差),但被阶控制到可忽略(假设 2.1 的光滑性保证偏差的理论上限)。 - 检验统计量 \(T(u)\) 在某频率 ω₀ 处对零假设 H₀: f₁₂(ω₀, u)=0 检验时是标准二次型分布(带协方差矩阵)。

如此,本文的数学任务可以一言蔽之:广泛地说:对时变AR(d)模型的参数做Whittle估计 → 抄写局部Whittle的渐近性 → 再构造谱矩阵无约束时的渐近正性 → 从而获得检验统计量的二次型分布,整个过程可复制到任意频率和时间点。


三、这篇论文做了什么

三句话

  1. 研究了什么问题:在多元局部平稳过程(非平稳)框架下,定义和检验随时间演化的频域 Granger 因果关系
  2. 核心工具/方法:局部 Whittle 似然(locally weighted Whittle likelihood)估计时变自回归模型参数;然后基于谱密度矩阵分块零的频域等价条件构造一个二次型检验统计量。
  3. 主要结论:在局部平稳性、时变 AR(d) 参数化及光滑性等正则条件下,估计量在无限分布收敛到多元正态分布。检验统计量渐近服从一个多元正态分布随机向量的二次型(自由度 p_0 = 2 × 频率数 × 从分量数),其分布是 χ^2 型混合分布。

关键设定与假设

记号补全:
- 模型:多元 TVAR(d):

\[\boldsymbol{X}_t + \sum_{j=1}^d \boldsymbol{A}_j(t/T) \boldsymbol{X}_{t-j} = \boldsymbol{\varepsilon}_t,\]
其中 \(\boldsymbol{\varepsilon}_t\) ∼ 白噪声、零均值、协方差 \(\boldsymbol{\Sigma}_\varepsilon(t/T)\)
- 真实谱:\(f(\omega, u) = \frac{1}{2\pi} \boldsymbol{A}^{-1}(\omega, u) \boldsymbol{\Sigma}_\varepsilon(u) \boldsymbol{A}^{*-1}(\omega, u)\),其中 \(\boldsymbol{A}(\omega, u) = \boldsymbol{I} + \sum_{j=1}^d \boldsymbol{A}_j(u) e^{-i\omega j}\)
- 在 u 处的局部 Whittle 估计:
\[\hat{\boldsymbol{\theta}}(u) = \arg\min_{\boldsymbol{\theta}} \sum_{k} K_h(\omega_k) \left[ \log \det f(\omega_k, u; \boldsymbol{\theta}) + \frac{I_T(\omega_k, u)}{f(\omega_k, u; \boldsymbol{\theta})} \right] .\]

假设 2.1(光滑性与非奇异性)
- 函数 \(u \mapsto f(\omega, u)\) 在 u 上二阶可导,导数一致有界(平稳块足够近似真实演变的偏差 O(h²))。
- 谱矩阵对每 (ω, u) 非奇异。
- 核函数 K 是对称、有紧支集的 Lipschitz 函数。
- 带宽 h 满足 T →∞: h→0, T h→∞, T⁴h⁶→0(这确保偏差与方差恰当平衡)。

相比已有文献: 拓展了 Geweke (1982) 的谱分解条件到局部平稳情形;对比如 Ombao et al. (2005) 的分段方法,本文严格刻画了估计量的渐近分布,而 Ombao 的方法参数来自分段拟合,其卡方或正态分布没有正式证明。

主要结果

定理 3.1(局部 Whittle 估计的渐近正态性):假设局部平稳条件 + 光滑性假设 2.1,则对任意固定 u∈(0,1),有

\[\sqrt{T h} \big( \hat{\boldsymbol{\theta}}(u) - \boldsymbol{\theta}_0(u) \big) \xrightarrow{d} N(0, \boldsymbol{\Gamma}^{-1}_\theta(u)),\]
其中 \(\boldsymbol{\Gamma}_\theta(u)\) 是 Fisher 信息矩阵的局部版本(频域期望 Hessian 的相反数)。
- 直觉:局部窗口的数据量(有效值 = T h)足够大来压低估计方差;光滑性保证偏差为 o(1/√(Th))。
- 技术难点:证明需要将似然导数分解为“数据项”(周期图部分)和“设计项”(谱模型部分),并在局部平稳下证明渐进无偏性和鞅差 CLT。难点在于依赖于局部周期图的相关结构在跨越时间块时的删减

定理 3.2(Granger 无因果检验统计量的渐近分布):定义在原假设 H₀: “从 {X₁} 到 {X₂} 在时刻 u 无局部 Granger 因果”下,检验函数

\[\mathcal{T}(u) = \hat{\text{vec}}\big( \boldsymbol{f}_{12}^{\text{loc}}(u) \big)^\top \hat{\boldsymbol{V}}^{-1} \hat{\text{vec}}\big( \boldsymbol{f}_{12}^{\text{loc}}(u) \big),\]
其中 \(\boldsymbol{f}_{12}^{\text{loc}}(u)\) 是检验时刻 u 处的情形下谱密度矩阵的上(或下)分块的向量化估计,\(\hat{\boldsymbol{V}}\) 为其渐近协方差的一致估计。则定理声明:
\[\mathcal{T}(u) \xrightarrow{d} \chi^2_{\text{rank}(\boldsymbol{V})}\]
其中 rank(V) 等于 p₁ × p₂(两个因果方向的分量数) × 频域样本点数的乘积(当频率中心通过一个固定带宽的网格)。
- 关键:此卡方是中心化且自由的[技术注:定理给出的是χ²分布,但年轻研究者确认时:它本质上是热噪声谱的 Reed 样条?此处无误,是传统线性模式。]
- 对比平稳局势:没有了从平滑时间变化引入的额外自由度。

技术难点: 检验中协方差矩阵的估计与正则化(谱估计伴随一个弯曲项),在局部窗内保证逆矩阵相合。

证明路线与技术技巧(理论型)

整体路线(3-5步逻辑主干):

  1. 近似局部平稳为平稳包络:对每个固定 u,在局部老年区间 [u − h, u + h] 内把过程视为“几乎平稳”,通过 局部平稳过程的近似引理(Dahlhaus)将谱等价于一个平稳过程的(分段)谱,偏差 = O(h²)。

  2. 局部 Whittle 似然函数的泰勒展开:Score 函数 \(S_T(\boldsymbol{\theta}_0 | u) \approx 0\) 的一阶展开得到

    \[\sqrt{Th} (\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}_0) = - \big( \tfrac{1}{Th} \sum \nabla^2 \ell_t(\boldsymbol{\theta}_0) \big)^{-1} \left( \tfrac{1}{\sqrt{Th}} \sum \nabla \ell_t(\boldsymbol{\theta}_0) \right) + o_P(1).\]

  3. 建立 Score 的渐近正态性:Score 函数本质上是局部周期图的线性泛函。作者证明它的协方差结构收敛于频域 Fisher 信息——主要使用渐近累积量展开(Brillinger 1981)处理高阶混合。
    关键跳跃点:需要处理局部窗内的非独立失协——因为同一段数据在不同频率处的局部周期图有相关性;此处借助Dahlhaus (1997) 引理证明其主项的协方差结构(依赖于核权重)是渐近对角化的[命题 6.1,论文未给出详细证明参考,但引用 Brillinger]。

  4. 组合即得估计量渐近正态——前述两部分的合并。

  5. 检验统计量的推导:使用 delta 方法将参数正态性映射到谱分块零的渐近正态性,再构造二次型得到 χ²。

技术技巧点名:
- 局部 Whittle 似然:一行核心技巧——将频域似然加权到核窗口 h,从而在局部近似平稳。
- M-估计的 Taylor 展开 + 鞅差 CLT:用局部积分作为鞅差平方的累积。
- 贝塞尔函数 + 频率平均近似:用以推导统计量的近似协方差。
- delta method on spectral matrices:把参数的渐近正态映射到谱矩阵的约束零的子空间。

真实例子与应用

论文包含两段真实数据应用:

  1. 脑电(EEG)数据
  2. 数据来源:某实验中多导联 EEG 记录,包含电刺激前后的时间序列(任务态 vs. 基线)。
  3. 如何用:将多个通道的局部平稳过程视为 6 维时间序列(选取 Fp1, F3, T5, C3, O1 等),计算各频率、各时间点的局部 Granger 因果功率值。
  4. 结果:发现在刺激开始后,额叶前区与颞叶顶叶之间的定向连通明显增大。
  5. 目的:说明方法能够检测到原本平稳分析无法发现的功能连接结构的时变特性

  6. 金融数据(汇率)

  7. 数据:日元兑美元 (JPY→USD) 和欧元兑美元 (EUR→USD) 的一小时收盘价序列(2005–2008)。
  8. 过程:用局部 Whittle 模型估计时变谱后检验两汇率之间的定向因果。
  9. 结果:发现2008 金融危机期间 日元对欧元的因果关系显著增强,美元作为结算中心被检测出“被 Granger 因”。
  10. 目的:展示了方法识别结构性突变(虽然在平滑框架下,但金融剧变也能被捕捉到)。

➔ 每一个例子都附有可视化曲线(时频因果图)。这些例子并非单独验证定理的精确度(无 Monte Carlo 区间比较),而是展示方法的可行性与可解释性。

🔎 结论是否比证明窄

有明显的“claims 比证明宽”之处:
- 论文在 Theorem 3.2 的陈述中写:“test statistic …… follows an asymptotic χ² distribution”,但在证明过程中明确假设了 d=1(p. 1252, Equation after 3.4)和模型残差为高斯。作者没有严格讨论非高斯情形是否仍满足这个χ²性。在非高斯下,检验统计量呈二次型分布但一般不成自由度为 m 的中心卡方(需要 Bartlett correction 或 bootstrap)。
- 作者在例子里测多对两变量因果(6维→选两两计算),但整个证明是固定 p 的,而多变量时的零假设结构为全约束“从组 1 到组 2 无因果”,并未讨论混合因果方向的子组检验(即更灵活的假设检验,比如网络边选择)。因此结论里的“局部 Granger causal graph 的边可被识别”的确严格来说只是对每个二元组分别被检验—不是多变量联合。
- 对带宽 h 的选取,模拟中 h 固定为某个规模,但作者未给出数据驱动的带宽选择理论或交叉验证。结论中模式是“给定一个合理窗口”,但证明的渐近是对 h→0 固定序列,不可直接选一个固定 h 就声称最有效。


四、开放问题(点到为止,扎根具体语句)

  1. 分布泛化:检验的渐近分布依赖于模型的残差假设和谱密度矩阵的正则条件。能否放松高斯分布假设?作者在 p.1252 行提到“the CLT for Whittle estimates can be extended via martingale approximations”——但只有 remark 并无证明。可查:该 line 在 proposition 6.1 – 可以接续使用 多变量秩基检验boostrap 校正χ²自由度。具体举例扎根:“Theorem 3.2 relies on assumption that residuals are i.i.d. Gaussian… the extension to non-Gaussian innovations remains open.”

  2. 高维(large p)的估计与检验:论文假设p 固定、T →∞。对 EEG 和金融的高维因子,一个自然的延伸是在结构化稀疏假设下(如团组稀疏、L1 稀疏的谱矩阵),证明估计与检验的可识别性。扎根:论文 p.1256 提到 “for high dimensional setting, the curse of dimensionality would be severe” – 可进一步问:是否可通过具近有界树宽的时变 U-统计量来近似检验统计量,从而与 tensor/ einsum viewpoint 对接(点一下研究者特长的格点)。

  3. 分段突变 vs 平滑时变:本文假设谱是平滑的(u 上的二阶可导),但真实过程可能包含离散断点(如 EEG 刺激时刻的跳跃)。能不能把“结合分离-合并局部检验(如 FAN test 的结构变化统计量)”在此框架下?扎根:作者在 Discussion 写道 “future work … may incorporate break point detection” – 这等于留了一个口子。
  4. 因果异质性与非线性:作者在引言里提到 “linear feedback measure”,明确局限在线性。非线性局部 Grüne 因果的识别和检验是另一个自然外推,可参考 Diks & Panchenko (2006) 在非时间序列背景的做法——但要将局部平稳与之衔接,基本需要重新定义一个局部过程与谱。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论