月: 2023年2月

前回に続いてSVARモデル(構造VARモデル)の話です。前回はモデルの数式の形を紹介しただけだったので、今回は実際に仮想的なデータを作ってみてPythonのコードで推定してみます。

サンプルとして、以下のようなモデルを考えました。

$${\scriptsize
\left[\begin{matrix}1&0&0\\0.3&1&0\\-0.2&-0.3&1\end{matrix}\right]\mathbf{y_t}=
\left[\begin{matrix}0.8\\-0.5\\0.3\end{matrix}\right]+
\left[\begin{matrix}-0.1&0.3&0.4\\0.2&-0.2&-0.3\\-0.3&0.4&0.2\end{matrix}\right]\mathbf{y_{t-1}}+
\left[\begin{matrix}0.1&-0.2&0.5\\-0.4&0.1&-0.4\\0.1&0.3&-0.2\end{matrix}\right]\mathbf{y_{t-2}}+
\boldsymbol{\varepsilon_t}
}$$
$${\scriptsize
\boldsymbol{\varepsilon}_t\sim W.N.\left(\left[\begin{matrix}0.2^2&0&0\\0&0.2^2&0\\0&0&0.2^2\end{matrix}\right]\right)
}$$

とりあえずこのモデルに従うサンプルデータを作らないといけないですね。

まず、各係数行列や定数項などを変数として格納しておきます。

import numpy as np
import pandas as pd


# 各係数を変数に格納
A = np.array(
    [
        [1, 0, 0],
        [0.3, 1, 0],
        [-0.2, -0.3, 1]
    ]
)
A1 = np.array(
    [
        [-0.1, 0.3, 0.4],
        [0.2, -0.2, -0.3],
        [-0.3, 0.4, 0.2]
    ]
)
A2 = np.array(
    [
        [0.1, -0.2, 0.5],
        [-0.4, 0.1, -0.4],
        [0.1, 0.3, -0.2]
    ]
)
c = np.array([0.8, -0.5, 0.3])

これを使って、データを作ります。行列$A$の逆行列を両辺にかけて、モデルを誘導形に変形して順番に計算したらOKです。最初の2時点のデータは適当に設定して、初期のデータを捨ててます。

# 最初の2点のデータy_0, y_1を仮置き
df = pd.DataFrame(
    {
        "y0": [1, 1],
        "y1": [1, 1],
        "y2": [1, 1],
    }
)

# Aの逆行列
A_inv = np.linalg.inv(A)

# 乱数固定
np.random.seed(1)
for i in range(len(df), 550):
    df.loc[i] = A_inv@(c+A1@df.iloc[-1].T+A2@df.iloc[-2].T+np.random.normal(size=3)*0.2)

# 最初の方のデータを切り捨てる。
df = df.iloc[50:]
df.reset_index(inplace=True, drop=True)

これでデータができました。グラフは省略していますが、plotしてみると定常であることがわかります。

データができたので、statsmodelsで推定してみましょう。ドキュメントはこちらです。
参考: statsmodels.tsa.vector_ar.svar_model.SVAR — statsmodels
推定結果についてはこちら。
参考: statsmodels.tsa.vector_ar.svar_model.SVARResults — statsmodels

これの使い方は結構特殊です。左辺の行列のうち、推定したいパラメーターを文字列”E”とした行列を作って一緒に渡してあげる必要があります。対角成分は1です。前回の記事で書きましたが、各過程が外生性が高い順に並んでると仮定しているんので、下三角行列になるようにしています。具体的には次のようなコードになります。

import statsmodels.api as sm


# 上式のAの中で求めたい要素を"E"とした行列を生成しモデルに渡す。
A_param = np.array([[1, 0, 0], ["E", 1, 0], ["E", "E", 1]])

svar_model = sm.tsa.SVAR(df, svar_type="A", A=A_param)
svar_result = svar_model.fit(maxlags=2)

はい、これでできました。推定されたパラメーターが変数svar_resultに入っていますので、順番に見ていきましょう。

# 左辺の行列A
print(svar_result.A)
"""
[[ 1.          0.          0.        ]
 [ 0.34089688  1.          0.        ]
 [-0.20440697 -0.31127542  1.        ]]
"""

# 右辺の係数はcoefsにまとめて入っている
print(svar_result.coefs)
"""
[[[-0.12822445  0.30295997  0.43086499]
  [ 0.28516025 -0.30381663 -0.40460447]
  [-0.25660898  0.38056812  0.17305955]]

 [[ 0.05884724 -0.18589909  0.52278365]
  [-0.40304548  0.11406067 -0.60635189]
  [-0.07409054  0.28827388 -0.23204729]]]
"""

# 定数項
print(svar_result.intercept)
"""
[ 0.86592879 -0.82006429  0.28888304]
"""

そこそこの精度で推定できていますね。

VARモデルの時と同様に、summary()メソッドで結果を一覧表示することもできます。ただ、statsmodelsの現在のバージョン(0.13.5)のバグだと思うのですが、k_exog_userってプロパティを手動で設定しておかないと、AttributeError: ‘SVARResults’ object has no attribute ‘k_exog_user’ってエラーが出ます。とりあえず0か何か突っ込んで実行しましょう。

svar_result.k_exog_user=0
svar_result.summary()
"""
  Summary of Regression Results   
==================================
Model:                        SVAR
Method:                        OLS
Date:           Mon, 27, Feb, 2023
Time:                     00:29:40
--------------------------------------------------------------------
No. of Equations:         3.00000    BIC:                   -9.36310
Nobs:                     498.000    HQIC:                  -9.47097
Log likelihood:           276.728    FPE:                7.18704e-05
AIC:                     -9.54065    Det(Omega_mle):     6.89230e-05
--------------------------------------------------------------------
Results for equation y0
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.865929         0.038253           22.637           0.000
L1.y0        -0.128224         0.042999           -2.982           0.003
L1.y1         0.302960         0.037971            7.979           0.000
L1.y2         0.430865         0.041409           10.405           0.000
L2.y0         0.058847         0.043158            1.364           0.173
L2.y1        -0.185899         0.040139           -4.631           0.000
L2.y2         0.522784         0.035344           14.791           0.000
========================================================================

Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const        -0.820064         0.043394          -18.898           0.000
L1.y0         0.285160         0.048778            5.846           0.000
L1.y1        -0.303817         0.043074           -7.053           0.000
L1.y2        -0.404604         0.046974           -8.613           0.000
L2.y0        -0.403045         0.048958           -8.233           0.000
L2.y1         0.114061         0.045533            2.505           0.012
L2.y2        -0.606352         0.040094          -15.123           0.000
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.288883         0.041383            6.981           0.000
L1.y0        -0.256609         0.046517           -5.516           0.000
L1.y1         0.380568         0.041078            9.265           0.000
L1.y2         0.173060         0.044797            3.863           0.000
L2.y0        -0.074091         0.046688           -1.587           0.113
L2.y1         0.288274         0.043423            6.639           0.000
L2.y2        -0.232047         0.038236           -6.069           0.000
========================================================================

Correlation matrix of residuals
            y0        y1        y2
y0    1.000000 -0.300511  0.090861
y1   -0.300511  1.000000  0.269623
y2    0.090861  0.269623  1.000000
"""

正直、SVARモデルを使う時って、左辺の係数行列$A$が一番注目するところだと思うのですが、その情報が出てこないってところがイマイチですね。このsummary()はVARモデルほどは使わないと思いました。バグが放置されているのものそのせいかな？

最後の結果出力だけイマイチでしたが、必要な値は各属性を直接見れば取れますし、そこそ手軽に使えるモデルではあったのでVARモデルで力不足に感じることがあったらこれも思い出してみてください。

久しぶりに時系列の話です。以前、ベクトル自己回帰モデル(VAR)というのを紹介しました。
参考: ベクトル自己回帰モデル

これは要するに、時系列データのベクトルを、それより前の時点のベクトルで回帰する(線形和として表現する)ことによって説明しようっていうモデルでした。

これを実際に業務で使おうとすると、非常に厄介な問題が発生します。それは、同じ時点での値どうしの間にも関係があるということです。

例えば、何かのECサイトの分析をしてて、サイトの訪問者数、会員登録数、売上、の3つの時系列データがあったとした場合、VARの観点で言うと、過去の訪問者数、過去の会員登録数、過去の売上、からその日の各値を予測しようって言うのがVARモデルです。

そうなった時に、いやいや、会員登録数は「当日の訪問者数」の影響を受けるし、売上は「当日の会員登録数」の影響を受けるでしょうとなります。VARモデルではそう言う影響が考慮できません。ここを改善したのが構造VARモデル(Structural Vector Autoregressive Model)です。

時系列分析でいつも引き合いに出している、沖本先生の「経済・ファイナンスデータの計量時系列分析」では、4.6節(99〜101ページ)でさらっと紹介されています。2ページくらいです。

具体的に、3変数で最大ラグが2のSVARモデルを書き出すと下のような形になります。

$$\left\{\begin{align}y_{0,t}&=c_0&&& +\phi_{0,0,1}y_{0,t-1}+\phi_{0,1,1}y_{1,t-1}+\phi_{0,2,1}y_{2,t-1}+\phi_{0,0,2}y_{0,t-2}+\phi_{0,1,2}y_{1,t-2}+\phi_{0,2,2}y_{2,t-2}+\varepsilon_0\\
y_{1,t}&=c_1&-\phi_{1,0,0}y_{0,t}&&+\phi_{1,0,1}y_{0,t-1}+\phi_{1,1,1}y_{1,t-1}+\phi_{1,2,1}y_{2,t-1}+\phi_{1,0,2}y_{0,t-2}+\phi_{1,1,2}y_{1,t-2}+\phi_{1,2,2}y_{2,t-2}+\varepsilon_1\\
y_{2,t}&=c_2&-\phi_{2,0,0}y_{0,t}&-\phi_{2,1,0}y_{1,t}&+\phi_{2,0,1}y_{0,t-1}+\phi_{2,1,1}y_{1,t-1}+\phi_{2,2,1}y_{2,t-1}+\phi_{2,0,2}y_{0,t-2}+\phi_{2,1,2}y_{1,t-2}+\phi_{2,2,2}y_{2,t-2}+\varepsilon_2\\\end{align}\right.$$

後ろの方、ブログ幅からはみ出しましたね。見ての通りそこそこ巨大なモデルになります。
この、$y_{1,t}$の予測に$y_{0,t}$が使われていたり、$y_{2,t}$の予測に$y_{0,t},y_{1,t}$が使われているのが、VARモデルとの違いです。

逆に、$y_{0,t}$の説明には$y_{1,y}$や$y_{2,t}$は用いられてはいません。

これは変数が外生性が高い順に並んでいることを仮定しているためです。これを仮定せず、相互に影響し合うようなモデルも研究されてはいるようなのですが、実データから係数を推定するのが非常に難しくなるので、SVARモデルを使うときはこの仮定を置いておいた方が良い、と言うよりこれが仮定できる時に利用を検討した方が良いでしょう。

通常は、時刻$t$の時点の項を左辺に移行して行列表記するようです。(そのため、上の例でも移項を想定してマイナスつけときました。)

移行して、行列、ベクトルを記号に置き換えていくと、下のような式になります。$\mathbf{D}$は$n\times n$の対角行列とします。要するに撹乱項はそれぞれ相関を持たないとします。

$$
\mathbf{B}_0\mathbf{y}_t=\mathbf{c}+\mathbf{B}_1\mathbf{y}_{t-1}+\cdots+\mathbf{B}_p\mathbf{y}_{t-p}+\boldsymbol{\varepsilon}_t,\ \ \
\boldsymbol{\varepsilon}_t\sim W.N.(\mathbf{D})
$$

この形の式を構造形(Structual form)と呼びます。VARモデルとの違いは、左辺に$\mathbf{B}_0$が掛かってることですね。

この構造形ですが、実際にデータがあった時にこのまま係数を推定することが難しいと言う問題があります。そのため、両辺に$\mathbf{B}_0$をかけて次の形の式を考えます。

$$
\mathbf{y}_t=\mathbf{B}_0^{-1}\mathbf{c}+\mathbf{B}_0^{-1}\mathbf{B}_1\mathbf{y}_{t-1}+\cdots+\mathbf{B}_0^{-1}\mathbf{B}_p\mathbf{y}_{t-p}+\mathbf{B}_0^{-1}\boldsymbol{\varepsilon}_t
$$

この形を、誘導形(reduced form)と呼びます。撹乱項にも逆行列がかかってるので、誘導形の撹乱項の各成分には相関が生まれている点に気をつけてください。

この誘導形はVARモデルと全く同じなので、VARモデルを推定するのと同じ方法でパラメーターを推定することができます。

そして、ここからが問題なのですが、誘導形から構造形を求めることはそう簡単ではありません。(逆に構造形が事前にわかっていた場合に、誘導形を求めることは簡単です。上でやった通り、両辺に逆行列かけるだけだからです。)

その難しさの原因を説明をします。上の方で、変数が外生性が高い順に並んでると仮定して$\mathbf{B}$の一部の成分が0であることを仮定していましたが、もしその仮定がなかったとしましょう。すると、構造形の方が誘導形よりパラメーターが多くなってしまうのです。

構造形の方は、$\mathbf{y}_0$の係数の行列が、対角成分は1とわかってるので残り$n(n-1)$個、右辺は定数項が$n$個、右辺の各行列が全部で$pn^2$個、撹乱項が$n$個ので、合わせて、$n(n-1)+n+pn^2+n=n(n+1)+pn^2$個のパラメータを持ちます。

一方で誘導形の方は、左辺のパラメーターこそ消えますが、右辺は定数項が$n$個、右辺の各行列が全部で$pn^2$個、撹乱項が$n(n+1)/2$個で、合計$n+pn^2+n(n+1)/2$個しかパラメーターを持ちません。

その差、$n(n-1)/2$個だけ、構造形がパラメーターが多いので、誘導形が定まった時に構造形が決められなくなってしまいます。

そこで、誘導形から構造形を一意に定めるために、$n(n-1)/2$個の制約を課す必要が発生します。その制約として、$\mathbf{B}_0$の上三角成分(対角成分より上)を0と決めてしまうのが、変数が外生性が高い順に並んでいると言う仮定です。

個人的には、この仮定があったとしても誘導形から構造形を導くのってそこそこ難しいように感じますが、それでも理論上は算出が可能になりますね。

具体的にデータを用意してPythonを使ってSVARのパラメーター推定をやってと記事を続ける予定だったのですがここまでで結構長くなってしまったのと、数式のせいでそこそこ疲れてしまったので、それは次回の記事に回そうと思います。

VARの欠点をクリアしたモデルではありますが見ての通り巨大なので、かなりデータが多くないと推定がうまくいかなかったりして、なかなか期待ほど活躍しないのですが、こう言うモデルもあるってことを認識しておくとどこかで役に立つのかなと思います。