分析ノート

思いのほか理論編が長くなってしまいすでにインパルス応答で４記事目ですが、ようやく実装の話です。
グレンジャー因果性と同様に、statsmodelsで学習したVARモデルはインパルス応答関数のメソッドも持っています。

メソッド名は irf です。
statsmodels.tsa.vector_ar.var_model.VARResults.irf

早速やってみましょう。
モデルは、直交化インパルス応答関数の計算例の記事と同じです。
サンプルデータの生成と、パラメーターの推定まで完了しているとします。
今回は次の状態になりました。

# 学習したモデルのパラメーター
result = model.fit(maxlags=1)
print(result.summary())
"""
  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Sun, 23, Feb, 2020
Time:                     21:45:36
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                   0.924001
Nobs:                     99.0000    HQIC:                  0.830357
Log likelihood:          -312.903    FPE:                    2.15278
AIC:                     0.766721    Det(Omega_mle):         2.02800
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.177695         0.918241            0.194           0.847
L1.y1         0.553623         0.112916            4.903           0.000
L1.y2         0.169688         0.169945            0.998           0.318
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.306800         0.408504            3.199           0.001
L1.y1         0.085109         0.050234            1.694           0.090
L1.y2         0.760450         0.075605           10.058           0.000
========================================================================

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

サマリーの最後に、撹乱項の相関行列がありますが、分散共分散行列がみたいのでそちらも表示しておきましょう。
sigma_uに入ってます。

# 分散共分散行列
print(result.sigma_u)
"""
[[4.12100608 1.1546159 ]
 [1.1546159  0.81561201]]
"""

さて、まずは非直行化インパルス応答関数を計算して可視化してみましょう。
ちなみに、直行化/非直行化を指定するのは可視化(plot)の引数であり、計算時点で両方とも算出されています。

import matplotlib.pyplot as plt

# 10期先までの非直交化インパルス応答関数のプロット
period = 10
irf = result.irf(period)
irf.plot()
plt.show()

結果がこちらです。

想定通りの結果になっていますね。

続いて、(直行化)インパルス応答関数です。
用語としては、単にインパルス応答関数と言ったら、直行化の方をさすのですが、ライブラリでは明示的に、
orth=Trueを指定しなければなりません。(デフォルトはFalse)

計算は完了しているので、可視化のみです。

irf.plot(orth=True)
plt.show()

結果がこちら。

0期先の値が特徴的で、 y2からy1へは影響していないのに、y1からy2へは影響が発生していますね。
また、もう一点注意しないといけない点があります。
非直行化の方は、y1からy1 および y2からy2への0期先の影響の値が1であることから分かる通り、
1単位のショックを与えた時のインパルス応答関数を計算されていました。
しかし、直行化の方は、そうではなく、1標準偏差のショックを与えた時のインパルス応答関数の値になっています。
(この辺の挙動はドキュメントのどこに書いてあるのだろう？)

結果をプロットするのではなく、値を利用したい場合もあると思います。
その場合、次の2変数にプロットされた値がそれぞれ入っていますのでこれが使えます。

irf.irfs
irf.orth_irfs

今回は直交化インパルス応答関数を実際に計算してみます。

例として扱うのは、非直交化インパルス応答関数の記事と同じ次の例です。
$$
\left\{\begin{matrix}\
y_{1t}&=& -1+ 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{1t}\\\
y_{2t}&=& 1 + 0.1 y_{1,t-1}+ 0.8 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
,\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\boldsymbol{\Sigma})
$$
$$
\boldsymbol{\Sigma} = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

さて、まずは$\boldsymbol{\Sigma}$を分解する必要があります。

$$
\begin{align}
\mathbf{A}=
\left(\begin{matrix}1&0\\0.3&1\end{matrix}\right),\\
\mathbf{D}= Var(\mathbf{u}_t)
\left(\begin{matrix}4&0\\0&0.64\end{matrix}\right)
\end{align}
$$
とすると、
$$
\boldsymbol{\Sigma} = \mathbf{ADA^{\top}}
$$
が成立します。
このため、$\boldsymbol{\varepsilon_t}$は次のように、
直交化撹乱項に分解できます。
$$
\left\{\begin{align}\
\varepsilon_{1t}&=u_{1t}\\\
\varepsilon_{2t}&=0.3u_{1t}+u_{2t}\
\end{align}\right.
$$
この撹乱項を用いることで、元のモデルがこのように変わります。

$$
\left\{\begin{matrix}\
y_{1t}&=& -1&+ &0.6 y_{1,t-1}&+ &0.3 y_{2,t-1}&+&u_{1t}\\\
y_{2t}&=& 1 &+ &0.1 y_{1,t-1}&+ &0.8 y_{2,t-1}&+&0.3u_{1t}+u_{2t}\
\end{matrix}\
\right.\
$$
$$
\mathbf{u}_t \sim W.N.(\mathbf{D})
$$
$$
\mathbf{D} = \left(\begin{matrix}4&0\\0&0.64\end{matrix}\right)
$$

この形に変形できれば、非直交化インパルス応答関数と同じ様に計算が出来ます。
例えば、$y_1$に1単位のショックを与えたときの同時点におけるインパルス応答は次になります。
$$
\begin{align}
IRF_{11}(0) = \frac{\partial y_{1t}}{\partial u_{1t}} = 1\\
IRF_{21}(0) = \frac{\partial y_{2t}}{\partial u_{1t}} = 0.3
\end{align}
$$
$IRF_{21}(0)\neq0$となるのが大きな特徴です。

ただし、逆に$y_2$に1単位のショックを与えても、$y_1$には影響がありません。
要するに次のようになります。
$$
\begin{align}
IRF_{12}(0) = \frac{\partial y_{1t}}{\partial u_{2t}} = 0\\
IRF_{22}(0) = \frac{\partial y_{2t}}{\partial u_{2t}} = 1
\end{align}
$$
これが三角分解の仮定により発生している現象であり、変数の並べ方が結果に影響してしまっている点です。
このためインパルス応答関数を使う時は、変数の並べ方に気を使う必要があります。

1期後先以降の値は、非直交化インパルス応答関数のときと同じ漸化式で逐次的に求めることが可能です。

前回の記事: 非直交化インパルス応答関数
の続きです。今回も引用元は沖本本です。

最後に書いたとおり、非直交化インパルス応答関数には、撹乱項の間の相関を無視しているという問題がありました。
これを解決するには、撹乱項を互いに無相関な撹乱項に分解して、それぞれにショックを与えたときの各変数の変化を調べることです。
しかし、一般に撹乱項を互いに無相関な撹乱項に分解することは不可能で、何かしらの仮定が必要になります。
そこで、撹乱項の分散共分散行列$\boldsymbol{\Sigma}$の三角分解を用いて、撹乱項を互いに無相関な撹乱項に分解できると仮定します。

定義 (直交化インパルス応答関数):
撹乱項の分散共分散行列の三角分解を用いて、撹乱項を互いに無相関な撹乱項に分解したとき、
互いに無相関な撹乱項は直交化撹乱項といわれる。また、$y_j$の直交化撹乱項に1単位または1標準偏差のショックを与えた時の、
$y_i$の直交化インパルス応答を時間の関数としてみたものは、$y_j$に対する$y_i$の直交化インパルス応答関数といわれる。

一般的に、インパルス応答(関数)というというと直交化インパルス応答(関数)を指すそうです。
(後日の記事でstatsmodelsによる実装紹介しますが、このモジュールでは非直交化のほうをデフォルトにしてるように見えます。
もしそうだとしたら気をつけて使わないといけないので、記事にする前にもう一回検証します。)

さて、もう少し具体的に定義します。
まず撹乱項 $\boldsymbol{\varepsilon}_t$ の分散共分散行列$\boldsymbol{\Sigma}$は正定値行列であるので、
$\mathbf{A}$を対角成分が1に等しい下三角行列、$\mathbf{D}$を対角行列として、
$$
\boldsymbol{\Sigma} = \mathbf{ADA^{\top}}
$$
と三角分解することが出来ます。
(と、本には書いてありますが、分散共分散行列が正定値なだけでなく、対称行列であることも考慮する必要あると思います。)

このとき、直交化撹乱項 $\mathbf{u}_t$は、
$$
\mathbf{u}_t = \mathbf{A}^{-1}\boldsymbol{\varepsilon}_t
$$
で定義され、各変数固有の変動を表すとみなされます。
なお、
$$
\begin{align}
Var(\mathbf{u}_t) &= Var(\mathbf{A}^{-1}\boldsymbol{\varepsilon}_t) = \mathbf{A}^{-1} Var(\boldsymbol{\varepsilon}_t)(\mathbf{A}^{-1})^{\top} \\
& =\mathbf{A}^{-1}\boldsymbol{\Sigma}(\mathbf{A}^{\top})^{-1}=\mathbf{D}
\end{align}
$$
であることから、$\mathbf{u}_t$が互いに無相関であることが分かります。
また、それぞれの分散も確認できます。

$y_j$のショックに対する$y_i$の(直交化)インパルス応答関数は、$u_{jt}$に1単位のショックを与えたときの、
$y_i$の反応を時間の関数と見たものであり、
$$
IRF_{ij}(k) = \frac{\partial y_{i,t+k}}{\partial u_{jt}}, \ \ \ k=0, 1, 2, \dots.
$$
と表されます。
1単位ではなくて、1標準偏差のショックを与えたときの変数の反応を見たい時は、
$IRF_{ij}(k)$に$\sqrt{Var(u_{jt})}$を掛けることで求まります。(非直交化のときと同じですね)

1標準偏差のショックを与えたときのインパルス応答関数を求める方法として、
三角分解の変わりにコレスキー分解を用いて撹乱項を分解して、
その分解した撹乱項に1単位ショックを与える方法も紹介されていますが、
撹乱項の絶対値が標準偏差倍違うだけの話なので省略します。

グレンジャー因果性で複数の過程の間の関係を調べる方法を紹介してきましたが、グレンジャー因果性には
関係の強さを測れないという問題がありました。
そこで、それを補うツールとしてインパルス応答関数(IRF、 impulse response function)というものがあります。

例によって、沖本先生の本を参照しながら紹介させていただきます。

インパルス応答関数はある変数にショックを与える(要するに少し値を変動させる)と、
それが他の変数に与える影響を分析することができるものです。
VARモデルにおいては、ショックの識別の仕方によって、複数のインパルス応答関数が存在しますが、
まずは非直交化インパルス応答関数について説明します。

一般的な$n$変量のVAR(p)モデルを考えます。
$$
\mathbf{y}_t = \mathbf{c}+\boldsymbol{\Phi}_1\mathbf{y}_{t-1}+\cdots+\boldsymbol{\Phi}_p\mathbf{y}_{t-p}+\boldsymbol{\varepsilon}_t,
\ \ \ \boldsymbol{\varepsilon}_t\sim W.N.(\boldsymbol{\Sigma})
$$
本では、$\boldsymbol{\Sigma}$は対角行列ではないと仮定されていますが、対角行列の場合も同じだと思います。
(対角行列だと、その直後に登場する直行化インパルス応答関数と全く同じになってしまうで説明の都合上仮定されているのかと。)

この時、非直交化インパルス応答関数は次のように定義されます。

定義 (非直交化インパルス応答関数):
$y_{jt}$の撹乱項$\varepsilon$だけに、1単位(または、1標準偏差)のショックを与えた時の、
$y_{i,t+k}$の値の変化は、$y_j$のショックに対する$y_i$の$k$期後の非直行化インパルス応答と呼ばれる。
また、それを$k$の関数としてみたものは、$y_j$のショックに対する$y_i$の非直行化インパルス応答関数と言われる。

変化量なので、偏微分を用いて次のように表記できます。(1単位のショックを与えた場合。)
$$
IRF_{ij}(k) = \frac{\partial y_{i,t+k}}{\partial \varepsilon_{jt}}, \ \ \ k=0, 1, 2, \dots.
$$
1標準偏差のショックを与えた場合は、$IRF_{ij}(k)$に$\sqrt{Var(\varepsilon_{jt})}$を掛けることで求まります。

さて、具体的な計算方法ですが、$k=0$から逐次的に計算していくことで簡単に求まります。

次のモデルを例に、計算例が乗っているのでやってみましょう。
$$
\left\{\begin{matrix}\
y_{1t}&=& -1+ 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{1t}\\\
y_{2t}&=& 1 + 0.1 y_{1,t-1}+ 0.8 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
,\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\Sigma)
$$
$$
\Sigma = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

$y_1$に１単位のショックを与えた時の非直交化インパルス応答関数を具体的に計算します。
k=0の場合は簡単で、そのまま微分するだけです。
$$
\begin{align}
IRF_{11}(0) = \frac{\partial y_{1t}}{\partial \varepsilon_{1t}} = 1\\
IRF_{21}(0) = \frac{\partial y_{2t}}{\partial \varepsilon_{1t}} = 0
\end{align}
$$
この後は、多変数関数の合成関数の連鎖律を使って計算していきます。$k=1$の場合、$y_1$,$y_2$それぞれについて計算すると次のようになります。
$$
\begin{align}
IRF_{11}(1) &= \frac{\partial y_{1,t+1}}{\partial \varepsilon_{1t}}
= 0.6\times \frac{\partial y_{1t}}{\partial \varepsilon_{1t}} + 0.3 \times \frac{\partial y_{2t}}{\partial \varepsilon_{1t}}\\
&= 0.6 IRF_{11}(0) + 0.3 IRF_{21}(0) =0.6
\end{align}
$$
$$
\begin{align}
IRF_{21}(1) &= \frac{\partial y_{2,t+1}}{\partial \varepsilon_{1t}}
= 0.1\times \frac{\partial y_{1t}}{\partial \varepsilon_{1t}} + 0.8 \times \frac{\partial y_{2t}}{\partial \varepsilon_{1t}}\\
&= 0.1 IRF_{11}(0) + 0.8 IRF_{21}(0) =0.1
\end{align}
$$

この後の $k=1,2,3,…$についても次の漸化式で逐次的に計算できます。
$$
\left\{\
\begin{align}\
IRF_{11}(k) &= 0.6IRF_{11}(k-1) + 0.3IRF_{21}(k-1)\\\
IRF_{21}(k) &= 0.1IRF_{11}(k-1) + 0.8IRF_{21}(k-1)\
\end{align}\
\right.\
$$
この式の中に元のモデルの式の定数項ができませんが、それはインパルス応答関数が変化量に着目しているため、
無関係(微分で消えてしまった)からです。

さて、これで非直交化インパルス応答関数が計算できましたが、
これには一つ問題点があります。
というのも、最初のモデルの式の$\boldsymbol{\Sigma}$を見るとわかるとおり、$\varepsilon_{1t}$と$\varepsilon_{2t}$には
相関があることが仮定されています。
分散共分散行列から相関係数を計算すると
$$
Corr(\varepsilon_{1t}, \varepsilon_{2t}) = \frac{1.2}{\sqrt(4\times 1)} = 0.6
$$
となり、そこそこ強い相関です。
そのため、$\varepsilon_{1t}$にだけショックを与えて、$\varepsilon_{2t}$はそのままという、
非直交化インパルス応答関数には少し不自然な点が残ります。
すでに長くなってきたので、この点を解消する方法につては次の記事で紹介させていただきます。