Rで計量時系列分析：単位根過程、見せかけの回帰、共和分、ベクトル誤差修正モデル

前回の記事ではVARモデルに基づく様々な計量時系列分析手法を取り上げました。今回はいよいよ現実世界の時系列データを扱う上では避けて通れない、単位根過程とそれにまつわる様々な問題とその解決策について触れてみようと思います。

ということでもはや毎回恒例になってますが、使用テキストはいつもの沖本本です。

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

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

ただしこの辺の数理科学的な原理・根拠が気になる人は、Hamilton本も脇に置いておいた方が良いでしょう。

Time Series Analysis

作者: James D. Hamilton
出版社/メーカー: Princeton Univ Pr
発売日: 1994/01/11
メディア: ハードカバー
購入: 1人クリック: 5回
この商品を含むブログ (8件) を見る

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

必要なRパッケージ

{tseries}, {vars}パッケージをインストールしておきましょう。なお{urca}は{vars}の依存パッケージなので{vars}インストール時に自動的に入ります。それから、演習用データを取ってきたいので{RFinanceYJ}もインストールしておいて下さい。

単位根過程

以前の記事でもちろっと触れましたが、今回ちょっとしっかり説明したいと思います。

単位根過程の基礎

これまでは、基本的には常に定常過程を扱ってきました。しかしながら、「トレンドがない」「平均に回帰する」という定常過程の代表的な性質に従う時系列データというのは、経済・ファイナンスの世界で言うと実はむしろ少数派だったりします。

例えば株価は常に「直前の値に対してランダムにブレる値を足していく」ことで得られる時系列で平均回帰性なんか到底あるとは思えません。またGDPや物価などは平均的に一定の割合で成長していくので、対数系列を見てもはっきりとしたトレンドを示すのが常です。

ということで、こういう「あるトレンドを伴って変動する時系列データ」を表現する過程が必要になる、というわけです。そこで登場するのが単位根過程。沖本本pp.105-106にそのものの定義と関連する定義群が載っています。

定義5.1（単位根過程）　原系列 $y_t$ が非定常過程であり、差分系列 $\Delta y_t = y_t - y_{t-1}$ が定常過程であるとき、過程は単位根過程(unit root process)といわれる。

ちなみに何故「単位根」と呼ばれるかというと、AR特性方程式が $z=1$ という解を1つ持つためです。またその性質上呼び名が沢山あり、差分定常過程、1次和分過程（ $y_t \sim I(1)$ と表記される：n次差分したときに定常過程になる場合は勿論n次和分過程 $I(n)$ と表記される）、さらに単位根過程の差分系列が定常かつ反転可能なARMA(p,q)過程となるときは、次数(p,1,q)の自己回帰和分移動平均過程もしくはARIMA(p,1,q)過程と呼ばれます。

代表的な和分過程として知られているのが、ランダムウォークです。沖本本p.106によれば、

定義5.3（ランダムウォーク）　過程 $y_t$ が
$y_t = \delta + y_{t-1} + \epsilon_t, \epsilon_t \sim iid(0,\sigma^2)$ 　　(5.1)
と表現されるとき、 $y_t$ はランダムウォーク(random walk)と呼ばれる。ただし、 $y_0 = 0$ とする。また定数項 $\delta$ はドリフト率と呼ばれ、(5.1)はドリフト率 $\delta$ のランダムウォークと呼ばれることもある。

このランダムウォークの差分系列は $\Delta y_t = \delta + \epsilon_t$ と表せるので、定常系列であることが容易にわかりますね。

なお、ランダムウォークをRで表現するのは非常に簡単です。

> rw<-cumsum(rnorm(300))
> plot.ts(rw,lwd=2)

f:id:TJO:20130816101244p:plain

平均回帰していないことが一目見て分かるかと思います。

単位根過程のトレンド

単位根過程の重要な性質として、まず線形トレンドを記述できるということがあります。まぁこれはそもそもの「和分」表現を見れば自明かと思いますが、実際に繰り返し代入していった結果が沖本本p.107にあって、

$y_t = \delta + y_{t-1} + \epsilon_t = 2\delta + y_{t-2} + \epsilon_{t-1} + \epsilon_t = \cdots = \delta t + \nu_t$

ただし $\nu_t \equiv \epsilon_1 + \epsilon_2 + \cdots + \epsilon_t$ であり、確率的トレンドと呼ばれます。これで線形トレンドを記述できるというわけです。

ところで、実は線形トレンドを記述できるモデルとして別にトレンド定常過程と呼ばれるものもあります。沖本本p.107によれば、

定義5.5（トレンド定常過程）　 $x_t$ を定常過程として、
$y_t = \delta t + x_t$
と表される過程はトレンド定常過程(trend stationary process)と呼ばれる。

単位根過程とこのトレンド定常過程、似ているようで実は全然似ていません。沖本本p.108に詳しく書かれている通り、トレンド定常過程はトレンドに沿って一定の範囲内でばらつくのに対し、単位根過程は確率的トレンドに従って不確実性を線形的に増大させる*1わけです。

この原理的な両者の差はかなり大きくて、例えば沖本本pp.110-111にもあるようにインパルス応答も単位根過程なら経過時間を問わず常に「1」になる一方で、トレンド定常過程においては普通の定常過程同様に減衰していくことが知られています。沖本本p.111ではGDPのモデル化にどちらを使うかによって、景気刺激策の可否の判断が180度逆になりかねないという極端なケースを挙げて、単位根過程であるかどうかの判定の大事さを説いています。

単位根検定

ここではひとまず単位根検定そのものに話を絞ります。コンセプトとしては、真の過程をAR(p)モデルと仮定し、過程が単位根AR(p)過程であるという帰無仮説を、過程が定常AR(p)過程であるという対立仮説に対して検定するというものです。

沖本本pp.111-118に詳しい解説が書かれているので、原理的側面についてはこのブログではほぼ全面的に割愛します*2。ただし重要なポイントとして理解しておくべきなのは、

単位根過程は平均回帰性がないので、「標本平均が真の平均に漸近していく」的な大数の法則や中心極限定理が使えず、そのためt検定も使えない
そこでランダムウォークの連続時間極限を表すブラウン運動を用いて漸近分布を作り、これを検定統計量として使う

という2点です。一応、ここでブラウン運動の定義を書いておきます。沖本本pp.114-115によれば、

定義5.6（標準ブラウン運動）　標準ブラウン運動(standard Brownian motion)は、以下の性質を満たす[0,1]で定義されたスカラー連続確率過程 $W(t)$ である。

(1) $W(0) = 0$

(2) 任意の時点 $0 \leq t_1 < t_2 < \cdots < t_k \leq 1$ において、 $W(t_2) - W(t_1), W(t_3) - W(t_2), \cdots, W(t_k) - W(t_{k-1})$ は互いに独立な正規分布に従い、 $W(s) - W(t) \sim N(0,s-t)$ である。

(3) 任意の与えられた実現値において、 $W(t)$ は確率1でtに関して連続である。

ちなみにその定義上、 $\int_0^1 W(r)dr$ のような表現を使う機会が多いことを覚えておいた方が良いと思います。

ともあれ、ブラウン運動に基づく漸近分布を用いた検定統計量によって、真の過程をAR(1)過程と仮定して単位根検定を行う手法をDicky-Fuller(DF)検定と呼び、これをAR(p)過程に拡張したものが拡張DF(Augmented DF: ADF)検定です。Rではadf.test(){tseries}関数として実装されています。

> x <- rnorm(1000)  # サンプルとして定常過程を使う
> adf.test(x)

	Augmented Dickey-Fuller Test

data:  x 
Dickey-Fuller = -9.5759, Lag order = 9, p-value = 0.01
alternative hypothesis: stationary 
# 帰無仮説が棄却された、すなわちこれは定常過程という判定

 警告メッセージ： 
In adf.test(x) : p-value smaller than printed p-value
> 
> y <- diffinv(x) # 単位根を含むように変えてみた
> adf.test(y)

	Augmented Dickey-Fuller Test

data:  y 
Dickey-Fuller = -2.0187, Lag order = 9, p-value = 0.5704
alternative hypothesis: stationary 
# 帰無仮説が採択された、すなわちこれは単位根過程

なお沖本本ではPhillips-Perron(PP)検定についても触れています。こちらはAR(p)過程までを守備範囲としていたADF検定とは異なり、差分時系列 $u_t$ に対してより一般的に自己相関や分散不均一性まで考慮した検定手法です。Rではpp.test(){tseries}関数で実装されています。

> x <- rnorm(1000)  # サンプルデータは定常過程
> pp.test(x)

	Phillips-Perron Unit Root Test

data:  x 
Dickey-Fuller Z(alpha) = -964.483, Truncation lag parameter = 7,
p-value = 0.01
alternative hypothesis: stationary 
# 定常過程という判定

 警告メッセージ： 
In pp.test(x) : p-value smaller than printed p-value
> 
> y <- cumsum(x)  # 単位根過程に変えてみた
> pp.test(y)

	Phillips-Perron Unit Root Test

data:  y 
Dickey-Fuller Z(alpha) = -4.6724, Truncation lag parameter = 7,
p-value = 0.8491
alternative hypothesis: stationary 
# 単位根過程だという判定になった

ところで、単位根AR過程のモデル推定に当たり「単位根を持つことを無視したらどうなるか？」ということが沖本本pp.120-122に書かれています。結論から言うと「そんなに問題はない」が、「非標準的な漸近理論に基づくポイントがあるためF検定が使えず」、「特にVARモデルに拡張した場合はGranger因果性検定が使えない」。

単位根VAR過程の場合は、差分過程に直さない限りGranger因果は使えないという問題があるという点に、特に注意が必要とみて良いでしょう。

見せかけの回帰

いよいよ来ました、見せかけの回帰です。かなり以前の記事で一度紹介していますが、あらためておさらいしてみましょう。

Rで実際に見てみる

せっかくなので、[twitter:@yokkuns]さんが先日CRANリポジトリに最新版を苦労してupして下さった{RFinanceYJ}パッケージを使って、日本の株式市場からサンプルデータを引っ張ってきます。ひとまず適当に5銘柄の株価時系列を取ってきて、これにある変数xを回帰させます。この変数xは、過去の値はもちろん、未来の値も欲しいだけ得られるものとします。すると。。。

> st1<-quoteStockTsData('4751.t',since='2012-08-01',date.end='2013-07-31')
> st2<-quoteStockTsData('2432.t',since='2012-08-01',date.end='2013-07-31')
> st3<-quoteStockTsData('3632.t',since='2012-08-01',date.end='2013-07-31')
> st4<-quoteStockTsData('4755.t',since='2012-08-01',date.end='2013-07-31')
> st5<-quoteStockTsData('3793.t',since='2012-08-01',date.end='2013-07-31')
> x<-********** # ここはまだ秘密

# 中略：適当にこれらを組み合わせてデータフレームyを作る

> names(y)<-c("x","stock1","stock2","stock3","stock4","stock5")
# 株価系列の名前をstock1-5とし、説明変数をxとする

f:id:TJO:20130816101701p:plain

こんな感じの株価時系列です。では、これに変数xを回帰させてみましょう。

> summary(lm(stock1~x,y))

Call:
lm(formula = stock1 ~ x, data = y)

Residuals:
   Min     1Q Median     3Q    Max 
-25104 -11186  -2742   9261  53998 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 169746.3     1179.3  143.94   <2e-16 ***
x            -1880.4      135.8  -13.84   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 15590 on 244 degrees of freedom
Multiple R-squared: 0.4399,	Adjusted R-squared: 0.4376 
F-statistic: 191.7 on 1 and 244 DF,  p-value: < 2.2e-16 

> summary(lm(stock2~x,y))

Call:
lm(formula = stock2 ~ x, data = y)

Residuals:
    Min      1Q  Median      3Q     Max 
-817.33 -270.85   26.47  261.11  937.86 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2506.090     26.984   92.87   <2e-16 ***
x              4.381      3.108    1.41     0.16    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 356.7 on 244 degrees of freedom
Multiple R-squared: 0.008078,	Adjusted R-squared: 0.004013 
F-statistic: 1.987 on 1 and 244 DF,  p-value: 0.1599 

> summary(lm(stock3~x,y))

Call:
lm(formula = stock3 ~ x, data = y)

Residuals:
    Min      1Q  Median      3Q     Max 
-324.62 -112.63    2.60   98.43  310.92 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1322.460     11.765  112.41   <2e-16 ***
x             19.879      1.355   14.67   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 155.5 on 244 degrees of freedom
Multiple R-squared: 0.4687,	Adjusted R-squared: 0.4665 # ←すごい！！！
F-statistic: 215.2 on 1 and 244 DF,  p-value: < 2.2e-16 

> summary(lm(stock4~x,y))

Call:
lm(formula = stock4 ~ x, data = y)

Residuals:
    Min      1Q  Median      3Q     Max 
-283.01 -166.04    4.32  124.78  431.80 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  831.215     13.095   63.48   <2e-16 ***
x            -16.917      1.508  -11.22   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 173.1 on 244 degrees of freedom
Multiple R-squared: 0.3402,	Adjusted R-squared: 0.3375 
F-statistic: 125.8 on 1 and 244 DF,  p-value: < 2.2e-16 

> summary(lm(stock5~x,y))

Call:
lm(formula = stock5 ~ x, data = y)

Residuals:
   Min     1Q Median     3Q    Max 
-44418 -23377  -7390   9450 162519 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  76901.3     2747.0  27.994   <2e-16 ***
x            -2819.5      316.4  -8.911   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 36310 on 244 degrees of freedom
Multiple R-squared: 0.2455,	Adjusted R-squared: 0.2424 
F-statistic: 79.41 on 1 and 244 DF,  p-value: < 2.2e-16

5系列中4系列で有意な回帰関係が得られた上に、stock3に対しては何と調整済みR^2値が0.467にも達するという結果になりました。いやいやこれはすごい。こんな夢の変数xがあれば、いくらでも株価の予測ができますね。今すぐ大儲けできちゃいますよ！

ではその変数xの正体は何か？ということですが。実は上で伏字になっていたところはこうなっていました。

> x<-cumsum(rnorm(246)) # おい！ランダムウォークかよ！！！

そう、実はただのランダムウォーク。本来なら変数xは完全にランダムでデタラメな系列のはずなのに、これを実際の日本市場の株価系列に回帰させたらほとんどが有意な回帰を見せた挙句、かなり大きなR^2値までもが得られてしまったというわけです。これは一体どういうことなんでしょうか？？？

単位根過程同士を回帰すると、本来存在しないはずの回帰関係（見せかけの回帰）が勝手に生じてしまう

既に以前の記事でも紹介済みですが、再掲しておきましょう。要は、単位根過程の被説明変数に対して同じく単位根過程の説明変数によって単なるOLS（最小二乗法）による推定をかけると、偏回帰係数 $\beta_1, \beta_2, \cdots, \beta_n$ に関するt統計量が勝手に発散してしまい、その結果としてp値も爆発的に小さくなり「有意な回帰がある」という誤った結果になってしまう現象がある、ということです。

一応、沖本本pp.126-127の記述を引用しておきます。ブラウン運動の知識があると、前回の記事の時よりは理解しやすいと思います。

見せかけの回帰の現象を定義するために、2つの独立なランダムウォークを考えよう。具体的には、 $\epsilon_1t$ と $\epsilon_2t$ を独立なiid系列として、 $x_t$ と $y_t$ を
$\begin{array} x_t=x_{t-1}+\epsilon_1t, \epsilon_1t \sim iid(0,\sigma^2_1) \\ y_t=y_{t-1}+\epsilon_2t, \epsilon_2t \sim iid(0,\sigma^2_1) \end{array}$
で定義する。このとき、 $x_t$ と $y_t$ は独立なランダムウォークであるので、
$y_t=\alpha + \beta x_t + \epsilon_t$ (6.3)
となる回帰モデルを考えたとすると、真の値は $\beta = 0$ となる。したがって、(6.3)をOLSで推定し、 $\beta = 0$ という仮説検定を行うと、 $\beta = 0$ が採択される確率が高いはずである。しかしながら、変数の非定常性が驚くべき結果をもたらすことが知られている。具体的には、(6.3)をOLSで推定したときのOLS推定量 $\hat{\alpha}$ と $\hat{\beta}$ について、
$\left( \begin{array} ~ T^{-1/2} \hat{\alpha} \\ \hat{\beta} \end{array} \right) \stackrel{L}\longrightarrow \left( \begin{array} ~ \sigma_1 h_1 \\ (\sigma_1 / \sigma_2) h_2 \end{array} \right)$ (6.4)
が成立することが知られている*3。ここで $W_1(\cdot)$ と $W_2(\cdot)$ を独立な標準ブラウン運動として、
$\left( \begin{array} ~ h_1 \\ h_2 \end{array} \right)$
$= \left( \begin{array}{cc} 1 & \int W_2(r)dr \\ \int W_2(r)dr & \int W_2^2(r)dr \end{array} \right)^{-1} \left( \begin{array} \int W_1(r)dr \\ \int W_2(r)W_1(r)dr \end{array} \right)$
である。(6.4)の結果は、 $\hat{\alpha}$ が $\sqrt{T}$ の速度で発散し、 $\hat{\beta}$ がある確率変数に収束することを示している。これは、 $\sqrt{T}$ の速度で真の値に収束していく標準的な場合や、 $T$ の速度で真の値に収束していく単位根検定統計量の場合とは非常に異なるものであり、 $\hat{\alpha}$ も $\hat{\beta}$ も一致推定量ではないことがわかる。また、このため、 $\hat{\alpha}$ と $\hat{\beta}$ に関するt統計量は発散することになる。つまり、 $T$ が大きいとき、 $\hat{\alpha}$ と $\hat{\beta}$ を用いてt検定を行うと、ほぼ確実に $\hat{\alpha}=0$ と $\hat{\beta}=0$ という帰無仮説は棄却されるのである*4。また、もう1つの驚くべき事実として、回帰の決定係数 $R^2$ が漸近的に1に収束することも知られている。

この奇怪な現象は、最初にGrangerとNewboldが当時始まったばかりの計算機によるモンテカルロシミュレーションによって発見し(Granger & Newbold, 1974)、その後Phillipsが解析的な証明を与えたものです(Phillips, 1986)。Phillipsの論文は非常に難解ですが、興味のある人は是非読んでみてください。

ということでまとめると、沖本本p.127では見せかけの回帰を以下のように定義しています。

定義6.1（見せかけの回帰）　単位根過程 $y_t$ を定数と $y_t$ と関係のない単位根過程 $x_t$ に回帰すると、 $x_t$ と $y_t$ の間に有意な関係があり、回帰の説明力が高いように見える現象は見せかけの回帰(spurious regression)といわれる。

見せかけの回帰を回避するためには？

以上のような理由により、時系列データに対して回帰分析を用いる際には、何よりもまず変数が単位根過程かどうかに注意しなければいけないというわけです。ちなみに先ほどの株価時系列たちについてADF検定で調べてみると、

> adf.test(y$stock1)

	Augmented Dickey-Fuller Test

data:  y$stock1 
Dickey-Fuller = -2.9476, Lag order = 6, p-value = 0.177
alternative hypothesis: stationary 

> adf.test(y$stock2)

	Augmented Dickey-Fuller Test

data:  y$stock2 
Dickey-Fuller = -2.8982, Lag order = 6, p-value = 0.1978
alternative hypothesis: stationary 

> adf.test(y$stock3)

	Augmented Dickey-Fuller Test

data:  y$stock3 
Dickey-Fuller = -3.2026, Lag order = 6, p-value = 0.08802
alternative hypothesis: stationary 

> adf.test(y$stock4)

	Augmented Dickey-Fuller Test

data:  y$stock4 
Dickey-Fuller = -1.949, Lag order = 6, p-value = 0.5973
alternative hypothesis: stationary 

> adf.test(y$stock5)

	Augmented Dickey-Fuller Test

data:  y$stock5 
Dickey-Fuller = -2.8346, Lag order = 6, p-value = 0.2246
alternative hypothesis: stationary

案の定全て単位根過程という結果になりました*5。ということで、こういう時の解決策が沖本本pp.127-129に書かれています。かいつまんでまとめると、

そもそも線形回帰ではなくVARモデルで係数推定する
1階差分して定常過程に直してから回帰分析する

1番目は明快で、要は「時系列データなんだから（見せかけの回帰など時系列データ固有の問題にやられやすい）線形回帰なんかやらずにVARモデリングしろやゴルァ」ということです。この辺は単位根過程の説明でも触れた通りなので何も問題ないように見えますが、この場合F検定が無効になるのでGranger因果性検定も使えなくなるというデメリットがあります。

2番目はシンプルなアイデアで、1階差分して定常過程にしてしまえば何やっても問題ないでしょ！ということ。・・・ところがですね、沖本本p.128でも指摘されていますが差分系列というのは原系列に比べると色々と情報が抜け落ちてしまうのです*6。おまけに、この後で紹介する共和分関係がある場合には、そもそも差分系列を用いると誤った結論に達する可能性すらあります。

そこで、その共和分関係とは何ぞや？について次節で触れてみましょう。

共和分、ベクトル誤差修正モデル(VECM)

沖本本p.129に載っている詳しい定義を下に引用してありますが、要は共和分というのは「単位根過程の線形和が定常過程になる」関係のことです。

定義6.2（共和分）　 $x_t$ と $y_t$ を単位根( $I(1)$ )過程とする。 $ax_t+by_t$ が定常( $I(0)$ )過程となるようなaとbが存在するとき、 $x_t$ と $y_t$ との間には共和分(cointegration)の関係がある、もしくは $x_t$ と $y_t$ は共和分している(cointegrated)といわれる。より一般的には、 $\mathbf{y_t}$ を $I(1)$ とする。 $\mathbf{a'y_t}$ が $I(0)$ 過程となるような $\mathbf{a}$ が存在するとき、 $\mathbf{y_t}$ には共和分の関係がある、もしくは $\mathbf{y_t}$ は共和分しているといわれる。また、このとき、 $(a,b)'$ や $\mathbf{a}$ は共和分ベクトル(cointegrating vector)と呼ばれる。

特に $x_t$ と $y_t$ の間に共和分関係が存在する場合、 $z_t = y_t - a x_t$ が定常過程となるような係数aが存在することになるわけですが、これは式の形から見てもわかる通り $z_t$ がある値を中心に平均回帰的に振る舞うことが予想されます。つまり両者の間には一種の均衡関係が成立するわけで、共和分の有無を調べることには両者の時系列同士の関係性を推測する上でも、その未来予測をしていく上でも大きな意味があるとも言えるわけです。

以前の記事で紹介した例をちょっと見てみましょう。

> adf.test(x_ci)

	Augmented Dickey-Fuller Test

data:  x_ci 
Dickey-Fuller = -1.488, Lag order = 5, p-value = 0.7905
alternative hypothesis: stationary 

> adf.test(y_ci)

	Augmented Dickey-Fuller Test

data:  y_ci 
Dickey-Fuller = -1.3165, Lag order = 5, p-value = 0.8624
alternative hypothesis: stationary 

> summary(lm(y_ci~x_ci))

Call:
lm(formula = y_ci ~ x_ci)

Residuals:
   Min     1Q Median     3Q    Max 
-9.548 -6.178 -3.125  6.718 13.851 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -12.2700     0.9531 -12.874   <2e-16 ***
x_ci         -0.2409     0.1904  -1.265    0.207    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 7.393 on 199 degrees of freedom
Multiple R-squared: 0.007983,	Adjusted R-squared: 0.002998 
F-statistic: 1.601 on 1 and 199 DF,  p-value: 0.2072 

> matplot(cbind(x_ci,y_ci),lwd=3,lty=1,type="l")

f:id:TJO:20130423124145p:plain

プロットを見ると何となく分かるかもですが、2つの時系列の間には微妙に均衡関係（両者の平均値が一定のトレンドを描いている）が見られます。こういうケースでこそ、共和分関係が成立すると見て良いでしょう。

ところで、共和分に関するRでの分析プロセスを見る前に、一つ大切なことを書いておきます。その数理的基礎が沖本本pp.134-137で紹介されているGranger表現定理なんですが、これを全部引用していると大変なので以下のように内容をまとめると、

共和分関係にある2つの単位根過程 $y_{1t}, y_{2t}$ を仮定した場合（これを共和分システムと仮に呼ぶ）、この二者を表すVARモデルも単位根を持つが、共和分システムの場合その差分系列のVMAモデルは単位根を持つため、反転不可能となる*7。

反転不可能な(V)MAモデルは(V)ARモデルでは表現不可能であるため、仮に共和分システムの差分系列に対してVARモデルを推定したとしても、それは本来ならば成立し得ない誤ったモデルとなる。

平たく言えば、共和分システムに対してVARモデルで表現することは不可能だということです。そこで、これを何とかして修正してうまく表現するモデルがないか？ということで登場したのが、ベクトル誤差修正モデル(VECM: vector error correction model)です。

沖本本pp.136-137にその導出過程が全部書いてありますが、これまた手間なので割愛します。その代わりp.137の定理6.1がここまでの良いまとめとなっているので、以下に引用しておきます。

定理6.1（VAR(p)表現を持つ共和分システム $\mathbf{y}_t$ の性質）　VAR(p)表現を持つ共和分システム $\mathbf{y}_t$ は以下の性質をもつ。
（※ただしここで $\mathbf{y}_t = (y_{1t},y_{2t},\cdots,y_{kt})'$ とする。 $y_{kt}$ は個々の時系列を表す）

(1) VAR(p)表現は単位根を持つ。

(2) $\Delta \mathbf{y}_t$ のVMA表現は反転不可能である。

(3) $\Delta \mathbf{y}_t$ のVAR表現は存在しない。

(4) $\mathbf{y}_t$ は以下のVECM(p-1)で表現できる。

$\Delta \mathbf{y}_t = \zeta_1 \delta \mathbf{y}_{t-1} + \zeta_2 \delta \mathbf{y}_{t-2} + \cdots + \zeta_{p-1} \delta \mathbf{y}_{t-p+1} + \mathbf{\alpha} + \zeta_0 \delta \mathbf{y}_{t-1} + \epsilon_t$
$= \zeta_1 \delta \mathbf{y}_{t-1} + \zeta_2 \delta \mathbf{y}_{t-2} + \cdots + \zeta_{p-1} \delta \mathbf{y}_{t-p+1} + \mathbf{\alpha} - \mathbf{BA'} \mathbf{y}_{t-1} + \epsilon_t$
ここで、 $\mathbf{A, B}$ はn×h行列であり、 $\mathbf{A'y}_t \sim I(0)$ はh個の共和分関係を表す。また、共和分ランクhは行列 $\zeta_0$ のランクによって決まる。

(5) VECMには定常な変数しか含まれない。また、
$|\mathbf{I}_n - \zeta_1 z - \zeta_2 z^2 - \cdots - \zeta_{p-1} z^{p-1}| = 0$
のすべての解の絶対値は1より大きい。

この $-\mathbf{BA'}_{t-1}$ という項は誤差修正項(error correction term)と呼ばれ、均衡からの乖離が大きくなった時に均衡に戻っていく力が働くことを表しています。そしてVECMには定常な変数しか含まれていません。これにより、VECMは共和分システムの定常過程表現になるわけです。

ベクトル誤差修正モデル(VECM)と共和分関係の推定

沖本本pp.138-143に、代表的なJohansenの手順(Johansen's procedure)をメインにVECM及び共和分関係推定の詳細が載っています。・・・が、ここでその数理的基礎を全て紹介しているととんでもなく手間なので、これも記事中では割愛させてもらいます。

端的にまとめると、

最尤法で共和分ベクトルを表すVECMの中の $\mathbf{A}$ を推定したいので...
正準相関の概念に基づいて $\mathbf{A}$ を求める
最適な推定結果を選ぶために共和分の検定を行う必要があるが、その方法論としてトレース検定と最大固有値検定の2種がある

という流れです。この辺は導出過程を理解していると分析ツール等を使って実際のデータに適用する際にも役立つと思われるので、気になる方は沖本本を参照するようにしてみて下さい。

推定したVECMは、原理的には差分系列のVARモデルに誤差修正項を加えたモデルであり、ベースはあくまでもVARモデルです。なので、VARモデル同様にインパルス応答や分散分解を計算することもできます（Granger因果は無理）。ということで、共和分関係があるケースでもVARモデル同様の分析・考察を加えることができて有難い限りですね。

ところで、Rではこの辺の一連の計算はca.jo(){urca}関数でやってくれます。まず、共和分ランクの推定まで。

> data(finland)
> sjf<-finland
> sjf.vecm<-ca.jo(sjf,ecdet="none",type="eigen",K=2,spec="longrun",season=4)
> summary(sjf.vecm)

###################### 
# Johansen-Procedure # 
###################### 

Test type: maximal eigenvalue statistic (lambda max) , with linear trend 

Eigenvalues (lambda):
[1] 0.30932660 0.22599561 0.07308056 0.02946699

Values of teststatistic and critical values of test:

          test 10pct  5pct  1pct
r <= 3 |  3.11  6.50  8.18 11.65
r <= 2 |  7.89 12.91 14.90 19.19 # ←ココまで5%点が採択
r <= 1 | 26.64 18.90 21.07 25.75 # ←ココで5%点が棄却
r = 0  | 38.49 24.78 27.14 32.14

Eigenvectors, normalised to first column:
(These are the cointegration relations)

           lrm1.l2    lny.l2    lnmr.l2    difp.l2
lrm1.l2  1.0000000  1.000000  1.0000000   1.000000
lny.l2  -0.9763252 -1.323191 -0.9199865   1.608739
lnmr.l2 -7.0910749 -2.016033  0.2691516  -1.375342
difp.l2 -7.0191097 22.740851 -1.8223931 -15.686927

Weights W:
(This is the loading matrix)

           lrm1.l2       lny.l2      lnmr.l2      difp.l2
lrm1.d 0.033342108 -0.020280528 -0.129947614 -0.002561906
lny.d  0.022544782 -0.005717446  0.012949130 -0.006265406
lnmr.d 0.053505000  0.046876449 -0.007367715  0.002173242
difp.d 0.005554849 -0.017353903  0.014561151  0.001531004

共和分ランクはr = 2ということになりました*8。Rでは{urca}と{vars}が連携しているので、このVECMモデルをvec2var(){vars}関数を用いてVARモデル様の表現に改めることができます。

> sjf.vec2var<-vec2var(sjf.vecm,r=2)
> print(sjf.vec2var)

Coefficient matrix of lagged endogenous variables:

A1:
          lrm1.l1      lny.l1     lnmr.l1    difp.l1
lrm1  0.855185363 -0.28226832 -0.09298924 -0.1750511
lny   0.036993826  0.33057494 -0.06731145 -0.1946863
lnmr -0.156875074 -0.01067717  0.76861874  0.4247362
difp  0.001331951  0.02850137  0.02361709  0.2063468


A2:
         lrm1.l2      lny.l2     lnmr.l2     difp.l2
lrm1  0.15787622  0.27655060 -0.10255593 -0.52017728
lny  -0.02016649  0.65497929 -0.08102873 -0.09357761
lnmr  0.25725652 -0.10358761 -0.24253117  0.26571672
difp -0.01313100 -0.01096218 -0.02802090  0.36002057


Coefficient matrix of deterministic regressor(s).

        constant         sd1          sd2         sd3
lrm1  0.03454360 0.039660747  0.037177941  0.10095683
lny   0.05021877 0.043686282  0.082751819  0.09559270
lnmr  0.22729778 0.008791390  0.012456612  0.02011396
difp -0.03055891 0.001723883 -0.007525805 -0.00835411

VARモデルの時に{vars}でやったように、predict()関数で予測したり、irf()関数でインパルス応答を計算することももちろん可能です。

> sjf.pred<-predict(sjf.vec2var,n.ahead=25,ci=0.95)
# predict()関数が使える
> plot(sjf.pred)
> sjf.irf<-irf(sjf.vec2var,n.ahead=25,ci=0.95)
# irf()関数も使える
> plot(sjf.irf,plot.type="multiple",lwd=2)

f:id:TJO:20130815235757p:plain

f:id:TJO:20130815235806p:plain

という感じで、VARモデルと同じように可視化までできます。改めて単位根過程VARモデルに関するおさらいを書いておくと、

VARモデルを推定する前に、まず個々の原系列に対して単位根検定を実施する
そもそもただの回帰関係が見たいだけなら、気にせずVARモデルを計算して良い
単位根過程が混ざっていなければ、そのままVARモデルを計算して因果性分析まで行って良い
単位根過程が混ざっていたら、まず共和分関係（ランク）を推定する
共和分関係がなければ、差分系列に対してVARモデルを推定して因果性分析まで行って良い
共和分関係があったら、VECMを推定してからVARモデルに変換して初めて因果性分析を行える

これで大体問題はないかと思います。Rを使った場合でもかなり手間のかかる一連の手順ですが、ここまでやれば間違いはないかと。

ただし、実はインパルス応答の推定に限って言えば必ずしもこのような手順を踏む必要はなく、特に何もせずいきなりVARモデル推定→インパルス応答推定と突っ走ってしまっても問題ないかもしれない、という研究結果が報告されているとのこと（こちらの公開論文にその理由の解説があり、その論拠となる先行研究も引用されています）。インパルス応答は単位根過程が混じることによる影響を受けにくいという主張があるそうです。

・・・あー、長かった（汗）。これでVARモデル周りは大体網羅できた感じですね。

最後に

沖本本だとこの次は第7章「GARCHモデル」に当たるんですが、GARCHはwebマーケティングで使う機会に乏しい*9ので、今回のシリーズ記事ではスキップします。。。ということで次回は最後の第8章「状態変化を伴うモデル」の話をする予定です。

*1:定義式を参照のこと

*2:全部書いてたら本が書けてしまう。。。

*3:OLS推定量の表し方の問題なので、別に調べることをお薦めします

*4:t統計量が発散（爆発的増大）してしまい、勝手に回帰関係が統計的に有意という判定になってしまうため

*5:これを見ても分かる通り、実際の経済・ファイナンスデータにはものすごーーーく単位根過程は多いのです

*6:例えば値の大小といった情報は失われてしまう

*7:特性方程式の全ての解の絶対値が1より大きいという条件を満たさない（単位根過程なので絶対値が1に等しい解を持つ）

*8:r≦1では帰無仮説が棄却されてr≦2では帰無仮説が採択されているため

*9:金融だと超重要です、念のため

渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前