VARモデル補遺（備忘録） - 渋谷駅前で働くデータサイエンティストのブログ

もう9年も前のことですが、沖本本をベースとした計量時系列分析のシリーズ記事を書いていたことがあります。その中で、密かに今でも自分が読み返すことがあるのがVAR（ベクトル自己回帰）モデル関連の記事です。

なのですが、仕事なり趣味なりでVARモデルを触っていると「あれ、これってどうなってたんだっけ」という事項が幾つか出てきて、しかも上記の自分のブログの過去記事を当たっても出てこないケースがちらほらあるんですね。

ということで、今回の記事では~~ネタ切れで新しく書くことが思い付かないので~~VARモデル周りで「最近になって調べて行き当たったこと」を備忘録的に補遺として書き留めておくことにします。とはいえ微妙に技術的な話題を含むので、いつもながらですが誤りなどありましたらコメントなどでご指摘くださると助かります。

CanadaデータセットでVARモデルを推定しておく

面倒なので、{vars}パッケージに同梱されているCanadaデータセットでVARモデルを推定しておきます。これはカナダの1980-2000年の各四半期における生産性(prod)・雇用(e)・失業率(U)・実質賃金(rw)の4経済指標を時系列データとしてまとめたものです。

library(vars)

data(Canada)

VAR.select <- VARselect(Canada, lag.max = 10, type = 'both', season = 4)
Canada.var <- VAR(Canada, p = VAR.select$selection[1], lag.max = 10,
                  type = 'both', season = 4)

Canada.pred <- predict(Canada.var, n.ahead = 8)
Canada.irf <- irf(Canada.var, impulse = 'U', n.ahead = 8, ci = 0.95, ortho = T)

特に理由はないですが、ついでに8期先予測と直交化インパルス応答関数（一旦インパルスはU即ち失業率とした）も推定しておきました。それぞれプロットすると以下のようになります。

plot(Canada.pred)
plot(Canada.irf)

この辺のプロットは以前の記事をお読みになったことのある方には馴染み深いかと思います。

インパルス応答関数推定における「撹乱項の1標準偏差」の求め方

ところで、上記の過去記事ではインパルス応答関数推定について以下のように書きました。

Rで主に用いることになる直交化インパルス応答関数を求める際に重要になるのは、撹乱項 $\mathbf{\epsilon}_t$ の分散共分散行列 $\mathbf{\Sigma}$ です。この $\mathbf{\Sigma}$ は、
$\mathbf{\Sigma} = \mathbf{ADA'}$
（ただし $\mathbf{A}$ は対角成分が1に等しい下三角行列、 $\mathbf{D}$ は対角行列）
と三角分解することができ、このとき直交化撹乱項 $\mathbf{u}_t$ を
$\mathbf{u}_t = \mathbf{A}^{-1} \epsilon_t$
と定義し、各変数固有の変動を表すとみなすことができます。なお、 $\mathbf{u}_t$ は互いに無相関です（ただし $Var(\mathbf{u}_t)=\mathbf{D}$ ：証明は沖本本を参照のこと）。

$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,\cdots$
と表されます。1単位ではなく1標準偏差のショックを与えたいときは、上の式にそのまま $\sqrt{Var(u_{jt})}$ を乗じれば良いだけです。

なのですが、Rでは三角分解の代わりにコレスキー分解を用いて撹乱項を分解し、その分解した撹乱項に1単位のショックを与えてインパルス応答関数を計算するという方法を用いています。これについてちょっと書いておきましょう。 $\mathbf{D}^{1/2}$ を(j,j)成分が $u_{jt}$ の標準偏差に等しい対角行列であるとすると、
$\mathbf{\Sigma} = \mathbf{ADA'}$
は、 $\mathbf{P} \equiv \mathbf{AD}^{1/2}$ として
$\mathbf{\Sigma} = \mathbf{AD}^{1/2} \mathbf{D}^{1/2} \mathbf{A'} = \mathbf{PP'}$
と書けます。この式は正定値行列 $\mathbf{\Sigma}$ のコレスキー分解になっていて、 $\mathbf{P}$ はコレスキー因子と呼ばれます。この $\mathbf{P}$ を用いれば、撹乱項を
$\mathbf{v}_t = \mathbf{P}^{-1}\mathbf{\epsilon}_t = \mathbf{D}^{-1/2}\mathbf{A}^{-1}\mathbf{\epsilon}_t = \mathbf{D}^{-1/2}\mathbf{u}_t$
と分解することもできます。なお $Var(\mathbf{u}_t)=\mathbf{D}$ であることから、 $v_{jt}$ は $u_{jt}$ をその標準偏差 $\sqrt{Var(u_{jt})}$ で割ったものになっています。これは言い換えると $v_{jt}$ における1単位の増加が $u_{jt}$ における1標準偏差の増加に等しいということであり、同じ要領で言えば $v_{jt}$ に1単位のショックを与えて計算したインパルス応答関数は、 $u_{jt}$ に1標準偏差のショックを与えて計算したインパルス応答関数に等しくなるということです。

一般には1標準偏差あたりのショックを用いることが多いことから、Rはじめ多くの統計分析ツールではコレスキー分解が使われているのだろう、とは沖本先生のコメントです。

ということで、vars::varを使って得られたVARモデルに対してvars::irfを用いると基本的には直交化インパルス応答関数の推定結果が得られます。ところが、往々にして実務データ分析ではその「撹乱項の1標準偏差」が一体どれくらいの大きさ(magnitude)なのか？という情報が必要になることがあります。これが意外と分からないんですよね。困ったなぁと思って色々ググっていたら、{vars}パッケージの作者に問い合わせた上で一通り試された方のブログ記事が見つかりました。

結論から言えば、今回のCanadaデータセットに対するVARモデルで失業率Uにショックを与えると仮定した場合の「撹乱項の1標準偏差」は、以下のようにして得られます。

shock_U <- summary(Canada.var)$varresult$U$sigma

9年前は全然気付いてなかったんですが、varestオブジェクトのvarresult以下には結構リッチな情報が詰まっているんですね。このことに気付いたことで、僕はもう一つ調べたいことが出てきたのでした。

VARモデルにおいて特定の変数に対する「他変数からの寄与分」を求める方法

きっかけは3年前にある学術集会でなされた研究報告で、要は「VARモデルでもMMM (Media Mix Model)同様に各変数からの寄与分を積み上げプロットとして可視化することができるはず」という話です。

$\mathbf{y_t}=\mathbf{c}+\mathbf{\Phi_1}\mathbf{y_{t-1}}+\cdots+\mathbf{\Phi_p}\mathbf{y_{t-p}}+\mathbf{\epsilon_t}, \mathbf{\epsilon_t} \sim W.N.(\mathbf{\Sigma})$
（ただし $\mathbf{c}$ は $n \times 1$ 定数ベクトル、 $\mathbf{\Phi_i}$ は $n \times n$ 係数行列）

というVARモデルの定義式を考えれば、それを計算するのは本来ならば容易なはずです。ところが、それを計算するためのデータやformulaってどこに格納されているんだっけ？とずっと思っていたのでした。

ここまで書けばもうお分かりの方も多いかと思いますが、それこそがvarestオブジェクトのvarresult以下なんですね。varest$varresultにはdatamatというデータ行列などが含まれているんですが、一旦summaryをかけてからvarresult$Uを見ると以下のようになっています。

summary(Canada.var)$varresult$U
#R> 
#R> Call:
#R> lm(formula = y ~ -1 + ., data = datamat)
#R> 
#R> Residuals:
#R>      Min       1Q   Median       3Q      Max 
#R> -0.57469 -0.18600  0.01721  0.15493  0.67160 
#R> 
#R> Coefficients:
#R>           Estimate Std. Error t value Pr(>|t|)    
#R> e.l1     -0.638553   0.127210  -5.020 4.39e-06 ***
#R> prod.l1  -0.110992   0.054195  -2.048 0.044665 *  
#R> rw.l1     0.008921   0.050866   0.175 0.861332    
#R> U.l1      0.623402   0.168646   3.697 0.000455 ***
#R> e.l2      0.539000   0.198172   2.720 0.008399 ** 
#R> prod.l2   0.081791   0.080207   1.020 0.311687    
#R> rw.l2     0.047954   0.066321   0.723 0.472280    
#R> U.l2     -0.076632   0.210502  -0.364 0.717025    
#R> e.l3     -0.068917   0.139933  -0.492 0.624052    
#R> prod.l3  -0.023500   0.055155  -0.426 0.671491    
#R> rw.l3    -0.015255   0.050986  -0.299 0.765762    
#R> U.l3      0.026064   0.175914   0.148 0.882678    
#R> const   165.486347  61.326740   2.698 0.008898 ** 
#R> trend     0.020595   0.013843   1.488 0.141729    
#R> sd1      -0.041149   0.100849  -0.408 0.684616    
#R> sd2       0.025118   0.102875   0.244 0.807889    
#R> sd3       0.030293   0.094624   0.320 0.749907    
#R> ---
#R> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#R> 
#R> Residual standard error: 0.2845 on 64 degrees of freedom
#R> Multiple R-squared:  0.9747,	Adjusted R-squared:  0.9684 
#R> F-statistic: 154.1 on 16 and 64 DF,  p-value: < 2.2e-16

そこでまずUを目的変数とするOLS線形モデルを推定します。これは以下のようにすれば得られます。

# Fit an individual model for U

datamat <- Canada.var$datamat[, -c(1:3)]
names(datamat)[1] <- 'y'
fit <- lm(y ~ -1 + ., data = datamat)

あとは、ちまちまと他変数それぞれの寄与分をpredictで算出するだけです。これが随分と面倒で、多分間違いなくずっともっと簡単にやれる方法があるはずなんですが、僕のRコーディング力がショボいだけなのでご容赦くださいorz

# Extract contribution from other variables

datamat_e <- datamat
datamat_e[, c(1, 3:5, 7:9, 11:18)] <- 0
ctrb_e <- predict(fit, newdata = datamat_e)

datamat_prod <- datamat
datamat_prod[, c(1:2, 4:6, 8:10, 12:18)] <- 0
ctrb_prod <- predict(fit, newdata = datamat_prod)

datamat_rw <- datamat
datamat_rw[, c(1:3, 5:7, 9:11, 13:18)] <- 0
ctrb_rw <- predict(fit, newdata = datamat_rw)

datamat_U <- datamat
datamat_U[, c(1:4, 6:8, 10:12, 14:18)] <- 0
ctrb_U <- predict(fit, newdata = datamat_U)

datamat_others <- datamat
datamat_others[, 1:13] <- 0
ctrb_others <- predict(fit, newdata = datamat_others)

# Stacked plot: is this right?

matplot(cbind(ctrb_e + ctrb_others,
              ctrb_e + ctrb_prod + ctrb_others,
              ctrb_e + ctrb_prod + ctrb_rw + ctrb_others,
              ctrb_e + ctrb_prod + ctrb_rw + ctrb_U + ctrb_others,
              datamat$y),
        type = 'l', lty = 1, ylab = '')

切片項（定数項＋トレンド項＋標準偏差）の扱いをどうするか迷ったんですが、一旦「先に足しておく」ようにしてあります（足さない方が正しい気もしますが）。ということで、面倒ですがこんな感じで実務データ分析でリクエストされがちなVARモデルの各種指標が求まりますよ、というお話でした。