Rで計量時系列分析：VARモデルの基礎（多変量時系列モデル） - 渋谷駅前で働くデータサイエンティストのブログ

前回の記事では単変量の時系列までを扱いました。今回は多変量（ベクトル）時系列を記述するVARモデルとその周辺のポイントを取り上げます。

ということでしつこいですが、使用テキストはいつもの沖本本です。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

作者: 沖本竜義
出版社/メーカー: 朝倉書店
発売日: 2010/02/01
メディア: 単行本
購入: 4人クリック: 101回
この商品を含むブログ (6件) を見る

以下タイトルにのっとってRで各モデルの挙動を見ながらやっていきます。

必要なRパッケージ

{vars}をインストールして展開して下さい。{forecast}や{tseries}などは今回は特に使いません。

多変量における自己相関、定常性など

単変量時系列過程の際にさんざん自己相関やら定常性やらうるさく言っておいて、まさか多変量にした時にガン無視ってわけにもいきませんので（笑）一応触れておきます。沖本本pp.74-76をご参照下さい。

まず、ベクトル $\mathbf{y}$ の期待値ベクトルは

$E(\mathbf{y}_t) = [ E(y_{1t}),\cdots, E(y_{nt}) ]'$

で定義されます。当たり前ながら、ベクトルの期待値は各成分の期待値から成るベクトルである、ということですね。次に自己相関はどうなるかというと、これは行列になります。k次自己共分散行列は

$Cov(\mathbf{y}_t,\mathbf{y}_{t-k}) = [Cov(y_{it},y_{j,t-k})]_{ij}$

$= \left( \begin{matrix} Cov(y_{1t},y_{1,t-k}) & \ldots & Cov(y_{1t},y_{n,t-k}) \\ \vdots & \ddots & \vdots \\ Cov(y_{nt},y_{1,t-k}) & \ldots & Cov(y_{nt},y_{n,t-k}) \end{matrix} \right)$

で定義されます。つまり $(i,j)$ 成分が $y_it$ と $y_{j,t-k}$ の共分散に等しい $n \times n$ 行列です。このとき、k次自己共分散行列の対角成分は各変数のk次自己共分散に等しくなっています。また単変量のときと同様に、自己共分散行列はkの関数として表されます。

期待値と自己共分散関数がtに依存しないとき*1、ベクトル過程は弱定常（共分散定常）と言われます。基本的にここからの議論はこの弱定常性を仮定していくものとします。

なお、単変量の時にも出てきたホワイトノイズはベクトル過程では以下のように定義されます。

$E(\mathbf{\epsilon}_t) = \mathbf{0}$

$E(\mathbf{\epsilon_t}\mathbf{\epsilon'_{t-k}}) = \left \{ \begin{array} \mathbf{\Sigma}, k=0 \\ \mathbf{0}, k \neq 0 \end{array} \right.$

この場合 $\mathbf{\epsilon_t}$ はベクトルホワイトノイズと呼ばれます。ちなみにここで一つ注意しなければならないのは、 $\mathbf{\Sigma}$ は $n \times n$ 正定値行列であり、対角行列である必要はないという点。 $\mathbf{\epsilon_t}$ は異時点においてはどの成分も相関を持たない一方で、同時点では各成分が相関を持つこともあることに気を付ける必要があります。なお、 $\mathbf{\epsilon_t}$ が分散共分散行列 $\mathbf{\Sigma}$ のベクトルホワイトノイズであるとき、以下 $\mathbf{\epsilon_t} \sim W.N.(\mathbf{\Sigma})$ と表記することにします。

VARモデル

・・・ふー、ようやくVARまで来た（笑）。ここからは沖本本pp.76-79の内容に入ります。

ところで何故VAR（Vector Autoregressive: ベクトル自己回帰）モデルと呼ぶのか？という点についてですが、理由は簡単でこれまで取り上げてきた単変量時系列過程 $y_t$ を単純に $n \times 1$ の列ベクトル $\mathbf{y} = (y_{1t}, y_{2t}, \cdots , y_{nt})'$ の形に並べることで表現する、というものだからです。

なので、VARと言ったら単純に「多変量時系列モデル」と脳内変換しても（通常であれば）差し支えないと思います。

VARモデルの定義

VAR(p)モデルを $\mathbf{y}_t$ を定数と自身のp期の過去の値に回帰したモデルとすると、

$\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$ 係数行列）

これでは分かりにくいので、例えば2変量VAR(1)モデルに書き下すとこんな感じになります。

$\left\{ \begin{array} ~y_{1t}=c_1+\phi_{11} y_{1,t-1}+\phi_{12} y_{2,t-1}+\epsilon_{1t} \\ y_{2t}=c_2+\phi_{21} y_{1,t-1}+\phi_{22} y_{2,t-1}+\epsilon_{2t} \end{array} \right.$
$\left ( \begin{array} ~ \epsilon_{1t} \\ \epsilon_{2t} \end{array} \right ) \sim W.N.(\Sigma)$

n変量VAR(p)モデルはn本の回帰式からなり、それぞれの回帰式は各変数を定数と全変数のp期間の過去の値に回帰した形になっています。そこでこのモデルが含むパラメータの個数を考えてみると、1本の回帰式が定数を含めてnp+1個の係数を含むので、n本になればn(np+1)個のパラメータになります。で、実は $\mathbf{\Sigma}$ がn(n+1)/2個のパラメータを持っているので、全体で見ると合計n(np+1)+n(n+1)/2個のパラメータを持つ、比較的大きなモデルになるとも言えるわけです。この「大きなモデル」というのが実は重要になってきます。

VARモデルの定常条件など

実は、ARモデルの際に用いたAR特性方程式であったりユール・ウォーカー方程式を行列版に拡張したもので、ほぼ同じような制約（全ての解の絶対値が1未満）を課した形で得ることができます。

VARモデルの推定

ARMAの時とは異なり、VARモデルの各方程式は同時点のその他の変数は含まないので、同時方程式モデル(simultaneous equation model)ではない*2とされます。ただし、各方程式は誤差項の相関を通じて関係しており、見かけ上無関係な回帰(SUR)モデル(seemingly unrelated regression model)と呼ばれるようです。

一般にSURモデルを推定するためには誤差項同士の相関を考慮に入れなければならず、全ての回帰式を同時に推定する必要があるそうですが、VARモデルは全ての回帰式が同一の説明変数を持つという特殊性があるため、実は各方程式を個別にOLSによって推定するだけで良いということが分かっています。しかも撹乱項 $\epsilon_t$ が多変量正規分布に従うと仮定した場合、OLS係数推定量は最尤推定量とも一致します。ということで、基本的にはOLSでただ解くだけでVARモデル推定が出来てしまうわけです*3。

なお、VARモデルの次数選択も基本的には情報量基準によります。ただ、一般にVARモデルが使いたい場面というのは割とグローバルなファイナンシャルデータだったりする上に、次数pの候補をどこまで遡るかによって容易にAICなりSICなりの極小値が変動してしまうという問題もあるので、必ずしも情報量基準に基づく決め方が万能とは限らないという点に留意しましょう。

RでVARモデルを推定してみる

あー、長かった（汗）。ここからは普通にRでVARモデルを推定していきます。手っ取り早く{vars}パッケージのサンプルデータであるCanadaを用いてみましょう。これはカナダの1980～2000年の主要なマクロ経済指数を、適切な変数変換を行って4つ抜き出してきたものです。prodが労働生産性、eが雇用、Uが失業率、rwが実質賃金を表しています。

基本的に、{vars}パッケージでVAR(p)モデル推定を行う場合の手順は

まず次数pをVARselect()関数で推定する
求めた次数pをVAR()関数に代入してVAR(p)モデル係数を推定する

という流れになっています。全然難しくないです。上記の制約条件を全て理解した上で適切にどちらもこなせれば、ですが（笑）。

実際にRで実行してみるとこんな感じになります。

> data(Canada)
> VARselect(Canada,lag.max=4) # 四半期ごとのデータなので最大次数を1年 = 4Q = 4とした
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      2      2      3 

$criteria
                 1            2            3            4
AIC(n) -5.70832549 -6.238366753 -6.359392786 -6.119193878
HQ(n)  -5.46956983 -5.808606566 -5.738628071 -5.307424635
SC(n)  -5.11281884 -5.166454767 -4.811075473 -4.094471239
FPE(n)  0.00332039  0.001960529  0.001750655  0.002258874

> Canada.var<-VAR(Canada,p=VARselect(Canada,lag.max=4)$selection[1])
# pに上記の推定結果をそのまま入力

> summary(Canada.var)
# 4変量VAR(3)モデルの詳細が以下に表示される

VAR Estimation Results:
========================= 
# VARモデル概要
Endogenous variables: e, prod, rw, U 
Deterministic variables: const 
Sample size: 81 
Log Likelihood: -150.609 
Roots of the characteristic polynomial:
1.004 0.9283 0.9283 0.7437 0.7437 0.6043 0.6043 0.5355 0.5355 0.2258 0.2258 0.1607
Call:
VAR(y = Canada, p = VARselect(Canada, lag.max = 4)$selection[1])

# eのモデル推定結果

Estimation results for equation e: 
================================== 
e = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 

          Estimate Std. Error t value Pr(>|t|)    
e.l1       1.75274    0.15082  11.622  < 2e-16 ***
prod.l1    0.16962    0.06228   2.723 0.008204 ** 
rw.l1     -0.08260    0.05277  -1.565 0.122180    
U.l1       0.09952    0.19747   0.504 0.615915    
e.l2      -1.18385    0.23517  -5.034 3.75e-06 ***
prod.l2   -0.10574    0.09425  -1.122 0.265858    
rw.l2     -0.02439    0.06957  -0.351 0.727032    
U.l2      -0.05077    0.24534  -0.207 0.836667    
e.l3       0.58725    0.16431   3.574 0.000652 ***
prod.l3    0.01054    0.06384   0.165 0.869371    
rw.l3      0.03824    0.05365   0.713 0.478450    
U.l3       0.34139    0.20530   1.663 0.100938    
const   -150.68737   61.00889  -2.470 0.016029 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.3399 on 68 degrees of freedom
Multiple R-Squared: 0.9988,	Adjusted R-squared: 0.9985 
F-statistic:  4554 on 12 and 68 DF,  p-value: < 2.2e-16 

# prodのモデル推定結果

Estimation results for equation prod: 
===================================== 
prod = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 

          Estimate Std. Error t value Pr(>|t|)    
e.l1      -0.14880    0.28913  -0.515   0.6085    
prod.l1    1.14799    0.11940   9.615 2.65e-14 ***
rw.l1      0.02359    0.10117   0.233   0.8163    
U.l1      -0.65814    0.37857  -1.739   0.0866 .  
e.l2      -0.18165    0.45083  -0.403   0.6883    
prod.l2   -0.19627    0.18069  -1.086   0.2812    
rw.l2     -0.20337    0.13337  -1.525   0.1319    
U.l2       0.82237    0.47034   1.748   0.0849 .  
e.l3       0.57495    0.31499   1.825   0.0723 .  
prod.l3    0.04415    0.12239   0.361   0.7194    
rw.l3      0.09337    0.10285   0.908   0.3672    
U.l3       0.40078    0.39357   1.018   0.3121    
const   -195.86985  116.95813  -1.675   0.0986 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.6515 on 68 degrees of freedom
Multiple R-Squared:  0.98,	Adjusted R-squared: 0.9765 
F-statistic: 277.5 on 12 and 68 DF,  p-value: < 2.2e-16 

# rwのモデル推定結果

Estimation results for equation rw: 
=================================== 
rw = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 

          Estimate Std. Error t value Pr(>|t|)    
e.l1    -4.716e-01  3.373e-01  -1.398    0.167    
prod.l1 -6.500e-02  1.393e-01  -0.467    0.642    
rw.l1    9.091e-01  1.180e-01   7.702 7.63e-11 ***
U.l1    -7.941e-04  4.417e-01  -0.002    0.999    
e.l2     6.667e-01  5.260e-01   1.268    0.209    
prod.l2 -2.164e-01  2.108e-01  -1.027    0.308    
rw.l2   -1.457e-01  1.556e-01  -0.936    0.353    
U.l2    -3.014e-01  5.487e-01  -0.549    0.585    
e.l3    -1.289e-01  3.675e-01  -0.351    0.727    
prod.l3  2.140e-01  1.428e-01   1.498    0.139    
rw.l3    1.902e-01  1.200e-01   1.585    0.118    
U.l3     1.506e-01  4.592e-01   0.328    0.744    
const   -1.167e+01  1.365e+02  -0.086    0.932    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.7601 on 68 degrees of freedom
Multiple R-Squared: 0.9989,	Adjusted R-squared: 0.9987 
F-statistic:  5239 on 12 and 68 DF,  p-value: < 2.2e-16 

# Uのモデル推定結果

Estimation results for equation U: 
================================== 
U = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 

         Estimate Std. Error t value Pr(>|t|)    
e.l1     -0.61773    0.12508  -4.939 5.39e-06 ***
prod.l1  -0.09778    0.05165  -1.893 0.062614 .  
rw.l1     0.01455    0.04377   0.332 0.740601    
U.l1      0.65976    0.16378   4.028 0.000144 ***
e.l2      0.51811    0.19504   2.656 0.009830 ** 
prod.l2   0.08799    0.07817   1.126 0.264279    
rw.l2     0.06993    0.05770   1.212 0.229700    
U.l2     -0.08099    0.20348  -0.398 0.691865    
e.l3     -0.03006    0.13627  -0.221 0.826069    
prod.l3  -0.01092    0.05295  -0.206 0.837180    
rw.l3    -0.03909    0.04450  -0.879 0.382733    
U.l3      0.06684    0.17027   0.393 0.695858    
const   114.36732   50.59802   2.260 0.027008 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.2819 on 68 degrees of freedom
Multiple R-Squared: 0.9736,	Adjusted R-squared: 0.969 
F-statistic: 209.2 on 12 and 68 DF,  p-value: < 2.2e-16 


# その他推定残差など

Covariance matrix of residuals:
            e     prod       rw        U
e     0.11550 -0.03161 -0.03681 -0.07034
prod -0.03161  0.42449  0.05589  0.01494
rw   -0.03681  0.05589  0.57780  0.03660
U    -0.07034  0.01494  0.03660  0.07945

Correlation matrix of residuals:
           e     prod      rw        U
e     1.0000 -0.14276 -0.1425 -0.73426
prod -0.1428  1.00000  0.1129  0.08136
rw   -0.1425  0.11286  1.0000  0.17084
U    -0.7343  0.08136  0.1708  1.00000

モデル推定の詳細を、可視化して見ることもできます。

plot(Canada.var)

f:id:TJO:20130725164222p:plain

f:id:TJO:20130725164255p:plain

f:id:TJO:20130725164349p:plain

f:id:TJO:20130725164417p:plain

また得られたVARモデルに基づいて、予測を行うこともできます。

> Canada.pred<-predict(Canada.var,n.ahead=20,ci=0.95)
> plot(Canada.pred)

f:id:TJO:20130725164143p:plain

なおシミュレーションを行いたい場合ですが、VAR.sim(){tsDyn}関数でVAR(p)過程をランダムに発生させることができます。

B1<-matrix(c(0.7, 0.2, 0.2, 0.7), 2) # VARモデルの係数を適当に決める
x<-VAR.sim(B=B1,n=100,include="const") # 代入してシミュレート
x.p<-VARselect(x.var,lag.max=10,include="const")$selection[1]
x.var<-VAR(x,p=x.p,type="const")

厳密には色々とVARモデルには制約があるのですが、今回はそこはスキップします*4。

なお、VARモデル同様にVMA / VARMA / VARIMAモデルも存在しますが、そもそもn変量VAR(p)モデルが合計n(np+1)+n(n+1)/2個ものパラメータを持つ大きなモデルでかなり複雑な振る舞いまで記述することが可能な上に、それらのもっと高度なモデルの推定が*5困難と言うこともあって、よほどのことがない限りはVARモデルで十分ということにされているようです。