• TOP
  • ブログ
  • 【後編】初心者のための時系列間の因果関係検出~らくらくノンパラメトリック・非定常・非線形・多変量に対...

ホットリンク 公式ブログ

データから見えてくる注目のトピックスをとりあげます。

テック

【後編】初心者のための時系列間の因果関係検出~らくらくノンパラメトリック・非定常・非線形・多変量に対応~

開発本部研究開発グループR&Dチームのサンティです。

今回のブログでは、エンジニア向けに時系列間の因果関係検出に関する論文「A Nonlinear Causality Estimator Based on Non-Parametric Multiplicative Regression」(以降[参考1])について前後編にわけてお伝えしています。本日はその後編です。サンプルコードは、githubにアップロードしてあります。

前編はこちらを参照ください。

 

実践編

ここからは実践編です。まず、悪名でよく知られる全く関係ない単位根過程の時系列同士の因果関係です。下記はそれぞれ単位根過程の時系列1000ポイントを二つ生成します。

In [20]:
# ここから必要なパッケージ
using DataFrames
using Gadfly
using Distributions
In [21]:
function dataset0()
    (cumsum(rand(Normal(),1000,1)),cumsum(rand(Normal(),1000,1)),nothing)
end
(X1,X2) = dataset0()
 
In [22]:
isCausal(X1,X2,d=3)
 
Cnpmr(X->Y): -0.0049280813575959545

 

Out[22]:
false

X1→X2の因果関係がありません。
 
In [23]:
isCausal(X2,X1,d=3)
 Cnpmr(X->Y): -0.0739435833211596
 
Out[23]:
false

X2→X1の因果関係もありません。とりあえず、基本中の基本、第一関門通過です。
 

次に[参考1]のdataset 1の検証です。2つの時系列は下記の式で生成します。 見ての通り、X2は独立で、X2からX1への因果性があります。しかも非線形です。

 

In [24]:
function dataset1(;N=1000)
    noise = rand(Normal(),N+1,2)
    x1 = zeros(N+1,1)
    x2 = zeros(N+1,1)
    for t = 2:(N+1)
        x2[t] = 0.6*x2[t-1] + noise[t,2]
        x1[t] = 0.8*x1[t-1] + 0.65*x2[t-1]^2 + noise[t,1]
    end
    (x1[2:end,:],x2[2:end,:],nothing)
end
(X1,X2) = dataset1()
widedf = DataFrame(x= collect(1:length(X1)), X1=X1[1:end],X2=X2[1:end])
longdf = stack(widedf,[:X1,:X2])
plot(longdf,ygroup="variable",x="x",y="value",Geom.subplot_grid(Geom.line))
 
 
Out[24]:
 

まずは、X1→X2の検証

 

In [25]:
isCausal(X1,X2,d=3)
 
Cnpmr(X->Y): -0.1306311369683256

Out[25]:
false

 値がマイナスなのでX1→X2の因果関係はありません。 次にX2→X1の検証
 
In [26]:
isCausal(X2,X1,d=3)
 Cnpmr(X->Y): 0.4527132517698962


Out[26]:

true
 

X2→X1の因果関係があります。しかも、Cnpmr値は[参考1]と一致範囲内です。 次に因果関係のある方向で予測要因の影響度を調べます。

 

In [27]:
sensitivity(X1,X2,d=3) #toY=X1, fromX=X2
 
Q_x1(t-1): 0.1270116553586665
Q_x1(t-2): 0.04839524137610466
Q_x1(t-3): 0.03347992177990148
 

見事[参考1]のdataset1の結果と一致で(t-1)のポイントが生成した関数の通り、飛び抜けた影響度となっています。

次にdataset2で多変量の事例を検証しましょう。

 

In [28]:
function dataset2(;N = 1000)
    noise = rand(Normal(),N+2,3)
    x1 = zeros(N+2,1)
    x2 = zeros(N+2,1)
    x3 = zeros(N+2,1)
    for t = 3:(N+2)
        x1[t] = 0.95*sqrt(2)*x1[t-1] - 0.9025*x1[t-2] + noise[t,1]
        x2[t] = -0.5*x1[t-1] + noise[t,2]
        x3[t] = 0.5*x3[t-1] - 0.5*x2[t-1] + noise[t,3]
    end
    (x1[3:end,:],x2[3:end,:],x3[3:end,:])
end
(X1,X2,X3) = dataset2()
widedf = DataFrame(x= collect(1:length(X1)), X1=X1[1:end],X2=X2[1:end], X3=X3[1:end])
longdf = stack(widedf,[:X1,:X2,:X3])
plot(longdf,ygroup="variable",x="x",y="value",Geom.subplot_grid(Geom.line))
 
Out[28]:
 
In [29]:
isCausal(X1,X2,d=3)
 
Cnpmr(X->Y): 0.671684662174414

Out[29]:
true

In [30]:
isCausal(X1,X3,d=3)
 
Cnpmr(X->Y): 0.1852700882985855

Out[30]:
true

In [31]:
isCausal(X2,X3,d=3)
Cnpmr(X->Y): 0.4915685081813289

Out[31]:
true
 

[参考1]のFIGURE1(D)(a)と一致しています。 次にX2の情報が分かる場合、X1→X3の因果性があるかどうかを検証(条件付き検証)します。

 

In [32]:
#conditional
isCausal(X1,X3,d=3,Z=X2)
Cnpmr(X->Y): -0.21581962879436334
 
 
Out[32]:
false
 

X2が分かる場合、X1→X3の因果関係がない結果になります。 これも[参考1]のFIGURE1(D)(b)と一致しています。生成式(10)の通りです。 影響度もチェックしましょう。

 

In [33]:
#sensitivity X1->X2
sensitivity(X2,X1,d=3)
 
Q_x1(t-1): 0.3535896399770177
Q_x1(t-2): 0.09410860094704274
Q_x1(t-3): 0.10359095901217734
 
 
In [34]:
#sensitivity X2->X3
sensitivity(X3,X2,d=3)
 
Q_x1(t-1): 0.34336832098832626
Q_x1(t-2): 0.140831503791619
Q_x1(t-3): 0.10430415578957795
 

[参考1]のパターンと一致して、(t-1)時点が飛び抜けて影響が大きいことが分かります。

次に、dataset3を検証しましょう。今度は線形と非線形の因果関係が混在しているケースです。

まずはデータの生成と可視化です。

 

In [35]:
    noise = rand(Normal(),N+1,3)
    x1 = zeros(N+1)
    x2 = zeros(N+1)
    x3 = zeros(N+1)
    for t = 2:(N+1)
        x1[t] = 3.4*x1[t-1]*(1-x1[t-1]^2)*exp(-x1[t-1]^2) + 0.4*noise[t,1]
        x2[t] = 3.4*x2[t-1]*(1-x2[t-1]^2)*exp(-x2[t-1]^2) +
                    0.5*x1[t-1]*x2[t-1]+0.4*noise[t,2]
        x3[t] = 3.4*x3[t-1]*(1-x3[t-1]^2)*exp(-x3[t-1]^2) +
                    0.3*x2[t-1]+0.5*x1[t-1]^2 + 0.4*noise[t,3]
    end
    (x1[2:end,:],x2[2:end,:],x3[2:end,:])
end
(X1,X2,X3) = dataset3()
widedf = DataFrame(x= collect(1:length(X1)), X1=X1[1:end],X2=X2[1:end], X3=X3[1:end])
longdf = stack(widedf,[:X1,:X2,:X3])
plot(longdf,ygroup="variable",x="x",y="value",Geom.subplot_grid(Geom.line))
 
 
Out[35]:
 

ランダムにしか見えませんね。

 

In [36]:
isCausal(X1,X2,d=3)
 
Cnpmr(X->Y): 0.051409672095286636

Out[36]:
true
 

X1→X2は有意です。

 

In [37]:
isCausal(X1,X3,d=3)
 
Cnpmr(X->Y): 0.03990845956605069

Out[37]:
true

In [38]:
isCausal(X2,X3,d=3)
 
Cnpmr(X->Y): 0.10255229695465998

Out[38]:
true
 

X1→X3, X2→X3も有意で、 結果は[参考1]のFIGURE2(D)(a)と一致しています。条件付きの場合も同様です。気になる方は是非実行してみてください。 次にdataset4を検証しましょう。今度はカオス時系列での検証です。

[参考1]の意図は時系列間のカップリングの強さしだいでライバルの手法より感度が高いことを見せたいだけなんです。ここではほんの一部でエッセンスだけしらべますのでカップリングの強さcを0.4にします。[参考1]によるとライバルのカーネルグレンジャー手法(K-GC)はc<0.5では因果関係が検出できないようです。

 

In [39]:
function dataset4(;c=0.4,N = 1000)
    x = zeros(N+2,1)
    y = zeros(N+2,1)
    for t = 3:(N+2)
        x[t] = 1.4 - x[t-1]^2 + 0.3*x[t-2]
        y[t] = 1.4 -
            (c*x[t-1] + (1-c)*(y[t-1]))*y[t-1] +
            0.1*y[t-2]
    end
    (x[3:end,:],y[3:end,:],nothing)
end
(X,Y) = dataset4()
widedf = DataFrame(x= collect(1:length(X)), X=X[1:end],Y=Y[1:end])
longdf = stack(widedf,[:X,:Y])
plot(longdf,ygroup="variable",x="x",y="value",Geom.subplot_grid(Geom.line))
 
 
 
Out[39]:
 
In [40]:
isCausal(X,Y,d=3)
 
Cnpmr(X->Y): 0.7581493700934339

Out[40]:
true
 

見事X→Yの因果関係があるという結果です。

 

In [41]:
isCausal(Y,X,d=3)
 
Cnpmr(X->Y): 0.007377674002540528

Out[41]:
true
 

Y→X因果関係は何回かやりますと有意であったり、有意でなかったりします。これは式上、フィードバック的な因果関係を持たしていないが、別の論文でこのようなカオス時系列について、なぜか因果関係があるようになってしまう解説があるようです。 いずれにしても因果関係の強さはX→Yの方向よりも低いと[参考1]の結果と一致しています。 検証の続きがありますが、長くなってしまいますので検証はここまでにしておきましょう。

 

考察

立役者はいうまでもなく、NPMRモデルにあります。本家の資料を見た限り利用シナリオは、モデル対象のデータが十分に多く、予測要因がお互いに影響しあって、その関係性がハッキリわかっていなくて、どうしようもない時に、最後のリソースとして利用するものとしています。ところで[参考1]では、NPMRモデルの適用可能な範囲が意外と広いと判明、めでたく因果関係検出に適していると報告されました。

実際にやってみたところ、おもちゃのようなテストケースは問題なく再現できました。あくまで、NPMRモデルの最適化をしなくても因果関係のパターンを正しく検出できる範囲で用意されたにすぎません。論文の第4節「Discussion」では、注意事項やベストプラクティス、さらには、課題となっているケースなどが念入りに語られています。初期の発見ではありますが、因果関係検出としての様々な可能性を見せています。アイディア次第でコスパのよいツールに拡張できそうな気がします。

[参考1]の手法が簡単で高度な知識が必要としませんでしたが、計算量はデータの長さの二乗オーダーでしたね。長い時系列の場合はちょっと時間がかかります。スケールしにくいかもしれませんが、処理工程の依存関係があまりなく、並列処理は簡単にできそうです。

最後に

前編後編全て読んでいただいた方、心から尊敬しております。弊社近くまで来られた際には、是非ランチでも行きましょう。

 

参考

  1. Nicolaou Nicoletta, Constandinou Timothy G.: "A Nonlinear Causality Estimator Based on Non-Parametric Multiplicative Regression", Frontiers in Neuroinformatics vol.10, 2016 https://doi.org/10.3389/fninf.2016.00019
  2. Ryan G. James, Nix Barnett, James P. Crutchfield: "Information Flows? A Critique of Transfer Entropies", Phys. Rev. Lett. 116, 238701 (2016)
  3. McCune, Bruce: "Non-parametric habitat models with automatic interactions", Journal of Vegetation Science 17: 819-830, 2006
 
▼ホットリンク では一緒に働く仲間を募集しています!
 プログラマー(募集人員4名)
 サーバサイドエンジニア(募集人員1名)