ゆるくJuliaで直交化インパルス応答関数の実装

16234 ワード
仮想通貨 株価 Julia 時系列分析 計量時系列分析 Julia テキストリンク

やりたいこと

VARモデルを使った時系列分析のツールは3つありました。グレンジャー因果性検定、インパルス応答関数、分散分解ですね。前回はそのうちのひとつ、グレンジャー因果性検定を実装しました。今回はインパルス応答関数を実装したいと思います。実装にあたって過去記事のコードを利用します。別ファイルのコードを呼び出す方法については、Juliaで自作のモジュールや関数をファイル単位で分割する方法をご覧ください。

ふわっと、インパルス応答関数

あなたはビットコインの価格を予想したいと考えています。いろいろ考えてみた結果、仮想通貨の値動きは互いに影響し合っている、という仮説を立てました。

そこで、いくつかの仮想通貨をピックアップしてそれらを変数とするVARモデル作ることに決めました。選んだ仮想通貨は、「ビットコイン」「イーサリアム」「ドージコイン」です。

あなたは選んだ３つの仮想通貨のティックデータをモデルに学習させ、その学習させたモデルを使ってインパルス応答関数を計算しました。

計算した結果、イーサリアムの価格がいくらか上昇したときに、翌営業日のビットコインの価格が5%上昇しその影響は2日間で消えることがわかりました。

あなたはイーサリアムの値動きに一喜一憂するのでした。めでたしめでたし。

カチッと、インパルス応答関数

インパルス応答関数のイメージは掴めたでしょうか？少し硬くいうと、インパルス応答関数とは、ある変数にショックを与えたときに他の変数にどれほどの影響を与えるかを計算します。

インパルス応答関数の定義に入る前に、少し確認しておきたいことがあります。インパルス応答関数は実は二つあって、非直交化インパルス応答関数と直交化インパルス応答関数があります。一般にインパルス応答関数といえば直交化インパルス応答関数を指すことのほうが多いみたいです。なので、本記事でインパルス応答関数というときは特に断りがない限り、直交化インパルス応答関数を指すものとします。

ではインパルス応答関数の定義をします。

直交化インパルス応答関数 沖本本p87

攪乱項の分散共分散行列の三角分解を用いて,攪乱項を互いに無相関な攪乱項に分解したとき,お互いに無相関な攪乱項は直交化攪乱項といわれる.また, $y_{i}$ の直交化攪乱項に1単位または1標準偏差のショックを与えたときの $y_{i}$ の直交化インパルス応答を時間の関数としてみたものは, $y_{j}$ に対する $y_{i}$ の直交化インパルス応答関数といわれる.

では、簡単にですが次の章でインパルス応答関数を求めてみたいと思います。

★非直交化インパルス応答と直交化インパルス応答関数の違い★ VARモデルの攪乱項は異なる時点間では相関は0でなければいけませんが、同時点では相関を持っていいです。つまり攪乱項の共分散行列が対角行列である必要はありません。そうすると、攪乱項にショックを与えていくらか変化させたとき、相関のある攪乱項も一緒に変化するはずです。 非直交化インパルス応答関数は同時点での攪乱項の相関を無視して、ショックの影響を計算してしまいます。これを修正したものが直交化インパルス応答関数です。

インパルス応答関数を求める

まずはVARモデルを定義します。簡単のために2変数VAR(1)を考えます。

$
\boldsymbol{y_{t}} = \boldsymbol{c} + \boldsymbol{\Phi_{1}} \boldsymbol{y_{t-1}} + \boldsymbol{\epsilon_{t}} \quad \boldsymbol{\epsilon_{t}} ～ W.N.( \Sigma ) \tag{1}
$

$c$は2×1のベクトル、$\Phi$は2×2の行列、$\epsilon$は2×1のベクトルです。攪乱項$\epsilon$を直交化攪乱項に変換します。共分散行列にコレスキー分解を用いると、

$$
\Sigma = PP'
$$

となります。この$P$と直交化攪乱項 $\boldsymbol{u_{t}}$ を用いると攪乱項は次のように表せます。

$$
\boldsymbol{\epsilon_{t}} = P \boldsymbol{u_{t}}
$$

式(1)に代入すると

$
\boldsymbol{y_{t}} = \boldsymbol{c} + \boldsymbol{\Phi_{1}} \boldsymbol{y_{t-1}} + P \boldsymbol{u_{t}} \quad \boldsymbol{u_{t}} ～ ~ W.N.( D ) \tag{2}
$

となります。$D$は対角行列です。

インパルス応答関数は式(2)を直交化攪乱項で微分することで求めることができます。
※IRF -> Impulse Response Function

$$
IRF(k) =
\frac{\partial \boldsymbol{y_{t+k}}}{\partial \boldsymbol{u_{t}}} \quad k=0,1,2, ... \tag{3}
$$

実際に微分してみましょう。まずは$IRF(0)$を求めます。$P$は下三角行列で成分は次の通りです。

\begin{align}
P &= 
\begin{bmatrix}
P_{11} & 0 \\
P_{21} & P_{22}
\end{bmatrix} \\
\quad \\
IRF(0) &= 
\frac{\partial \boldsymbol{y_{t}}}{\partial \boldsymbol{u_{t}}}  \\
&= \begin{bmatrix}
\frac{\partial \boldsymbol{y_{1,t}}}{\partial \boldsymbol{u_{1,t}}} & \frac{\partial \boldsymbol{y_{1,t}}}{\partial \boldsymbol{u_{2,t}}} \\
\frac{\partial \boldsymbol{y_{2,t}}}{\partial \boldsymbol{u_{1,t}}} & \frac{\partial \boldsymbol{y_{2,t}}}{\partial \boldsymbol{u_{2,t}}}
\end{bmatrix} \\
&= \begin{bmatrix}
P_{11} & 0 \\
P_{21} & P_{22}
\end{bmatrix} \\
&= \begin{bmatrix}
IRF(0)_{11} & IRF(0)_{12} \\
IRF(0)_{21} & IRF(0)_{22}
\end{bmatrix} \\
\end{align}

添え字がたくさん出てきますが、そこまで複雑なことはしていません。どこかの成分を実際に微分すればわかりやすいと思います。
ここで、$IRF_{12}$は$y_{2}$の攪乱項にショックを与えたときの、$y_{1}$の反応です。

次にIRF(1)を求めます。

\begin{align}
IRF(1) &= 
\frac{\partial \boldsymbol{y_{t+1}}}{\partial \boldsymbol{u_{t}}}  \\
&= \begin{bmatrix}
\frac{\partial \boldsymbol{y_{1,t+1}}}{\partial \boldsymbol{u_{1,t}}} & \frac{\partial \boldsymbol{y_{1,t+1}}}{\partial \boldsymbol{u_{2,t}}} \\
\frac{\partial \boldsymbol{y_{2,t+1}}}{\partial \boldsymbol{u_{1,t}}} & \frac{\partial \boldsymbol{y_{2,t+1}}}{\partial \boldsymbol{u_{2,t}}}
\end{bmatrix} \\
&= \begin{bmatrix}
\Phi_{11} \frac{\partial \boldsymbol{y_{1,t}}}{\partial \boldsymbol{u_{1,t}}} 
+ \Phi_{12} \frac{\partial \boldsymbol{y_{2,t}}}{\partial \boldsymbol{u_{1,t}}}
& \Phi_{11} \frac{\partial \boldsymbol{y_{1,t}}}{\partial \boldsymbol{u_{2,t}}} 
+ \Phi_{12} \frac{\partial \boldsymbol{y_{2,t}}}{\partial \boldsymbol{u_{2,t}}} \\
\Phi_{21} \frac{\partial \boldsymbol{y_{1,t}}}{\partial \boldsymbol{u_{1,t}}} 
+ \Phi_{22} \frac{\partial \boldsymbol{y_{2,t}}}{\partial \boldsymbol{u_{1,t}}} 
& \Phi_{21} \frac{\partial \boldsymbol{y_{1,t}}}{\partial \boldsymbol{u_{2,t}}} 
+ \Phi_{22} \frac{\partial \boldsymbol{y_{2,t}}}{\partial \boldsymbol{u_{2,t}}} 
\end{bmatrix} \\
&= 
\begin{bmatrix}
\Phi_{11} IRF(0)_{11} + \Phi_{12} IRF(0)_{21} & \Phi_{11} IRF(0)_{12} + 
\Phi_{12} IRF(0)_{22} \\
\Phi_{21} IRF(0)_{11} + \Phi_{22} IRF(0)_{21} & \Phi_{21} IRF(0)_{12} + 
\Phi_{22} IRF(0)_{22}
\end{bmatrix} \\
&= 
\begin{bmatrix}
\Phi_{11} & \Phi_{12} \\
\Phi_{21} & \Phi_{22}
\end{bmatrix}
\begin{bmatrix}
IRF(0)_{11} & IRF(0)_{12} \\
IRF(0)_{21} & IRF(0)_{22}
\end{bmatrix} \\
&= 
\boldsymbol{\Phi_{1}}
\begin{bmatrix}
IRF(0)_{11} & IRF(0)_{12} \\
IRF(0)_{21} & IRF(0)_{22}
\end{bmatrix} 
\end{align}

どうやら $IRF(k)$ はひとつ前の $IRF(k-1)$ と係数の積で表すことができそうです。

IRF(k) = \boldsymbol{\Phi_{1}} IRF(k - 1) \tag{4}

式(4)を実装することを目標にコードを書きます。直接の実装はcalcirf( n ) = n == 0 ? C.L : W * calcirf( n - 1 )でしています。

ソースコード

orth_irf.jl

# 過去に作ったグレンジャー因果性検定のコードを呼び出す
include("granger_causality.jl")

# 直交化インパルス応答関数
function irf( m::var, periods=1 )  
    # 定数は微分すると消える　=> 沖本本 p89
    W = m.params[2:end, :]'
    D = size(W)[1]
    # コレスキー分解
    C = cholesky( m.Σ )

    # IRFの逐次計算を行う関数
    calcirf( n ) =  n == 0 ? C.L : W * calcirf( n - 1 )

    # IRFを再帰的に求める
    res = zeros( Float64, D, D, periods + 1 )
    for n in 0:periods
        res[:,:,1+n] += calcirf( n )
    end

    # periods * D * D にreshape
    res_sorted = zeros( Float64, periods + 1, D, D )
    s = 1
    for i in 1:D
        for j in 1:D
            res_sorted[:, j, i] += res[s:D*D:end]
            s += 1
        end
    end
    res_sorted
end

# 直交化インパルス応答関数のプロット
function plot_irf( irfs )
    periods = size( irfs )[1] - 1 
    D = size( irfs )[2]

　　p = []
    for i in 1:D
        for j in 1:D
            tmp = plot( 0:periods, irfs[:,j,i], xticks=(0:periods), label="irf"*"$j"*"$i" )
            push!( p, tmp )
        end
    end 
    plot( p..., layout=(D, D) )
end

使い方

サンプルを用意しました。Xには適当な時系列データを入れてください。画像はplot_irfの出力例です。

usecase.jl
# size( X ) = ( 100, 3 )

# 3変数var(2)
params = var( 2, rand(2, 3), rand( 2, 2 ) )
# 学習
fit!( params, X )

irfs = irf( params, 10 )
plot_irf( irfs )

【重要】再帰的構造

なんで最後に書いたかというと、インパルス応答関数の出力例の画像を見てほしかったからです。

まずは一番左の列のx軸が0の部分を見てください。この列は1つ目の変数にショックを与えたときの各変数への影響をプロットしたものです。$IRF_{12}, IRF_{13}$ は0になっていますね？これは直交化攪乱項を計算するときに下三角行列を用いたことによって起きるものです。式(2)を実際に書き下すとわかるのですが、$y_{1,t}$ に含まれる攪乱項は$u_{1, t}$ だけなのに対して、$y_{2,t}$ に含まれる攪乱項は｛ $u_{1, t}, u_{2, t}$ ｝の二つです。

これが何を意味するかというと、$y_{1}$ の攪乱項は $y_{2}$ に影響を与えるのに、$y_{2}$ 特有の攪乱項は $u_{2}$は $y_{1}$ に影響を与えないのです。こういう構造のことを再帰構造といいます。

仮想通貨の例をまた持ち出すと、変数を、ビットコイン→イーサリアム→ドージコインと並べるとき、同時点ではイーサリアムのショックはビットコインに影響を与えることはできないという意味です。

インパルス応答関数はこのような仮定を置いているので、変数の並び替えには気を付けなければいけません。

参考文献