statsmodels でグレンジャー因果性の確認

前回の記事でグレンジャー因果性という用語を紹介したので、今回はPythonで実際に計算してみましょう。
参考: グレンジャー因果性

データは、VARモデルのパラメーター推定を行った時のように、乱数で生成します。
参考: Python(statsmodels)でVARモデルのパラメーター推定

コードの重複する部分は省略して、変数resultに、推定したモデルが格納されているとします。
(今回の記事用に実行し直したので、パラメーターは前回と違います。)

この推定結果のオブジェクトが、 test_causality というメソッドを持っており、
それを使って、グレンジャー因果性の検証ができます。
ドキュメント: statsmodels.tsa.vector_ar.var_model.VARResults.test_causality


statsmodelsにはこの他にも、二つの過程のグレンジャー因果性を確認するツールとして、
grangercausalitytests というのが用意されているようなのですが、今回はこちらは使いません。

それでは早速やってみましょう。


# 推定結果のVARモデルの確認。(以前の記事のコードで生成されている前提)
print(result.summary())
"""
  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Sun, 09, Feb, 2020
Time:                     16:37:19
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                    1.41978
Nobs:                     99.0000    HQIC:                   1.32614
Log likelihood:          -337.444    FPE:                    3.53438
AIC:                      1.26250    Det(Omega_mle):         3.32953
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.263791         0.705621            1.791           0.073
L1.y1         0.604660         0.124431            4.859           0.000
L1.y2         0.267836         0.145104            1.846           0.065
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const        -3.289405         0.292191          -11.258           0.000
L1.y1         0.601292         0.051526           11.670           0.000
L1.y2         0.263057         0.060086            4.378           0.000
========================================================================

Correlation matrix of residuals
            y1        y2
y1    1.000000  0.596470
y2    0.596470  1.000000
"""

# y2 から y1 へのグレンジャー因果性の確認
test_result = result.test_causality(caused="y1", causing="y2")
print(test_result.summary())
"""
Granger causality F-test. H_0: y2 does not Granger-cause y1. Conclusion: fail to reject H_0 at 5% significance level.
==============================================
Test statistic Critical value p-value    df
----------------------------------------------
         3.407          3.890   0.066 (1, 192)
----------------------------------------------
"""

# y1 から y2 へのグレンジャー因果性の確認
test_result = result.test_causality(caused="y2", causing="y1")
print(test_result.summary())
"""
Granger causality F-test. H_0: y1 does not Granger-cause y2. Conclusion: reject H_0 at 5% significance level.
==============================================
Test statistic Critical value p-value    df
----------------------------------------------
         136.2          3.890   0.000 (1, 192)
----------------------------------------------
"""

当然データによって結果が変わる(乱数で生成しているので、モデルが同じでも実行するたびに変わります。)のですが、
y1からy2へのグレンジャー因果性がないという過程は棄却されているので、y1はy2に影響しているようです。
y2からy1の向きは、そうはならなかったですね。

さて、上のコード例のように、 test_causality というメソッドに、
caused と causing の引数で変数名を渡してあげるだけで、テスト結果を返してくれました。
例では変数名を指定しましたが、 0 や 1など、何番目の変数なのかを指定してもOKです。
これが2変数ではなくもっと多くの変数のVARモデルで、多変数から多変数へのグレンジャー因果性を調べる場合は、
causedやcausingをリストで指定すれば、同じようにテストできます。

さて、 test_causality は CausalityTestResults というオブジェクトを返してくれますが、
上の例のように、 summary() というメソッドを持っているので、人間の目で見る場合は、これをprintするのが良いでしょう。

p値などのプログラムで使いたい場合は、それぞれ変数に持っているのでそちらを使えます。


# 帰無仮説
print(test_result.h0)
# H_0: y1 does not Granger-cause y2

# 統計量
print(test_result.test_statistic)
# 136.18235785962648

# p値
print(test_result.pvalue)
# 3.9579973656016236e-24

# 結論
print(test_result.conclusion_str)
# Conclusion: reject H_0

グレンジャー因果性

ここ数回VARモデルの話を書いていますが、VARモデルの目的の一つは複数の時系列データ間の関連を調べることです。
普通は、各データにその背景があって、それらに関する理論からあの値がこうなったらどういう影響があってこっちの値がどうなる、と調べるのですが、
データだけから因果性の有無を判別することができると非常に便利です。
そのような考えを元に、 Grangerが提案したのがグレンジャー因果性です。

今回の記事も、沖本先生の経済ファイナンスデータの軽量時系列分析を元に各定義を紹介していきます。

定義(グレンジャー因果性)
現在と過去の$x$の値だけに基づいた将来の$x$の予測と、現在と過去の$x$と$y$の値に基づいた将来の$x$の予測を比較して、
後者のMSEの方が小さくなる場合、$y_t$から$x_t$へのグレンジャー因果性(Granger causality)が存在すると言われる。

これは、$x,y$の2つの過程に対する定義ですが、もっと多くの変数が存在する場合に拡張できます。
それが次の一般的なグレンジャー因果性です。

定義(一般的なグレンジャー因果性)
$\mathbf{x}_t$と$\mathbf{y}_t$をベクトル過程とする。また、$\mathbf{x}$と$\mathbf{y}$の現在と過去の値を含む、時点$t$において利用可能な情報の集合を
$\boldsymbol{\Omega}_t$とし、$\boldsymbol{\Omega}_t$から現在と過去の$\mathbf{y}$を取り除いたものを、$\tilde{\boldsymbol{\Omega}}_t$とする。
この時、$\tilde{\boldsymbol{\Omega}}_t$に基づいた将来の$\mathbf{x}$の予測と、$\boldsymbol{\Omega}_t$に基づいた将来の$\mathbf{x}$の予測を比較して、
後者のMSEのほうが小さくなる場合、$\mathbf{y}_t$から$\mathbf{x}_t$へのグレンジャー因果性が存在すると言われる。
ここで、MSEの大小は行列の意味での大小であることに注意されたい。

沖本本によると、グレンジャー因果性が頻繁に用いられるようになったたった一つの理由は、VARの枠組みでは、
グレンジャー因果性の有無が比較的簡単に明確な形で検証できるからだそうです。

2変量VAR(2)モデルを用いた説明が紹介されているのでみていきましょう。
モデルを具体的に書くと次のようになります。
$$
\left\{\begin{matrix}\
y_{1t}=c_1+\phi_{11}^{(1)}y_{1,t-1}+\phi_{12}^{(1)}y_{2,t-1}+\phi_{11}^{(2)}y_{1,t-2}+\phi_{12}^{(2)}y_{2,t-2}+\varepsilon_{1t},\\\
y_{2t}=c_2+\phi_{21}^{(1)}y_{1,t-1}+\phi_{22}^{(1)}y_{2,t-1}+\phi_{21}^{(2)}y_{1,t-2}+\phi_{22}^{(2)}y_{2,t-2}+\varepsilon_{2t}\
\end{matrix}\
\right.\
$$

ここで、$y_{2t}$から、$y_{1t}$へのグレンジャーレンジャー因果性が”存在しない”ということjは、
$\phi_{12}^{(1)}=\phi_{12}^{(2)}y_{2,t-2}=0$ということ同値です。
一般のVAR(p)において、$y_{2t}$から、$y_{1t}$へのグレンジャー因果性が存在しないとは、VARの$y_1$の式において、$y_2$に関係する係数が全部$0$ということになります。
その結果、VARの枠組みだと、F検定を用いてグレンジャー因果性を検定できます。
(複数の回帰係数をまとめて検定するので、t検定ではなくF検定だそうです。この辺り、自分も理解が浅いので別途勉強し直そうと思います。)

手順としては次のようになります。
$H_0: \phi_{12}^{(1)}=\phi_{12}^{(2)}y_{2,t-2}=0$ を検定すれば良いので、まずは、
$$
y_{1t}=c_1+\phi_{11}^{(1)}y_{1,t-1}+\phi_{12}^{(1)}y_{2,t-1}+\phi_{11}^{(2)}y_{1,t-2}+\phi_{12}^{(2)}y_{2,t-2}+\varepsilon_{1t}
$$
をOLSで推定して、残渣平方和を$SSR_1$とします。
次に制約を課したモデル、
$$
y_{1t}=c_1+\phi_{11}^{(1)}y_{1,t-1}+\phi_{11}^{(2)}y_{1,t-2}+\varepsilon_{1t}
$$
をOLSで推定して、その残差平方和を$SSR_0$とします。
この時、F統計量が次の式で定義されます。
$$
F\equiv \frac{(SSR_0-SSR_1)/2}{SSR_1/(T-5)}
$$
そうすると、$2F$は漸近的に$\chi^2(2)$に従うことが知られています。
なので、$2F$の値を$\chi^2(2)$の$95%$店と比較して、$2F$の方が大きければ、$y_2$から$y_1$へのグレンジャー因果性が存在しないという
帰無仮説を棄却して、$y_2$は$y_1$の将来を予測するのに有用であると結論づけられます。

これは1つの過程から別の1つの過程への影響を調べるものですが、もっと一般には次の手順になります。
(ここも沖本先生の本の丸写し)

n変量VAR(p)モデルにおけるグレンジャー因果性検定の手順
(1) VARモデルにおける$y_{kt}$のモデルをOLSで推定し、その残差平方和を$SSR_1$とする。
(2) VARモデルにおける$y_{kt}$のモデルに制約を課したモデルをOLSで推定し、その残差平方和を$SSR_0$とする。
(3) F統計量を
$$
F\equiv \frac{(SSR_0-SSR_1)/r}{SSR_1/(T-np-1)}
$$
で計算する。ここで、$r$はグレンジャー因果性検定に必要な制約の数である。
(4) $rF$を$\chi^2(r)$の$95%$点と比較し、$rF$のほうが大きければ、ある変数(群)から$y_{kt}$へのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論する。

あら?帰無仮説が棄却できなかった時に、存在しないって結論していいのかな?ここの記述はちょっと怪しいですね。

最後に、グレンジャー因果性の長所と短所がまとめられています。(沖本本で株式市場の例が紹介された後。)
長所は定義が明確であることt、データから容易に検定できることです。

短所は、まずはグレンジャー因果性が通常の因果性と異なっていることです。
グレンジャー因果性は因果性が存在する必要条件ではありますが、十分条件でありません。
なので、グレンジャー因果性があったからといって因果性があるとは結論づけられません。
また、グレンジャー因果性の方向と、通常の因果性の方向が同じになると限らず、
極端な例では逆向きだけ存在することもあるそうです。
そのため、何か別の理論等により因果関係の方向がはっきりしている時以外は、
定義通りに、予測に有用かどうかだけの観点で解釈することが良さそうです。

また、グレンジャー因果性は訂正的概念で、関係の強さは測れません。
(F統計量が大きいほど関係が強いみたいな考えは誤りということですね。)

Python(statsmodels)でVARモデルのパラメーター推定

今回は2変量のVARモデルにしたがうデータから、Pythonを使ってモデルの係数を推定する方法を紹介します。
答え合わせをしたいので、使うデータは特定のVARモデルから生成します。

モデルはいつもの沖本先生の本のP102ページの演習問題から拝借しました。

$$
\left\{\begin{matrix}\
y_{1t} &=& 2+ 0.5 y_{1,t-1}+ 0.4 y_{2,t-1}+\varepsilon_{1t},\\\
y_{2t} &=& -3 + 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\mathbf{\Sigma})
$$
$$
\Sigma = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

まずはデータを生成します。


import numpy as np
import matplotlib.pyplot as plt

mean = np.array([0, 0])
cov = np.array([[4, 1.2], [1.2, 1]])

# 少し長めの系列でデータを生成する。
data = np.zeros([110, 2])
epsilons = np.random.multivariate_normal(mean, cov, size=110)
print(np.cov(epsilons, rowvar=False))
"""
[[3.96508657 1.1917421 ]
 [1.1917421  1.01758275]]
"""

for i in range(1, 110):
    data[i, 0] = 2 + 0.5*data[i-1, 0] + 0.4*data[i-1, 1] + epsilons[i, 0]
    data[i, 1] = -3 + 0.6*data[i-1, 0] + 0.3*data[i-1, 1] + epsilons[i, 1]

# 初めの10項を切り捨てる
data = data[10:]
# data.shape == (100, 2)

fig = plt.figure(figsize=(12, 5), facecolor="w")
ax = fig.add_subplot(111)
ax.plot(data[:, 0], label="$y_1$")
ax.plot(data[:, 1], label="$y_2$")
ax.set_title("サンプルデータ")
plt.legend()
plt.show()

生成されたデータをプロットしたのがこちら。

さて、いよいよ statsmodelsでVARモデルのパラメーター推定をやりますが、
方法は非常に簡単で、時系列モデルのapiからVARモデルを呼び出してデータを渡し、学習するだけです。

ドキュメント:
Vector Autoregressions tsa.vector_ar


from statsmodels.tsa.api import VAR
model = VAR(data)
result = model.fit(maxlags=1)
print(result.summary())
"""
  Summary of Regression Results
==================================
Model:                         VAR
Method:                        OLS
Date:           Mon, 03, Feb, 2020
Time:                     01:10:12
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                    1.23044
Nobs:                     99.0000    HQIC:                   1.13680
Log likelihood:          -328.071    FPE:                    2.92472
AIC:                      1.07316    Det(Omega_mle):         2.75521
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.879511         0.616917            3.047           0.002
L1.y1         0.464346         0.124135            3.741           0.000
L1.y2         0.462511         0.132200            3.499           0.000
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const        -2.683660         0.305923           -8.772           0.000
L1.y1         0.561468         0.061557            9.121           0.000
L1.y2         0.385671         0.065556            5.883           0.000
========================================================================

Correlation matrix of residuals
            y1        y2
y1    1.000000  0.583971
y2    0.583971  1.000000
"""

今回はモデルの形が事前にわかっているので、 maxlags=1 としましたが、 ここを少し大きめの値にしておけば自動的にlagの値を選んでくれます。
選び方は 引数”ic”で指定できます。 (未指定ならAIC)

ic : {‘aic’, ‘fpe’, ‘hqic’, ‘bic’, None}
Information criterion to use for VAR order selection.
aic : Akaike
fpe : Final prediction error
hqic : Hannan-Quinn
bic : Bayesian a.k.a. Schwarz

さて、結果の係数ですが、summaryを見ると大体次の値になっています。(四捨五入しました)
$$
\left\{\begin{matrix}\
y_{1t} &=& 1.89+ 0.46 y_{1,t-1}+ 0.46 y_{2,t-1}\\\
y_{2t} &=& -2.68 + 0.56 y_{1,t-1}+ 0.39 y_{2,t-1}\
\end{matrix}\
\right.\
$$
とても正確というわけではありませんが、概ね正しく推定できていそうです。

Pythonで多変量正規分布に従う乱数を生成する

ベクトル自己回帰のダミーデータを生成するために、多変量正規分布に従う乱数が必要なので、
Pythonで生成する方法を紹介します。

numpyとscipyにそれぞれ用意されています。
同じ名前の関数だったので、どちらかの実装をもう一方がラップしているのかと思っていたのですが、
引数の微妙な違いなどあり、どうやら個別に実装されているようです。

ドキュメントはそれぞれ次のページにあります。

numpy
numpy.random.multivariate_normal
(この記事の numpy は version 1.16を使っています。 numpy 1.17.0 のリリースノートを見ると、random moduleに変更が加えられており、どうやらこの関数にも影響が出てるようなのでご注意ください。)

scipy
scipy.stats.multivariate_normal

さて、実際に期待値と分散共分散行列を指定してそれぞれ乱数を生成してみましょう。


import numpy as np
from scipy.stats import multivariate_normal

# 期待値と分散共分散行列の準備
mean = np.array([3, 5])
cov = np.array([[4, -1.2], [-1.2, 1]])

# numpy を用いた生成
data_1 = np.random.multivariate_normal(mean, cov, size=200)

# データ型の確認
print(data_1.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_1, axis=0))
# [3.00496708 4.94669956]

# 分散共分散の確認
print(np.cov(data_1, rowvar=False))
"""
[[ 3.86542859 -1.31389501]
 [-1.31389501  0.93002097]]
"""

# scipyで生成する方法
data_2 = multivariate_normal(mean, cov).rvs(size=200)

# データ型の確認
print(data_2.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_2, axis=0))
# [2.81459692 5.10444347]

# 分散共分散の確認
print(np.cov(data_2, rowvar=False))
"""
[[ 4.46151626 -1.28084696]
 [-1.28084696  1.06831954]]
"""

それぞれきちんと生成できたようです。

分散共分散行列の正定値性のバリデーションなど細かなオプションを持っていますが、
あまり使う機会はなさそうです。
(きちんと行う場合も、事前に固有値を求めて確認しておけば大丈夫だと思います。)

ベクトルホワイトノイズ

前回の記事で弱定常過程の定義をベクトルに拡張しましたが、今回はホワイトノイズをベクトルに拡張した
ベクトルホワイトノイズを定義します。
参考: 定常過程の例としてのiid過程とホワイトノイズ

ベクトル過程 $\boldsymbol{\varepsilon}_t$ が全ての時点 $t$ において次の二つの式を満たすとします。
$$
\begin{align}
E(\boldsymbol{\varepsilon}_t) &= \mathbf{0}&\\
E(\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_{t-k}^{\top}) &=
\left\{
\begin{matrix}
\boldsymbol{\Sigma},& k=0\\\
\mathbf{0},& k\neq 0 \
\end{matrix}
\right.
\end{align}
$$

この時、 $\boldsymbol{\varepsilon}_t$ をベクトルホワイトノイズと呼びます。
このベクトルホワイトノイズは弱定常であり、自己相関を持たないことはすぐわかります。
$k\neq 0$の時、$k$次の自己共分散行列が$\mathbf{0}$ですから。

沖本先生の本で注意されている通り、同時点での各変数は相関を持つことができます。

ただ、一点、次の記載があるのには少し引っかかります。(初版第13刷 P.76)

$\boldsymbol{\Sigma}$ は$n\times n$正定値行列であり、対角行列である必要はないことには注意が必要である.

$n\times n$正定値行列であることは確かなのでしが、これ、必ず対角行列になるような気がします。
というのも、$\boldsymbol{\Gamma}_k = \boldsymbol{\Gamma}^{\top}_{-k}$に、$k=0$ を代入したら、
$\boldsymbol{\Gamma}_0 = \boldsymbol{\Gamma}^{\top}_{0}$ となり、これは対角行列ですよね。

僕が何か勘違いしているのかもしれませんがおそらく誤植か何かのように思います。
(出版社サイトの今日時点の正誤表には入ってないようです)

弱定常ベクトル過程

ベクトル自己回帰モデルについて理解していくにあたって、弱定常過程やホワイトノイズの概念をベクトルに拡張したものを考える必要があるので、
ここで各用語の定義を整理しておきます。(前回の記事と順番前後してしまいました。)
出典は今回も沖本先生の本です。

また、1変数の場合については次の記事をご参照ください。
参考: 確率過程の弱定常性と強定常性

まず、 $\mathbf{y}_t = (y_{1t}, y_{2t}, \ldots, y_{nt})^{\top}$ は $n\times1$の行列(列ベクトル)としておきます。
そして、これの期待値ベクトルは普通に要素ごとの期待値をとって、次のように定義されます。
$$
E(\mathbf{y}_t) = (E(y_{1t}), \ldots, E(y_{nt}))^{\top}.
$$

次に、自己共分散を拡張します。ここからがポイントで、$y_{it}$に対して$y_{i,t-k}$との共分散を考えるだけではなく、
$\mathbf{y}_t$の他の要素の過去の値$y_{j,t-k}$との共分散も考慮します。そのため、結果は次のように行列になります。

そのため、$k$次自己共分散行列は次の$n\times n$行列で定義されます。

$$
\begin{align}
Cov(\mathbf{y}_t, \mathbf{y}_{t-k}) &= [Cov(y_{it},y_{j,t-k})]_{ij}\\
&=\left(\
\begin{matrix}\
Cov(y_{1t},y_{1,t-k}) & \cdots & Cov(y_{1t},y_{n,t-k}) \\\
\vdots & \ddots & \vdots \\\
Cov(y_{nt},y_{1,t-k}) & \cdots & Cov(y_{nt},y_{n,t-k}) \
\end{matrix}\
\right)
\end{align}.
$$

この行列の対角成分は、各変数の$k$次自己共分散に等しくなります。
また、$k=0$の時、$k$次自己共分散行列は、分散共分散行列に等しくなり、対称行列になりますが、
$k\neq 0$の時は、一般には対称行列にならないので注意が必要です。

期待値も$k$次自己共分散行列も定義式の中に$t$が入っているので、
これらは$t$の関数なのですが、この値がどちらも$t$に依存しなかった時、
1変数の場合と同じように、ベクトル過程は弱定常と言われます。

弱定常性を仮定する時は、期待値を$\boldsymbol{\mu}$、$k$次自己共分散行列を$\boldsymbol{\Gamma}_k$と書きます。
なお、$\boldsymbol{\Gamma}_k = \boldsymbol{\Gamma}^{\top}_{-k}$が成り立ちます。(転地を等号が成り立たない点に注意が必要です。)

さて、自己共分散行列は、異なる時点での各変数の関係を表したものですが、各要素の値の大小は単位に依存してしまい、
これだけ見ても関係が強いのか弱いのかわかりません。
そこで、自己共分散行列を標準化した自己相関行列を考えます。
それは次の式で定義されます。
$$
\boldsymbol{\rho}_k = Corr(\mathbf{y}_t, \mathbf{y}_{t-k}) = [Corr(y_{it},y_{j,t-k})]_{ij}.
$$

また、$\mathbf{D} = diag(Var(y_1), \ldots, Var(y_n))$、($\mathbf{y}_t$の分散を対角成分に持つ対角行列) とおくと、
次のように書けます。
$$
\boldsymbol{\rho}_k = \mathbf{D}^{-1/2} \boldsymbol{\Gamma}_k \mathbf{D}^{-1/2}.
$$

自己相関行列も自己共分散行列と同じように、対角成分は各変数の自己相関になり、
$\boldsymbol{\rho}_k = \boldsymbol{\rho}^{\top}_{-k}$ が成立します。

長くなってきたので一旦記事を切ります。(ベクトルホワイトノイズは次回。)

ベクトル自己回帰モデル

久しぶりに時系列分析の話題です。
今回はベクトル自己回帰(VAR)モデルを紹介します。
参考書はいつもの沖本先生の経済・ファイナンスデータの計量時系列分析です。
(VARモデルは第4章)

ベクトル自己回帰(VAR)モデルは、自己回帰(AR)モデルをベクトルに一般化したものです。
自然数$p$に対して、ベクトル $\mathbf{y}_t$ を定数(のベクトル)と、自身の$p$期以内の過去の値に回帰したモデルになります。
つまり、次のようなモデルになります。
$$
\mathbf{y}_t = \mathbf{c}+\mathbf{\Phi}_1\mathbf{y}_{t-1}+\cdots+\mathbf{\Phi}_p\mathbf{y}_{t-p}+\mathbf{\varepsilon}_t,
\ \ \ \mathbf{\varepsilon}_t\sim W.N.(\mathbf{\Sigma})
$$
($\mathbf{\varepsilon}$がなぜか太字にならない。自分の環境だけかな。)

ここで、$\mathbf{c}$ は $n\times 1$定数ベクトルで、$\mathbf{\Phi}_i$は$n\times n$行列です。
$\mathbf{\varepsilon}_t$はベクトルホワイトノイズと呼ばれるものです。
(単純に各要素がホワイトノイズのベクトルというわけではありません。このブログでは新出語なので後日別の記事で取り上げます。)

2変量のVAR(1)モデル、要するに実質的なVARモデルとしては一番小さいモデルを具体的に表記すると次のようになります。

$$\left\{\begin{matrix}\
y_{1t} &=&c_1+\phi_{11}y_{1,t-1}+\phi_{12}y_{2,t-1}+\varepsilon_{1t},\\\
y_{2t} &=&c_2+\phi_{21}y_{1,t-1}+\phi_{22}y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\mathbf{\Sigma})$$

$$
\Sigma = \left(\begin{matrix}\
\sigma_1^2 & \rho\sigma_1\sigma_2\\\
\rho\sigma_1\sigma_2 & \sigma_2^2\
\end{matrix}\right)\
$$

さて、具体的にモデルの形を書き出したのでパラメーターの数を数えてみたいと思います。
回帰式の係数(4個)と定数(2個)で計6個、さらに分散共分散行列部分が3個で、合計9個のパラメーターを持っていることがわかります。

一般的に$n$変数のVAR(p)モデルの場合、
回帰式の定数が$n$個、係数が$n^2p$個で合計$n(np+1)$個、分散共分散行列部分が$n(n+1)/2$個のパラメーターを持ちます。
そのため、比較的大きなモデルになってしまいます。

自分の経験でも、モデルに組み込みたい変数の値の数($n$)は10個も20個もあって、ラグ($p$)は3期くらい前まで見たい、
ということがありましたが、集められたデーターが全然足りず、推定が全然できないということがありました。
(そのときは結局変数をかなり絞り込んだモデルをいくつも作って個別に検証しました。)

さて、ベクトル自己回帰モデルを使う目的ですが、主に二つあります。
一つは、予測の精度向上です。
自己回帰モデルでは、ある値が自身の過去の値で回帰できることこを仮定して予測しますが、
ある値が自分の過去の値以外の外部要因に一切影響を受けない、ということは実用上あまりなく、
大抵は複数の値が相互に関係しているので、複数の種類の値を用いて予測することにより精度の向上が期待できます。

もう一つは、複数の変数間の関係を分析することです。
モデルの形からも明らかなように、それぞれの変数が他の変数にどれだけ寄与しているかをみることができます。

pandasで縦横変換(pivot_table)

前回の更新でPrestoでデータの縦横変換をする方法を紹介しましたが、
クエリで処理を完結させる必要がないときは、一旦pandasのデータフレームに格納してから処理をするのも便利です。

その場合、 pandas.pivot_table を使います。

使い方は簡単で、pd.pivot_tableに、変換したいデータフレーム、
列にするカラム、行にするカラム、集計する値、集計に使う関数を指定するだけです。
fill_value引数で欠損値を埋めるなどの細かい設定もできます。
ドキュメントの例を使ってやってみます。


import pandas as pd
# データ作成
df = pd.DataFrame(
    {
        "A": ["foo", "foo", "foo", "foo", "foo", "bar", "bar", "bar", "bar"],
        "B": ["one", "one", "one", "two", "two", "one", "one", "two", "two"],
        "C": ["small", "large", "large", "small",
              "small", "large", "small", "small", "large"],
        "D": [1, 2, 2, 3, 3, 4, 5, 6, 7],
        "E": [2, 4, 5, 5, 6, 6, 8, 9, 9]
    }
)
print(df)
"""
     A    B      C  D  E
0  foo  one  small  1  2
1  foo  one  large  2  4
2  foo  one  large  2  5
3  foo  two  small  3  5
4  foo  two  small  3  6
5  bar  one  large  4  6
6  bar  one  small  5  8
7  bar  two  small  6  9
8  bar  two  large  7  9
"""

table_0 = pd.pivot_table(
                df,
                values="D",
                index="A",
                columns="C",
                aggfunc="sum",
                fill_value=0,
        )
print(table_0)
"""
C    large  small
A
bar     11     11
foo      4      7
"""

# 行や列、集計関数は配列で複数指定することもできる
table_1 = pd.pivot_table(
                df,
                values="D",
                index=["A", "B"],
                columns="C",
                aggfunc=["sum", "count"],
                fill_value=0,
        )
print(table_1)
"""
          sum       count
C       large small large small
A   B
bar one     4     5     1     1
    two     7     6     1     1
foo one     4     1     2     1
    two     0     6     0     2
"""

PrestoのMap型を使った縦横変換

前回の記事で紹介したPrestoのMap型ですが、これを使うとデータの縦横変換(ピボット)がスマートに行えます。

参考: ピボットテーブル&チャート

上記リンク先のトレジャーデータの記事中の画像のテーブルを例に説明します。

縦持ちのテーブル (vtable)
uid key value
101 c1 11
101 c2 12
101 c3 13
102 c1 21
102 c2 22
102 c3 23

このようなテーブルを次のように変換したいとします。

横持ちのテーブル
uid c1 c2 c3
101 11 12 13
102 21 22 23

この変換は、MAP_AGGを使って、key列とvalue列の対応のMapを一旦作成し、
それぞれのMapから各Keyの値を取り出すことで実現できます。
具体的にクエリにしたのが次です。


SELECT
    uid,
    kv['c1'] AS c1,
    kv['c2'] AS c2,
    kv['c3'] AS c3
FROM (
    SELECT
        uid,
        MAP_AGG(key, value) AS kv
    FROM
        vtable
    GROUP BY
        uid
) AS t

個人的には服問い合わせはあまり好きではなく、PrestoではWITHを使って書きたいので次のようにすることが多いです。


WITH
    t AS (
        SELECT
            uid,
            MAP_AGG(key, value) AS kv
        FROM
            vtable
        GROUP BY
            uid
    )
SELECT
    uid,
    kv['c1'] AS c1,
    kv['c2'] AS c2,
    kv['c3'] AS c3
FROM
    t

Prestoではない、(MAP型のない)通常のSQLでは次のように書かないといけないのですが、
上記のMap型を使ったものの方が随分すっきりかけているように見えます。
(MAP_AGGを知らない人には読めないのが難点ですが)


SELECT
    uid,
    MAX(
        CASE WHEN key = 'c1' THEN
            value
        ELSE
            NULL
        END
    ) AS c1,
    MAX(
        CASE WHEN key = 'c2' THEN
            value
        ELSE
            NULL
        END
    ) AS c2,
    MAX(
        CASE WHEN key = 'c3' THEN
            value
        ELSE
            NULL
        END
    ) AS c3
FROM
    vtable
GROUP BY
    uid

PrestoのMap型について

Prestoには Mapというデータ型があります。
これは Pythonのdictと似たようなデータ型で、 キーと値のペアで構成されるものです。

ドキュメントを見ると、
MAP(ARRAY['foo', 'bar'], ARRAY[1, 2]) というMAPを作る関数の例が紹介されていますが、実行すると次のような結果を得ることができます。


-- クエリ
SELECT
    MAP(ARRAY['foo', 'bar'], ARRAY[1, 2])

-- 以下結果
{"bar":2,"foo":1}

以前紹介した、 TD_PARSE_AGENT という関数がありますが、この結果もMapです。
参考: TreasureDataのTD_PARSE_AGENT関数が便利

Map から、キーをしていて値を取得するときは、上のTD_PARSE_AGENTの記事でやっているように、map に [‘key’]をつけて取得するか、
ELEMENT_AT(map(K, V), key) を使います。
個人的にはPythonなどと近い書き方の [ ] を使う方法の方が好きです。

そのほかの、Mapの使い方ですが、ドキュメントの
6.17. Map Functions and Operators
というページにまとまっているのでこちらがわかりやすいです。

さて、 Tableのある列の値を キーとし、別の列の値を値とするMapを作りたくなる場合があります。
そのようなときは、 MAP_AGG という関数が使えます。
ドキュメントでは別のページにあるので探しにくいのですが、こちらにあります。
6.14. Aggregate Functions
の Map Aggregate Functions。

キーに指定した列に重複があったらどうなるかが気になったのですが、
試してみたところ、バリューの列のどれかの値が何らかの規則で一つだけ選ばれて採用されるようです。
(どのような基準で選ばれているのかは結局わかりませんでした。使うときは注意が必要ですね。)

キーに重複があり、バリューを(配列の形で)全部残したいときは MULTIMAP_AGG が使えます。

MAP_AGG で MAPを作って、 map[key] で値にアクセスする、ということだけ覚えておけば、
特に問題なく使うことができると思います。