PCMCI: Causal Discovery in Multivariate Crypto Time Series

Jakob Runge（2018、2020）が開発した PCMCI は、多変量時系列のための因果発見アルゴリズムです。グラフィカルモデル学習用の PC アルゴリズムと、条件付け集合が共通ドライバー、間接経路、ソース変数の自己相関という3つを同時に無力化するよう構成された瞬間条件独立性（MCI）検定を組み合わせています。最後の要素は丁寧に読む価値があります。これこそが、PCMCI がモメンタムの強い系列で素朴なリード・ラグ検定を破綻させる偽陽性率を引き継がない理由です。

この記事では、2段階の構成を第一原理から分解し、tigramite ライブラリで実装し、この手法を取引の入力にする前に必要な実データ研究を整理します。

既存の内容に加わるもの

このブログではすでに相関を扱ってきました。暗号資産の連動は大部分が共通因子によるものです。ペア間のシグナル相関は、10ペアのブック全体の分散の65%を PC1 だけで吸収できることを示しています。ペアごとの高い相関は取引可能な関係ではありません。そのため統計的裁定取引とペア取引では共和分を使います。しかも依存性そのものも安定していません。DCC-GARCH 動的相関は、静的な標本相関が依存性の定常性、一定の限界ボラティリティ、方向間の対称性を仮定している理由を列挙していますが、いずれも成り立ちません。

これらの手法はいずれも有向グラフを生成しません。BTC、ETH、SOL がともに動き、その大部分を共通因子が説明していることは示せても、BTC が SOL を直接駆動しているのか、それとも関係が完全に ETH を介しているのかは示せません。方向を探る通常の手段はグレンジャー因果性ですが、標準的なグレンジャー検定は二変量であり、観測されていない共通ドライバーと真の因果性を分離できません。

ここには、見た目以上に重要な第三の失敗モードがあります。各暗号資産のリターン系列には正の自己相関があり、自己相関はそのような2系列間で行う依存性検定の統計量を膨らませます。つまり、互いに独立な自己相関過程でも、有限標本では有意な関係があるように見えます。相関の分析ではこれを標本サイズ補正として扱います（シグナル相関の記事で説明した effective_N）。因果発見ではこれは補正ではなく中心的な設計問題であり、まさに第2段階の条件付け集合が解決する問題です。

PCMCI アルゴリズム

PCMCI は2段階で動作します。第1段階では各変数の因果的な親の候補集合を特定し、第2段階では慎重に構成した条件付き独立性検定で各候補リンクを検定します。

第1段階：PC 安定条件選択（PC1）

第1段階の目的は、各変数 $X_t^j$ について真の因果的親を含む上位集合 $widehat{\mathcal{B}}(X_t^j)$ を見つけることです。この集合は後で第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$ の独立性を、$X_{t-\tau}^i$ 自身を除く、残った親のうち関連が最も強い $p$ 個を条件として検定します。

$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}$ は、$\widehat{\mathcal{B}}(X_t^j) \setminus \{X_{t-\tau}^i\}$ に含まれ、$X_t^j$ との関連が最も強い $p$ 個の変数の集合です。

「安定」とは、1回の反復 $p$ における削除の判断が前の反復の親集合に基づくことを意味します。これにより結果が処理順序に依存しません。再現性にとって重要な性質です。

この段階の有意水準 $\alpha_{\text{PC}}$ は意図的に緩く設定します（通常は0.2以上、または tigramite の pc_alpha=None オプションで自動決定）。ここでの目的は偽陽性を制御することではなく、真の親をすべて残しながら無関係な変数をできるだけ刈り込むことです。$\widehat{\mathcal{B}}$ に少数の偽陽性があっても許容できますが、偽陰性は許容できません。

第2段階：瞬間条件独立性（MCI）検定

すべての変数について推定された親集合 $\widehat{\mathcal{B}}(X_t^j)$ が得られたら、第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)$

核心となる革新は条件付け集合です。次の2つを同時に条件付けします。

1. ターゲットの親 $\widehat{\mathcal{B}}(X_t^j) \setminus \{X_{t-\tau}^i\}$：共通ドライバーと $X_t^j$ へのその他の間接経路の影響を取り除きます。
2. ソースの親 $\widehat{\mathcal{B}}(X_{t-\tau}^i)$：ソース変数の自己相関を取り除き、検定統計量の膨張を防ぎます。

両方の集合を同時に条件付けることで、MCI はすべての交絡因子、媒介変数、自己相関の影響を差し引いた、特定のラグ $\tau$ における $X_{t-\tau}^i$ から $X_t^j$ への直接的で瞬間的な因果効果を分離します。

検定統計量には有効な条件付き独立性検定を使えます。線形関係には偏相関が標準的な選択です。非線形依存には、$k$ 最近傍法で推定した条件付き相互情報量（CMI）や GPDC（ガウス過程距離相関）検定を使えます。

デフォルト検定としての偏相関

ほとんどの金融アプリケーションでは、偏相関検定（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})$

つまり偏相関は、条件変数 $\mathbf{Z}$ の線形影響を取り除いた後の $X$ と $Y$ の線形関係を測定します。ガウスデータにおける条件付き独立性の帰無仮説の下では、検定統計量は既知の分布に従うため、$p$ 値を正確に計算できます。

偏相関の値は $[-1, 1]$ に収まり、因果リンクの強さを自然に順位付けできます。両方の親集合で条件付けした、ラグ1の BTC とラグ0の ETH の偏相関が0.3なら、「特定されたすべての交絡因子と自己相関の影響を取り除くと、BTC リターンへの1標準偏差のショックが、1期間後の ETH リターンの0.3標準偏差の変化を予測する」という意味です。

検定結果から因果グラフ（DAG）へ

PCMCI の出力は2つの行列です。

• 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$ が検定間の任意の依存を扱える補正です。デフレート Sharpe と多重検定は、これを具体化する帰無分布の偽発見率を測定します。公開するグラフには、検定したリンク数と補正後の閾値を必ず併記してください。

得られる構造は時系列グラフ（TSG）です。これは次のような有向グラフです。

• ノードは特定の時間ラグにおける変数を表します。
• 有向エッジは対応するラグを持つ因果リンクを表します。
• 自己ループは自己回帰効果を表します。

この TSG は、ラグをまたいで集約した変数間の因果リンクの存在と方向だけを示す要約グラフに縮約できます。解釈にはこちらの方が実用的な場合が多いでしょう。

Tigramite による実装

Tigramite は PCMCI のリファレンス実装で、ドイツ航空宇宙センター（DLR）に所属する Jakob Runge のグループが開発・保守しています。データ処理用の整理された API、複数の条件付き独立性検定、いくつかの PCMCI 変種、組み込みの可視化機能を提供します。

インストール

pip install tigramite

構造回復の健全性チェック（研究結果ではない）

市場データで実行する前に、パイプラインが自分で埋め込んだ構造を見つけられることを確認してください。これは実装の単体テストであり、それ以上のものではありません。以下の真の構造はデータ生成ループに手で書き込んだものなので、回復できてもコードが動くことを示すだけで、暗号資産については何も語りません。そのように読んでください。

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

埋め込んだ構造は、重要な否定的結果も含めて回復されます。この標本では BTC と SOL が強く相関しているにもかかわらず、BTC -> SOL の直接エッジはありません。MCI はその関係を ETH という媒介変数だけに帰属させます。これは二重の条件付け集合が生み出すべき挙動であり、このブロックの目的はそれを確認することです。これは単体テストの合格であって、暗号資産についての発見ではありません。

実際の市場データを扱う

PCMCI には弱定常な入力が必要なので、価格ではなく標準化した対数リターンを渡します。ADF の仕組みと生の価格が不適切な理由は、統計的裁定取引とペア取引で説明しています。標準化リターン行列が得られたら、tigramite 固有の部分は次の2行です。

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

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

通常のリターン処理からそのまま引き継げない tigramite 固有の注意点が2つあります。

1. 欠損データ：Tigramite は欠損観測を扱うマスク配列をサポートします。dataframe.mask で欠損箇所を示してください。作り物の観測値が、それを含むすべての条件付け集合に伝播するため、前方補完はしないでください。
2. 頻度の整合：すべての系列を同じ時間グリッドにそろえる必要があります。暗号資産では取引所が同期した OHLCV データを提供するため、通常は簡単です。

この手法に必要な研究

ここまでの内容はすべて仕組みです。このブログで記事を公開できるのは、実データに適用して、答えが「安定したものは何もない」だった場合も含め、何が見つかったかを報告してからです。必要な実行内容は次のとおりです。

• 設定：指定した期間の実際の資産バスケットの時間足リターンに PCMCI+ を適用し、tau_max を領域知識で固定し、$N^2 \times \tau_{\max}$ 本の全リンクに Benjamini-Yekutieli 補正を施した alpha_level = 0.01 を使う。デフレート Sharpeの基準に従い、検定したリンク数と補正後の閾値を報告します。
• 報告する結果：補正後に残るエッジ、その偏相関の大きさ、アウト・オブ・サンプルで残るエッジの数。
• 少なくとも1つの安定性測定：これがなければ、この記事は暗号資産ティッカーを付けただけの tigramite README です。ローリングウィンドウ間でエッジがウィンドウからウィンドウへ持続する割合、tau_max と pc_alpha に対する回復グラフの感度、または同じデータでの ParCorr と CMIknn の不一致を測定します。

不安定なグラフも公開可能な成果であって、失敗ではありません。「主要5資産での PCMCI：因果グラフはローリングウィンドウを維持できない」は、このブログがすでに扱っている正直な否定的結果の流れに合いますし、動作するだけのチュートリアルより有用な記事です。

適切なパラメータを選ぶ

最大ラグ $\tau_{\max}$

このパラメータは因果発見の時間的な範囲を制限します。低すぎるとゆっくり伝播する効果を見逃し、高すぎると計算コストと多重検定の負担が増えます。

時間足の暗号資産データでは、$\tau_{\max} = 6$ から $24$ が妥当な範囲です。日次データでは $\tau_{\max} = 5$ から $10$ で、ほとんどのリード・ラグ関係を捉えられます。ここでは領域知識が重要です。資金調達率の効果が伝播するのに8時間かかると分かっているなら、$\tau_{\max} \geq 8$ に設定します。

条件選択の有意水準 $\alpha_{\text{PC}}$

pc_alpha=None にすると、tigramite は推奨されるデフォルトである赤池情報量規準を使ってこのパラメータを自動選択します。手動で設定するなら、0.1から0.4の値が適切です。低い値ほど第1段階が積極的になり（残す親が少なくなり）、計算コストは下がりますが、真の親を削除するリスクが高まります。

最終有意水準 $\alpha$

これは標準的な仮説検定の閾値です。探索的分析なら $\alpha = 0.05$ で問題ありません。取引判断に使うものなら $\alpha = 0.01$ 以下を使い、先に説明した依存性を考慮した補正を適用した後に使ってください。名目水準は実際に検定している水準ではありません。

PCMCI+：同時リンクを追加する

標準 PCMCI はラグ付き因果リンク（$\tau \geq 1$）だけを発見します。しかし、情報が数秒で資産間を伝わる暗号市場では、時間足のサンプリング頻度のため、多くの因果効果が同時（$\tau = 0$）に見えます。

PCMCI+（Runge、2020）は PCMCI を拡張し、ラグ付きと同時の因果リンクを両方発見します。同時リンクはデフォルトでは無向です（時間の順序だけでは同じ時点の原因と結果を区別できないため）が、一部は標準 PC アルゴリズムの向き付け規則（コライダー検出、非巡回性制約）で向きを付けられます。

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

時間足以下の頻度で暗号資産を扱う場合、多くのクロスアセット効果はサンプリング間隔より速く発生するため、PCMCI+ の方が標準 PCMCI より適切なことが多いでしょう。

非線形拡張

金融時系列は非線形依存を示すことが多くあります（ボラティリティ・クラスタリング、レジームに依存するリード・ラグ効果など）。Tigramite は複数の非線形条件付き独立性検定を提供します。

ガウス過程距離相関（GPDC）

from tigramite.independence_tests.gpdc import GPDC

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

GPDC はガウス過程回帰で条件付け集合の影響を取り除き、その後残差に距離相関を適用します。非線形効果の検出では 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$ 最近傍法で条件付き相互情報量を推定します。完全にノンパラメトリックで、任意の関数依存を検出できます。代償は計算コストと、統計的検出力を得るためにより多くのデータが必要なことです。

ほとんどの暗号資産取引では、まず ParCorr から始めてください。非線形の因果メカニズムを示す具体的な証拠があり、十分なデータ（通常 $T > 5000$）がある場合にだけ、非線形検定へ切り替えます。

因果グラフを組み込める場所

因果グラフが相関行列には与えられないもの、それが出次数です。20〜50資産で PCMCI を実行すると、出ていくエッジが多いノードは情報の先行者になります。その動きはブック全体に対する予測情報を持ち、対称な相関行列では構造上表せない方向性の主張です。因子分解は BTC が分散を支配していることを示しますが、出次数の順位は BTC が先行しているために支配的なのか、それとも BTC を含むすべてが同じマクロショックにさらされているためなのかを示します。

安定したグラフが回復できたと仮定すると、このブログの既存研究に接続できる場所は明確です。グラフトポロジーのローリングウィンドウによる変化は構造変化のシグナルであり、HMM によるレジーム検出のテーマです。未解決の問題は、因果グラフが状態モデルより早く遷移を検出できるかです。複数の取引所にまたがる単一資産に PCMCI を適用すれば、スマート注文ルーティングですでにミリ秒精度で測定した取引所間の価格先行性を扱えます。時間足のグラフはそれを言い換えるのではなく、上回る必要があります。また、資金調達率裁定取引で扱った、8時間の資金調達メカニズムが作る周期的なパープ・スポットリンクは、答えが分かっているテストケースです。PCMCI は独立に検証できる因果構造を回復するはずであり、発見ではなく検証の対象になります。

これらはいずれも発見ではありません。上で述べた研究を待つ仮説です。

制限事項と注意点

PCMCI は強力なツールですが、実務者が理解すべき重要な制限があります。

1. 因果十分性の仮定：PCMCI は関連するすべての変数が観測されていると仮定します。隠れた共通ドライバー（クジラの取引活動や未公開ニュースなど）が2つの観測資産に影響すると、それらの間に直接的な因果リンクがあると PCMCI が誤って報告する可能性があります。LPCMCI 変種は潜在的な交絡因子を許すことでこの問題に部分的に対処しますが、向き付けされたエッジは少なくなります。

2. 定常性の仮定：因果構造は分析ウィンドウ全体で一定だと仮定されます。実際には暗号市場のダイナミクスは急速に変化します。ローリングウィンドウ分析で構造変化を検出してください。

3. 線形と非線形：ParCorr では線形の因果効果だけを検出します。「BTC が5%を超えて下落したときだけ BTC が ETH を下落させる」のような非線形メカニズムは、線形検定では見えません。

4. サンプリング頻度は重要：サンプリング頻度より速く起きる因果効果は、標準 PCMCI では同時（$\tau = 0$）リンクとして現れ、方向が曖昧になる可能性があります。PCMCI+ を使い、より高頻度のデータを検討してください。

5. 多重検定：$N^2 \times \tau_{\max}$ 個の依存した検定を、上で説明したように補正します。依存構造が明らかな選択肢を排除する理由は、デフレート Sharpe と多重検定を参照してください。

6. 標本サイズの要件：信頼できる因果発見には十分なデータが必要です。目安として、$N < 10$ 変数で ParCorr を使うなら $T > 500$、非線形検定またはより大きな変数集合なら $T > 2000$ を目指します。

PCMCI と他手法の比較

グレンジャー因果性や移送エントロピーに対する PCMCI の主な利点は、多変量設定で自己相関と共通ドライバーを正しく考慮する MCI 検定です。数十の相関し自己相関も持つ資産が同時に相互作用する暗号市場は、まさにこの状況です。

結論

持ち帰るべき貢献は2段階の構成です。親集合を制限する疎な条件選択を行い、その後、条件付け集合にソースとターゲットの両方の親を含める MCI 検定を行います。後半こそが PCMCI をグレンジャーや移送エントロピーから分ける部分であり、暗号資産リターンのように自己相関のある系列で手法を実行可能にする部分です。

この手法に証拠は付属していません。回復されたグラフは情報の流れについての仮説であってシグナルではなく、このブログでは、指定されたバスケット、指定された期間、依存性を考慮した補正、そしてウィンドウを動かしたときにグラフが維持されるかを測る少なくとも1つの指標を用いて実データで実行されるまで、何の意味も持ちません。その実行が存在するまでは、ここにあるものをすべてツールとして扱ってください。

グラフが安定しないと分かったなら、それが記事になります。動くチュートリアルより価値のある否定的結果は、このブログでも初めてではありません。

参考文献

• 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