多変量の時系列データの良いところは、データから因果関係わかることです。その一つが今回テーマとするインパルス応答関数です。ざっくりいうと、変数の変化が他の変数に与える影響を見ることで因果関係を確認できるようにします。
今回は下記の図のようにインパルス応答関数を行い、それをグラフ化します。このグラフを見ると、各変数が未来に関してどれだけ影響を与えているかがわかりますね。まさしく、見聞色の覇気を鍛えすぎた状態を体感できるのです。
兎にも角にも実践してみることをモットーとし、まず先にPythonを使ってインパルス応答関数を行い、その後に理論的背景を見ていきます。今回のコードはこちらに置いてあります。
インパルス応答関数の実装
多変量時系列データの因果関係がわかることは素晴らしいですよね。インパルス応答関数を使えばそれができます。どうやっているかを簡単に言うと、誤差（撹乱項）の１単位の変化による先々の変数の変化量を表したものです。言葉で見ると、簡単そうに見えますね。実は、少ないコードで実行できるというお手頃関数なのです。詳しいことは後ほど述べますので、とにかく実装していきましょう。
データの準備
それでは早速実践していきますが、一にも二にもデータが無いと始まりません。そこで、多変量時系列データを取得しましょう。とはいっても秒でGetできます。
データの取得
今回使用するデータは、前回の予測誤差分散分解の実装で使ったものと同じデータである、アメリカの経済データです。これは、statsmodelsモジュールに組み込まれているため、簡単に獲得できます。
import pandas as pd import statsmodels.api as sm from statsmodels.tsa.api import VAR macro_data = sm.datasets.macrodata.load_pandas().data
データの加工
モデルで扱えるデータにするため取得したデータのデータフレームのインデックスに、時期を表すクオーターを代入し、今回扱う変数である「realgdp」「realcons」「realinv」だけを取り出します。経済データの因果関係がわかれば、財布の紐が緩んでしまうかもしれませんね。
- realgdb：実質国内総生産
- realcons：実質個人消費支出
- realinv：実質国内総投資
macro_data['year'] = macro_data['year'].astype(int) macro_data.index = pd.date_range( str(macro_data['year'].min()), periods=len(macro_data), freq='Q' ) columns = ['realgdp', 'realcons', 'realinv'] target_data = macro_data[columns].copy() target_data = np.log(target_data) - np.log(target_data.shift(1)) target_data.dropna(inplace=True) target_data *= 100
VARのモデル構築の記事で確認しましたが、今回使用するデータは非定常的です。つまり、だんだん増加したり（減少したり）といったトレンドを持っています。そこで、対数差分を行い定常的なデータに変換しています。時系列データにおける対数差分は、前期からの比率を表すので、このあとの結果がわかりやすくなるように100をかけています。
VARモデルの構築
データの加工が終わったのでVARモデルを構築します。たった数行でできます。
model = VAR(target_data) model_result = model.fit(maxlags=15, ic='aic') model_result.summary()
AIC基準でラグ３のVARモデルを選択しています。VARモデルができたので、このモデルを使いインパルス応答関数を実装します。
インパルス応答関数の実装
インパルス応答関数は、statsmodelsモジュールを使えば簡単に実行できます。下記では、10期先までの結果を表示しています。インパルス応答関数は、非直交化と直交化の２種類あり、plotのorth引数で切り替えることができます。orthをTrueにすると直交化された結果が表示されます。
irf = model_result.irf(10) irf.plot(orth=True)
上の図を見ると、各変数の因果関係がわかりますね。早速、結果を分析していきますが、インパルス応答関数の良いところは、１ショックに対して変化量の値もわかることです。ただし、statsmodelsのirfの直交化バージョンでは、１ショックを１標準偏差とした結果が出てきます。このときの、１標準偏差は、realgdp=1.0、realcons=1.24、realinv=1.92です。なぜこの値になるかは、後ほど、少し下の内容を見てください。また、VARに入れたデータは100を掛けてあるので、表示された結果の値は％表示です。
- 国内総生産の1%の増加において、直近の変化が一番大きいのは国内総投資であり、同時期には3%増加。ただし、１クオーター先までは影響あるが、その後の影響は少ない
- 個人消費支出の1.24%の増加において、国内総生産は１クオーター先で0.4%増加するが、その後は影響がなくなっていく。国内総投資に対しては、同時期に-1.56%の影響だが、１クオーター先では1.96%の増加
- 国内総投資の1.92%の増加において、他の変数への影響は小さい
ざっくりまとめると、個人消費が増えたときに１クオーター遅れて、国内総生産と投資が増加するようですね。また、国内総生産が増加すると、１クオーター先の個人消費が上がるという結果ですが、個人消費が国内総生産を支えているという関係性から、国内総生産が上がっている状態だと、１クオーター後も個人消費が増加する状態が続くと解釈することが良さそうですね。
直交化インパルス応答関数における分解した撹乱項の標準偏差
先程、直交化インパルス応答関数において、１標準偏差の値が出てきたので、その値の求め方を見ていきます。
sigma = model_result.sigma_u L = np.linalg.cholesky(sigma) L_inv = np.linalg.inv(L) L_t_inv = np.linalg.inv(L.T) D = L_inv * sigma * L_t_inv new_sigma_sqrt = np.sqrt(np.diag(D)) new_sigma_sqrt
なぜ、この計算式で直交化インパルス応答関数の撹乱項の標準偏差が求まるかは、以下の章を見てもらえばわかるかと思います。上記の結果から、realgdpは1.0、realconsは1.24、realinvは1.92が１標準偏差であることがわかります。
インパルス応答関数の理論
インパルス応答関数を使えば、多変量時系列のデータにおいて各変数の因果関係がわかります。では、なぜそれができるのかに関して見ていきます。大雑把に言うと、インパルス応答関数とは、「ある変数の変化が他の変数にどれだけ影響を与えるかを定量化」したものです。
変数 \(x\) から変数 \( y \) に対しての因果関係があるならば、\( x \) の値が変われば \( y \) の値も必ず変わることになります。なので、\( x \) の変化に対して \( y \) が変化するかを見ることで因果関係がわかることになります。時系列界隈では、変数を変化させることをショックを与えると呼んでいるようです。
変数にショックを与えるやり方として、非直交化インパルス応答関数と直交化インパルス応答関数の２通りのやり方があります。まずは、非直交化インパルス応答関数について見ていきます。
非直交化インパルス応答関数
インパルス応答関数を行うために最初にやることは、VARモデルを構築し、各変数を時系列モデルで定式化することです。VARモデルは以下の式で表されます。
- \( n \) : 変数の数
- \(p\) : ラグ
- \( \mathbf{y}_t \) : 時刻 \(t\) のときの \(n\) 個の時系列データのベクトル
- \( \mathbf{\phi_t} \) : 時刻 \( t \) の変数に対する \( n \times n \) の係数行列
- \( \mathbf{c} \) : \(n \times 1 \) 定数ベクトル
- \( \mathbf{\epsilon}_t \) : 時刻 \( t \) のときの \( n \) 個の誤差ベクトル
- \( \mathbf{\Sigma} \) : \( n \times n \) の共分散行列
各変数を定式化できたので、ショックを与えたときの各変数の変化を見ることで因果関係がわかります。では、どこにショックを与えるかと言うと、それは撹乱項の \( \epsilon_{i,t}\) です。 \( \epsilon_{i,t} \) を１単位（あるいは１標準偏差）変化させたときの \( y_{j,t+k} \) の変化を見ます。ここで、\( k \) が出てきているように、直後の \( y_{i,t+1} \) の変化だけではなく、\( k \) 期先の変化を確認できます。
\( k \) 期先の \( y_i \) のショックに対する \( y_j \) の変化を偏微分で表せば、以下の形になります。
上記の式を \( k = 0 \) から逐次的に解いていくことで、指定した期間先までのインパルス応答関数が求まります。定式化だけを見ると意外とあっさりしていることがわかります。それはVARモデル自体が線形的でシンプルな構造をしており、応用がしやすいからですかね。つまり、シンプルに考えることが大事なのですね。
ここで具体例を書くのが筋なのかもしれませんが、そこまですると紙面を大幅に割いてしまうので、省略させていただき、この方法の問題点を取り上げていきます。
問題点は何かというと、撹乱項の \( \epsilon_{i,t}\) と \( \epsilon_{j,t} \) が相関していたら、撹乱項の \( \epsilon_{i,t} \) にショックを与えたときに \( \epsilon_{j,t} \) も変化し、\( y_{j,t+k} \) が \( y_{i,t} \) 以外の影響を受けてしまうことです。
この問題を解決するために、直交化を行い無相関な撹乱項に分解します。そのため、次に見ていくインパルス応答関数は直交化インパルス応答関数と呼ばれます。これに対比する形で、上記の方法は非直交化インパルス応答関数と呼ばれています。
直交化インパルス応答関数
非直交化インパルス応答関数の問題点は、撹乱項同士の相関を無視していることでした。そこで、互いに無相関な撹乱項に分解し、その無相関な撹乱項にショックを与えることで、インパルス応答関数を行います。まず確認すべきことは、撹乱項 \( \epsilon_t \) の 共分散行列 \( \mathbf{\Sigma} \) です。この行列の対角要素以外の要素により、撹乱項同士が相関しているかがわかります。つまり、この行列が対角行列ならば、撹乱項同士は無相関です。
何をすればよいかと言うと、撹乱項 \( \mathbf{ \epsilon_t }\) を分解して作った新たな撹乱項の共分散行列が、対角行列になれば良いということです。対角行列になれば、New撹乱項同士は無相関になり、インパルス応答関数の問題は解決です。
では、どうやって対角行列の新しい撹乱項を作るかですが、行列演算を使い実は難なく？できます。順番にやり方を見ていきます。まず注目するのは、撹乱項 \( \epsilon_t \) の共分散行列が正定値行列であることです。いきなり斜め上からそれっぽい行列用語が出できて、即閉じしかけたかもしれませんが、ここでは正定値行列について深追いしません。正定値行列であるため、誤差分散行列は以下の形に分解できます。
- \( \mathbf {L} \) := 下三角行列
- \( \mathbf {D} \) := 対角行列
上記の正定値行列の分割の仕方を修正コレスキー分解と呼びます。なんと、対角行列がさくっと現れました。共分散行列を対角行列にしたいという明確な目標があるので、右辺を対角行列 \( \mathbf {D} \) だけの形になるように式変形してみましょう。
\[
\begin {eqnarray}
\mathbf { \Sigma } = \mathbf {L} \mathbf {D} \mathbf {L}^T \\
\mathbf {L}^{-1} \mathbf { \Sigma } (\mathbf {L}^T)^{-1} = \mathbf {D}
\end {eqnarray}
\]
共分散行列 \( \mathbf { \Sigma } \) を下三角行列の逆行列 \( \mathbf{ L }^{-1} \) で挟むと対角行列になることがわかります。つまり、New撹乱項の相関はなくなります。
このときに出てきた下三角行列の逆行列 \( \mathbf{ L }^{-1} \) が重要な役割を持っています。下三角行列の逆行列に撹乱項 \( \mathbf {\epsilon}_t \) をかけたものを \( \mathbf {u}_t \) とします。
\[
\mathbf {u}_t = \mathbf {L} ^ {-1} \mathbf { \epsilon }_t
\]
なぜ上の式が急に出てきたかと言うと、それは答えを知っているからではありますが、上の式の分散を求めてみるとしっくりくるでしょう。
\[
\begin {eqnarray}
\rm {Var} ( \mathbf {u}_t ) & = & \rm {Var} ( \mathbf {L} ^ {-1} \mathbf { \epsilon }_t ) \\
& = & \mathbf {L}^{-1} \rm {Var} ( \mathbf {\epsilon}_t ) (\mathbf {L}^{-1})^T \\
& = & \mathbf {L}^{-1} \rm {Var} ( \mathbf {\epsilon}_t ) (\mathbf {L}^t)^{-1} \\
& = & \mathbf {D}
\end {eqnarray}
\]
なんと、\( \mathbf {u}_t \) の分散は対角行列である \( \mathbf { D } \) になりました。ということは、New撹乱項である \( \mathbf {u}_t \) の各要素は相関していません。これがまさしく求めていた、新しい撹乱項です。\( u_{i,t} \) にショックを与えて \( y_{ j,t+k } \) の変化を見ることで、撹乱項の相関問題を解決したインパルス応答関数を求めることができます。
直交化インパルス応答関数における、\( k \) 期先の \( y_i \) のショックに対する \( y_j \) の変化を偏微分を使って表すと以下の形になります。
これで直交化インパルス応答関数の計算方法がわかりましたね。話が込み入ったので、\( \mathbf {u}_t \) は何だったかを思い出すと、撹乱項 \( \mathbf {\epsilon}_t \) の誤差分散行列 \( \mathbf {\Sigma} \) を修正コレスキー分解したときに出てくる下三角行列 \( \mathbf {L} \) の逆行列に撹乱項 \( \mathbf {\epsilon}_t \) を掛けたものでした。つまり、計算できますね。
行列の分散
上記の一連の計算式において、頭を一番悩ませたのは行列の分散でした。まぁ、知らなかっただけすが。備忘録ということで行列の分散の式展開を書いておきます。ここで、\( \mathbf{A} \) は行列であり、\( \mathbf{ x^T }\) は列ベクトルです。
\[
\begin {eqnarray}
\rm {Var} ( \mathbf {A} \mathbf {x}^T ) & = & \frac {(\mathbf {A} \mathbf {x}^T - \overline { \mathbf {A} \mathbf {x}^T}) (\mathbf {A} \mathbf {x}^T - \overline { \mathbf {A} \mathbf {x}^T})^T } {n-1} \\
& = & \frac {(\mathbf {A} \mathbf {x}^T - \mathbf {A} \overline { \mathbf {x}^T}) (\mathbf {A} \mathbf {x}^T - \mathbf {A} \overline { \mathbf {x}^T})^T } {n-1} \\
& = & \frac {(\mathbf {A} \mathbf {x}^T - \mathbf {A} \overline { \mathbf {x}^T}) (\mathbf {x}^T \mathbf {A}^T - \overline { \mathbf {x }^T }\mathbf {A}^T) } {n-1} \\
& = & \frac {\mathbf {A} ( \mathbf {x}^T - \overline { \mathbf {x}^T}) (\mathbf {x}^T - \overline { \mathbf {x}^T }) \mathbf {A}^T } {n-1} \\
& = & \mathbf {A} \frac {( \mathbf {x}^T - \overline { \mathbf {x}^T}) (\mathbf {x}^T - \overline { \mathbf {x}^T } )} {n-1} \mathbf {A}^T \\
& = & \mathbf {A} \rm {Var} (\mathbf {x}^T) \mathbf {A}^T \\
\end {eqnarray}
\]
終わりに
今回は多変量時系列データにおける因果関係を求める方法に関して書きました。途中の説明で力尽きて書ききれなかったこともあるので、もっと具体的に知りたいよって方はこの本をぜひ御覧ください。この本は時系列分析全般の学習において必ずおすすめされますね。
