前回に続いて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モデルで力不足に感じることがあったらこれも思い出してみてください。