plmパッケージとpanelvarパッケージ
最近、パネルデータを扱うことが増えてきたのだが、パネルデータで動学的な（つまりt-1期とかのラグ項が出てくる）分析をやろうとすると最小二乗法ではなくGMM推定量を用いる必要がある。
備忘として、動学的パネルデータ分析（ダイナミックパネル分析）の基本的な考え方とRのパッケージの使いかたをここにメモしておこうと思う。といっても自分自身の理解もだいぶあやふやで、色々間違いもありそうなので、お読みになった方から指摘いただけると大変助かります…（汗）
Rの場合、結論から言うとまずは{plm}パッケージを使うのがいいと思う。いわゆる「パネルVAR」の形で分析したいなら、2018年に開発されたらしい{panelvar}パッケージを使うことができるのだが、後述するとおり{panelvar}のほうにはまだ不便なところもあって、今のところ、なるべく{plm}でできる分析をやったほうがいいと思う。まぁ、どちらも使い方自体は簡単である。なお、いずれにしても分析は線形モデルの範囲内になる。
{plm}パッケージについては日本語でもネット上に少しだけ解説があるが、たぶん{panelvar}のほうはまだ無いんじゃないかな？
なお理論的な背景については、幸い日本語でも、
千木良・早川・山本（2011）『動学的パネルデータ分析』
https://www.amazon.co.jp/dp/4862851029/
という非常に読みやすく実務的な利用にも配慮された教科書がある（神）。
また、クロスセクションと時系列の基本が頭に入ってる人なら、この奥井亮氏のスライドを読めばだいたいのストーリーは分かる（神）。
奥井亮：動学的パネルデータモデル
また、パッケージを使うときに、引数の設定や結果の解釈をするためにはヘルプ（plm・panelvar）を確認する必要があり、さらに意味が分からないところは作者の解説文書（plm・panelvar）もみることになると思われる。
しかし、とりあえずはさっきの奥井氏のスライドが大まかに理解できたら、あとは手を動かしながら並行して細かい解説を見ていったほうが早いのではないだろうか。
OLS（最小二乗法）の利用を妨げる諸問題
まず、社会調査データを扱う場合に発生する、OLSでは対処できない問題の種類、原因、対処法を大まかに知っておく必要がある。解説はネットで調べても山のように見つかるのでイメージだけまとめておく。
OLSで問題が生じる、ありがちなパターンとして、
- 目的変数yの誤差項と説明変数xが相関している（内生性）
- 目的変数yの誤差項どうしで相関（自己回帰の関係）がある（系列相関）
- 目的変数yの誤差項の分散が一定じゃない（不均一分散）
といったケースがある。
内生性があるとOLSではパラメータの一致推定量が得られない*1。系列相関は、よほど大きくなければ一致推定は得られるとされるが、パラメータの標準誤差が正確に推定できず仮説検定（パラメータが有意かどうかの判断）が有効でなくなる。不均一分散の場合も、基本的に一致推定は得られるものの、仮説検定が有効でなくなる。
内生性は、「欠落変数」（交絡変数）や、「同時性」（逆の因果）や、セレクションバイアスや、測定誤差等によって引き起こされるもので、とくに、欠落変数による見せかけの相関や、同時性による逆の因果については、査読や研究発表などでも指摘されることは多いだろう。内生性の問題には、基本的には、操作変数（yの誤差項とは相関しないがxとは相関する変数）を用いた方法で対処することになる。有名なのは二段階最小二乗法（2SLS）*2で、それを一般化したものと言えるのがGMM（一般化モーメント法、一般化積率法）*3。とりあえず、「操作変数ってスゴい」ということを理屈として実感するには、まず二段階最小二乗法を理解するのがいいと思う。
系列相関は、欠落変数により引き起こされる場合もあるし、自己相関的な状態依存（惰性や反動など）によって生じる場合もある。分析後に検定して系列相関がみられる場合、可能ならモデルを変えたほうがいいとは思うが、誤差項の分散共分散構造を理論的に仮定できるのであれば、GLS（一般化最小二乗法）で対処する場合もある*4。パラメータの標準誤差が正確に推定できない問題については、ロバスト標準誤差*5を用いて対処する。
不均一分散は、たとえば、規模がかなり異なる会社や国のデータが混じっている場合や、あるいは急速に成長したりしている場合に、元の数字が大きければ分散も大きくなるといった具合で発生する。その場合、変数変換して、たとえば人口やGDPや売上規模で割るとか、あるいはべき乗分布している変数の対数を取ったりすることで、分散が均一に近づく場合がある*6。ただ、変数変換により情報が失われる面もあるから、そもそもモデルの関数型を変えるべきではないかと言われることもある。不均一の原因（たとえば人口規模の違いとか）が分かっているのであれば、その原因変数をGLSの一種である重み付き最小二乗法（WLS）のウェイトに使うことでも対処できる。標準誤差は、ロバスト標準誤差を使う必要がある。
パネルデータにおけるバイアスへの対処
ところで、パネルデータを用いると、欠落変数については、それが「個体によって異なるが時間に対して不変（iに依存するがtに依存しない）な変数」や「時間によって変化するが個体間を通じて共通（tに依存するがiに依存しない）な変数」であれば、固定効果として考慮することで消去できる。
固定効果の考慮とは、具体的には、個体ダミー・年次ダミーを入れたり、階差を取ったり、元データに対して平均除去（demeaning、縦断方向と横断方向の平均をそれぞれ引く。中心化 centeringともいう。）を行ったりした上でOLS推定すればよい。ただしパラメータの標準誤差についてはロバスト標準誤差を用いた上で、仮説検定をする必要がある。
ちなみに、時間と個体の両方に依存するなど、複雑な動きをする欠落変数の影響はこの方法では消去できないことを忘れてはいけない。たとえば、100地域×10期分のデータがあり、目的変数と説明変数が１つずつあって双方に天候が影響していると思われるのだが天候のデータが手元になく、かつ地域と期のどちらでみても天候が似るような観測単位ではない場合は、困ったことになる。*7
なお、固定効果モデルのほかにランダム効果モデルというのもある。あまりポピュラーではない気がするがこの記事で説明されていて、「内生性が考えられるが系列相関はないと言える」場合なら固定効果モデル、「内生性はないが系列相関や分散不均一性はあるかも知れない」場合ならランダム効果モデル*8を用いるとされている。
さて、このように、欠落変数バイアスにある程度対処できるというのが、パネルデータを用いる大きな利点であると言える。・・・が、しかしである。
上の奥井氏のスライドのp.37〜38あたりの説明が分かりやすいが、パネルデータモデルにラグ変数を導入して動学的な関係を組み込んだ場合、固定効果モデルと同様の方法では、固定効果を除くときに、別の内生性が新たに発生してしまう！！
たとえば一階階差を取る場合でいうと、ΔYtの誤差項に含まれているYt-1の誤差が、説明変数側のΔYt-1にも含まれているからである。
そこで、たとえばラグを1期分まで含むモデルなら、その変数からみて1期以上前の値を操作変数に用いることで、内生性に対処することになる。奥井氏のスライドのp.37でいうと、(y_it−1 − y_it−2)とε_it−1のあいだに相関があって困るという話なのだが、yの2期ラグであるy_it−2なら、(y_it−1 − yit−2)とは相関をもつけどε_it−1とは相関をもたないと考えられるから*9、これを操作変数に使うわけである。
で、具体的には回帰係数をGMM推定量として求めるのだが、細かいことは上記の千木良・早川・山本（2011）などに書いてあるからいいとして、イメージを言えば、目的変数と説明変数と操作変数の関係をいくつかのモーメント条件として記述した上で、ある最適化問題に帰着させてこれを解くと、いい感じにパラメータの一致推定量が得られるというわけである（雑過ぎ？）。
いい感じと言っても、注意点はある。
- GMMは操作変数を大量に使う（使える）方法であり、「過剰識別制約」の検定（Sargenの検定やHansenのJ検定）を行って、要するに無駄に操作変数を取り入れ過ぎていないか（操作変数と撹乱項の直交が確保されてるか）を確認する必要がある。*10
- 千木良・早川・山本（2011）でも注意喚起されているが、系列相関が無いことが前提になっているので、推定後に一応、系列相関がないかどうかArellano-Bond検定（{plm}にはあるけど{panelvar}では提供されていない）等でチェックする。
GMMの種類
千木良・早川・山本（2011）ではGMMのバリエーションについて説明されている。まとめると・・・
一階階差GMM：最も一般的。後述のシステムGMMと違って、初期条件に仮定が必要ないところが最大の長所。短所として、撹乱項に系列相関があると一致性がなくなることが挙げられるが、系列相関は事後的に検定できるし、操作変数の工夫で回避できる（どういう工夫なのかは説明されてない）。それより、一番深刻な問題は、「弱操作変数」の問題である（弱操作変数の問題に対して頑健なのは後述のシステムGMM）。なお、1-step GMMと2-step GMMがあって、2ステップというのは1ステップ目の推定時の残差を用いて最適なウェイト行列を得てからもう1回推定するというやつなのだが、ふつうは2ステップのほうを使っておけばよい。*11
FOD-GMM：変数の変換の仕方が一階階差（first-difference）ではなくFOD（forward orthogonal deviation）になるだけで、利用可能な操作変数を全部使う場合は一階階差GMMと同じ推定値が得られる。操作変数を絞る場合は、少し結果に違いが生まれる。長所として、T≦N（クロスセクション方向の長さが時系列方向の長さを下回らない）であれば、FOD-GMMのほうが一階階差GMMよりも効率的な（分散が小さい）推定値を得られるらしい。
レベルGMM：モデルの説明変数から個別効果を取り除くのではなく、操作変数のほうから個別効果を取り除く方法。短所として、初期条件が平均定常でなければならない。しかし、レベルGMMはα（自己相関項の係数）が1に近づいても弱操作変数の問題が発生しにくいという長所がある。
システムGMM：一階階差GMMとレベルGMMをあわせて一つのシステムとみなして推定する（なお2ステップ推定が必要）。長所としては、αが1に近い場合でも弱操作変数の問題が生じないこと（αが1に近づくとレベルGMMのウェイトを増やすことで弱操作変数の問題を回避）。もうひとつは、一階階差GMMやレベルGMMよりも効率的であること。短所としては、レベルGMMのモーメント条件を含んでいることから、初期条件が平均定常でなければならない点。また、小標本ではバイアスが大きくなる可能性も指摘されている。
おそらく一階階差GMMかシステムGMMを使うのが基本なのだろうが、奥井氏の説明では、どういう場合にどっちが良いのかはいろいろ研究が進められてる段階で、議論が尽きないらしい。
{plm}パッケージの使い方
以上が前置きで、ここから具体的にRでの計算を実行していく。まずは{plm}からだ。
このパッケージの使いかたについては、以下の日本語資料で実行例が紹介されているので読むといいと思う。（ちなみにplmは「動学的」パネルデータ分析に限らず、静学的な固定効果モデルとかにも使えるもので、前半はその説明になっている。）
長倉大輔：R によるパネルデータモデルの推定
{plm}のヘルプの、pgmm()関数のところに載っている例を実行してみる。データは、{plm}に入っている練習用データのEmplUKを使う。
> data(EmplUK) > head(EmplUK, 20) firm year sector emp wage capital output 1 1 1977 7 5.041 13.1516 0.5894 95.7072 2 1 1978 7 5.600 12.3018 0.6318 97.3569 3 1 1979 7 5.015 12.8395 0.6771 99.6083 4 1 1980 7 4.715 13.8039 0.6171 100.5501 5 1 1981 7 4.093 14.2897 0.5076 99.5581 6 1 1982 7 3.166 14.8681 0.4229 98.6151 7 1 1983 7 2.936 13.7784 0.3920 100.0301 8 2 1977 7 71.319 14.7909 16.9363 95.7072 9 2 1978 7 70.643 14.1036 17.2422 97.3569 10 2 1979 7 70.918 14.9534 17.5413 99.6083 11 2 1980 7 72.031 15.4910 17.6574 100.5501 12 2 1981 7 73.689 16.1969 16.7133 99.5581 13 2 1982 7 72.419 16.1314 16.2469 98.6151 14 2 1983 7 68.518 16.3051 17.3696 100.0301 15 3 1977 7 19.156 22.6920 7.0975 95.7072 16 3 1978 7 19.440 20.6938 6.9469 97.3569 17 3 1979 7 19.900 21.2048 6.8565 99.6083 18 3 1980 7 20.240 22.1970 6.6547 100.5501 19 3 1981 7 19.570 24.8714 6.2136 99.5581 20 3 1982 7 18.125 24.8447 5.7146 98.6151
これは、ArellanoとBondがGMMによる動学パネルの分析方法を提案した1990年の論文で使っているデータで、firmが企業を表しているのだが、企業番号1〜140、年次は1977年から1983年までのデータがあるというものだ。
データフレームのレイアウトについて、{plm}パッケージでは個体やグループを表す変数は数字にしておく必要がある点に注意が必要だ。で、基本的には個体変数を1列目、時間変数を2列目に書くのだが、列が別の場所にあったとしても後述のindex引数で明示的に指定すればOK。
以下では、Arellano & Bond(1990)のTABLE4の(c)列と同じ結果が得られる分析をしてみよう。元論文では変数の対数を取っていて、上の長倉氏の実行例のようにformulaの中でlog()を挟めばいいのだが、コードが見づらくなるので先に全部対数を取っておく。
# 先に対数をとっておく for(i in 4:7){ EmplUK[,i] <- log(EmplUK[,i]) }
そしてパネルモデルの推定を実行！
> # パネルGMMを実行 > # Arellano & Bond(1990)のTABLE4の(c)列と同じ結果を出す > > z1 <- pgmm(emp ~ lag(emp, 1:2) # 自己相関項（ラグ1〜2） + + lag(wage, 0:1) # ラグ付き説明変数 + + capital # ラグ無し説明変数 + + lag(output, 0:1) # ラグ付き説明変数 + | lag(emp, 2:99), # 操作変数（最大限使う） + data = EmplUK, # データフレームを指定 + index=c('firm','year'), # 個体と時間変数の指定 + effect = "twoways", # 固定効果の設定 + model = "twosteps", # 2-stepGMMならtwostepとする + transformation = 'd' # 'ld'ならシステムGMM、'd'なら一階階差GMM + )
ラグを含める変数はlag(変数名,期）という書式で表現する。
期のところを0にするとラグではなく当期の値になる。つまり、0:2なら当期から2期前まで、1:2なら1期前から2期前までの変数がモデルに入るという意味である。もちろん「0」とだけ書いたらラグなし変数になる。
注意点として、{dplyr}パッケージの中にlagという関数があり、{dplyr}を読み込んでると変な誤作動が起きて「`n` must be a nonnegative integer scalar, not an integer vector of length...」みたいなメッセージが出たりする。（参考）
なので、dplyrの利用を避けるか、前処理で使ってしまっている場合は（というかパネルデータの前処理でdplyrのgroup_byとかが便利なので使うことが多いと思う*12）、plmを用いる直前で下記のように封印しておかなければならない。
detach(package:dplyr)
で、モデルの式を書いたあとに|を挟んで操作変数を定義する。これも書き方は同じで、上記ではlag(emp, 2:99)と書いているので、emp（雇用）の2期目から99期目までを操作変数に使うと指定していることになる。
ただここの99には意味がなく、これは要するに「使える操作変数は全部使う」みたいな意味合いである（もちろんTが100期以上あるならもっと大きくする。使えない期数のものは自動的に無視される）。ヘルプの書き方がそうなっているので、それにあわせている。操作変数の数を減らすときは、明示的に少ない数を与える。
操作変数は、右辺にある内生変数がラグ1期なら、ラグ2期以前の値を使う必要がある。内生変数がラグ1〜3期まであるなら、少なくとも操作変数としてはラグ２期〜4期（双方の個数が一致する時「丁度識別」という）のものは必ず含めておく必要がある。先程も言ったように、操作変数の数をどこまで入れるかというのは難しい問題で、分析例をみてると使える最大期数まで使っているものは多いのだが、後述する過剰識別制約の検定結果をみてアウトだったら減らす必要がある。
ところで、empの自己相関項以外の説明変数が上記では全部外生変数として考慮されてるわけだが、それらに内生性があることが懸念される場合は、それらに対する操作変数も与えればよい。具体的には、操作変数のところを、
| lag(emp, 2:99) + lag(wage, 2:99)
とかにすればよく、ヘルプにもそういう使用例が載っている。*13
indexは、何も指定しなければ、データフレームの1列目が個体変数、2列目が時間変数と解釈される。明示的に指定するときは、列名称又は列番号のベクトルで、c(個体変数,時間変数)の順に与える。
effectのところは、要するに時間効果（タイムダミー）を入れるかどうかの設定である。個体効果（ここでは企業=firmの個別効果）はもともと入っててそれが一階階差を取るときに消えるという話なんだが、ここを'twoways'にすると、時間効果もタイムダミーの形で考慮されるようになる。
modelのところは2-step GMMでやるのが普通なので、'twosteps'としておけばよいだろう。
transformationは、'd'にしておくと単なる一階階差GMM、'ld'にするとシステムGMMになる。
以下、推定結果。
> summary(z1, robust = F, time.dummies = F) Twoways effects Two steps model Call: pgmm(formula = emp ~ lag(emp, 1:2) + lag(wage, 0:1) + capital + lag(output, 0:1) | lag(emp, 2:99), data = EmplUK, effect = "twoways", model = "twosteps", transformation = "d") Unbalanced Panel: n = 140, T = 7-9, N = 1031 Number of Observations Used: 611 Residuals: Min. 1st Qu. Median Mean 3rd Qu. Max. -0.6190677 -0.0255683 0.0000000 -0.0001339 0.0332013 0.6410272 Coefficients: Estimate Std. Error z-value Pr(>|z|) lag(emp, 1:2)1 0.474151 0.085303 5.5584 0.0000000272221654 *** lag(emp, 1:2)2 -0.052967 0.027284 -1.9413 0.0522200 . lag(wage, 0:1)0 -0.513205 0.049345 -10.4003 < 0.00000000000000022 *** lag(wage, 0:1)1 0.224640 0.080063 2.8058 0.0050192 ** capital 0.292723 0.039463 7.4177 0.0000000000001191 *** lag(output, 0:1)0 0.609775 0.108524 5.6188 0.0000000192269973 *** lag(output, 0:1)1 -0.446373 0.124815 -3.5763 0.0003485 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Sargan test: chisq(25) = 30.11247 (p-value = 0.22011) Autocorrelation test (1): normal = -2.427829 (p-value = 0.01519) Autocorrelation test (2): normal = -0.3325401 (p-value = 0.73948) Wald test for coefficients: chisq(7) = 371.9877 (p-value = < 0.000000000000000222) Wald test for time dummies: chisq(6) = 26.9045 (p-value = 0.0001509)
coefficientsのところは分かりやすいと思うが、他に注目するところとしては、下の各種検定結果だ。
とくに、Sargen testのところは過剰識別制約の検定になっており、ここのp-valueが有意（例えば<.05）だと、操作変数入れすぎ。
またAutocorrelation testは系列相関のテストである。一階階差を取っていることから、誤差項に1次の自己相関があるのは自然なことなので、(1)のほうはp-valueが有意でも問題ない。(2)が有意（例えばp<.05）だと、系列相関があることになるので、その場合はモデルを変えたりする必要がある。
ところで、上記の論文ではロバスト標準誤差を使ってなくて、かつタイムダミーの係数の推定結果を載せてないのだが、それらを使う場合はsummary()の中で指定すればいいで。上のコードがそうなってるが、ここでタイムダミーを表示するかどうかも設定できる（表示するかどうかが違うだけで、計算自体は、effect = "twoways"としてあるならタイムダミーありで推定される）。
なおこれらは、pgmm()の引数としても設定できるのだが、summary()と矛盾する場合にsummary()側の設定が優先されるみたいなので、summary()で指定することにしとけばいいと思う。
> # ロバスト標準誤差を用い、かつタイムダミーのパラメータも表示 > summary(z1, robust = T, time.dummies = T) Twoways effects Two steps model Call: pgmm(formula = emp ~ lag(emp, 1:2) + lag(wage, 0:1) + capital + lag(output, 0:1) | lag(emp, 2:99), data = EmplUK, effect = "twoways", model = "twosteps", transformation = "d") Unbalanced Panel: n = 140, T = 7-9, N = 1031 Number of Observations Used: 611 Residuals: Min. 1st Qu. Median Mean 3rd Qu. Max. -0.6190677 -0.0255683 0.0000000 -0.0001339 0.0332013 0.6410272 Coefficients: Estimate Std. Error z-value Pr(>|z|) lag(emp, 1:2)1 0.4741506 0.1853985 2.5575 0.0105437 * lag(emp, 1:2)2 -0.0529675 0.0517491 -1.0235 0.3060506 lag(wage, 0:1)0 -0.5132048 0.1455653 -3.5256 0.0004225 *** lag(wage, 0:1)1 0.2246398 0.1419495 1.5825 0.1135279 capital 0.2927231 0.0626271 4.6741 0.000002953 *** lag(output, 0:1)0 0.6097748 0.1562625 3.9022 0.000095304 *** lag(output, 0:1)1 -0.4463726 0.2173020 -2.0542 0.0399605 * 1979 0.0105090 0.0099019 1.0613 0.2885484 1980 0.0246512 0.0157698 1.5632 0.1180087 1981 -0.0158019 0.0267313 -0.5911 0.5544275 1982 -0.0374420 0.0299934 -1.2483 0.2119056 1983 -0.0392888 0.0346649 -1.1334 0.2570509 1984 -0.0495094 0.0348578 -1.4203 0.1555141 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Sargan test: chisq(25) = 30.11247 (p-value = 0.22011) Autocorrelation test (1): normal = -1.53845 (p-value = 0.12394) Autocorrelation test (2): normal = -0.2796829 (p-value = 0.77972) Wald test for coefficients: chisq(7) = 142.0353 (p-value = < 0.000000000000000222) Wald test for time dummies: chisq(6) = 16.97046 (p-value = 0.0093924) >
{panelvar}パッケージの使い方
さて、こんどは{panelvar}パッケージを使ってみよう。
これはその名の通り、パネルデータでVARモデルを組むのに使えるもので、ふつうの時系列データのVARモデルのように、インパルス応答関数を出す機能も提供されている。
Stataのxtabond2をRに移植する目的で開発されたものらしいので、あとで暇ができたら、Stataの出力とも比較しておこうかなと。
さて、上記の{plm}と同じデータをつかって推定してみよう。変数を全部内生変数にして、4変数のVARにしてみる。
pvar.1<- pvargmm(dependent_vars = c('emp','wage','capital','output'), lags = 2, # 内生変数のラグ次数 transformation = "fod", # fdなら一階階差、fodならFOD変換 #predet_vars = c('',''), # 先決内生変数 #exog_vars = c('',''), # 外生変数 data = EmplUK, # 分析するデータフレーム panel_identifier=c('firm','year'), # 個体＆時間変数の指定 steps = c("twostep"), # 1-stepか2-stepか system_instruments = T, # TRUEだとシステムGMM max_instr_dependent_vars = 3, # 内生変数の操作変数の最大ラグ max_instr_predet_vars = 3, # 先決変数の操作変数の最大ラグ min_instr_dependent_vars = 2L, # 内生変数の操作変数の最小ラグ min_instr_predet_vars = 2L # 内生変数の操作変数の最小ラグ ) summ.1 <- summary(pvar.1)
dependent_varsのところには、内生変数を与える。列名称か、列番号。ここにあたえた変数でVARモデルが構成される。lagは内生変数のラグ次数。
transformationが'fd'なら一階階差変換でのGMM、'fod'ならFOD変換をしたGMMになる。
上の例では空っぽにしてあるのだが、predet_varsのところには先決変数を指定できる。先決変数というのは、文字通り先に決まるという意味で、yの誤差項と相関はしてないと考えられるが、yの誤差項が将来のxとは相関すると考えられるようなxのこと。その場合、一階階差GMMの処理でラグを取ったときに内生性が生まれてしまうのでこれらの変数にも操作変数を与える必要があるわけである。
exog_varsはそのまんま、外生変数を与える。つまり外生変数付きVARにできるわけだ。
panel_identifierは、{plm}のindexと同じ。
stepsは'twostep'で良いだろう。（plmと違って単数形で書くみたい）
system_instrumentsをT（TRUE）にしてるとシステムGMMになる。Fにすれば単なる一階階差GMMとかになる。
maxとminがついているところは、操作変数を何期前から何期前まで使うかを指定する引数である。
minを2にしてmaxを4にしていると、t-2からt-4までの値が操作変数に使われる。
ここもさっきと同様、ヘルプの使用例では最大限使うという意味でmaxは99にされていた。
minを「2L」というふうにLをつけてるのは、明示的に整数値を与えるという意味なのだが、なんでそうする必要あるのかはわからない。ただ、ヘルプでそうなってたのでそうしている。
推定結果をみてみる。
> summary(pvar.1) --------------------------------------------------- Dynamic Panel VAR estimation, two-step GMM --------------------------------------------------- Transformation: Forward orthogonal deviations Group variable: firm Time variable: year Number of observations = 611 Number of groups = 140 Obs per group: min = 4 avg = 4.364286 max = 6 Number of instruments = 308 ================================================================ emp wage capital output ---------------------------------------------------------------- lag1_emp 1.0963 *** 0.0690 0.1847 -0.0245 (0.0963) (0.0707) (0.1019) (0.0326) lag1_wage 0.1564 0.7162 *** -0.1312 -0.0792 (0.1431) (0.0873) (0.1271) (0.0478) lag1_capital 0.1722 ** -0.0056 1.1135 *** 0.0206 (0.0563) (0.0375) (0.0508) (0.0214) lag1_output 0.5364 *** -0.6185 *** 0.8483 *** 1.1223 *** (0.1191) (0.0932) (0.1548) (0.0550) lag2_emp -0.0960 -0.0434 -0.1164 0.0437 (0.0721) (0.0538) (0.0906) (0.0330) lag2_wage -0.0597 0.0874 0.1626 0.1107 ** (0.1013) (0.0640) (0.1013) (0.0367) lag2_capital -0.1623 *** -0.0153 -0.1574 ** -0.0329 (0.0477) (0.0369) (0.0500) (0.0307) lag2_output -0.6652 *** 0.4143 *** -0.8420 *** -0.5837 *** (0.1564) (0.1232) (0.1860) (0.0717) const 0.2682 1.5356 ** -0.2375 2.0043 *** (0.7296) (0.5465) (0.6384) (0.2880) ================================================================ *** p < 0.001, ** p < 0.01, * p < 0.05 --------------------------------------------------- Instruments for equation Standard GMM-type Dependent vars: L(2, 3) Collapse = FALSE --------------------------------------------------- Hansen test of overid. restrictions: chi2(272) = 132.37 Prob > chi2 = 1 (Robust, but weakened by many instruments.) >
一番下に、Hansenの過剰識別制約検定の結果が載っている。一番右の「= 1」というところが要するにp値で、これが有意でなければ過剰識別ではないのでOKということだが、1となってるので問題ないのだろう。
ちなみに、以下のようにそれだけを求める関数もある。
> hansen_j_test(pvar.1) $statistic [1] 132.3731 $p.value [1] 1 $parameter [1] 272 $nof_instruments [1] 308 $method [1] "Hansen-J-Test"
ところで、lagの次数を何期にすべきかの選びかたが難しい（{vars}パッケージのように選択するための機能が提供されていない）。推定後のpvarモデルをAndrews_Lu_MMSC()関数に与えるとAIC等が得られるので、ラグを変えながらちまちまAICを確認していけば良いのかも知れない。ちなみにpvargmmの1回の推定にかなり時間がかかるのでとても大変・・・。
> # HQは予めパラメータの設定がいるが2.1にしておけばいいと思う（dn=2loglogn） > Andrews_Lu_MMSC(pvar.1, HQ_criterion = 2.1) $MMSC_BIC [1] -1785.741 $MMSC_AIC [1] -465.6269 $MMSC_HQIC [1] -1034.676
インパルス応答関数は以下のようにして求められる。
# インパルス応答関数（結果はfrom=impulse側の変数ごとにまとめられている） irf.1 <- oirf(pvar.1, n.ahead=10) plot(irf.1)
なお、irf.1の中身は見づらいのだが、リストがimpulse側の変数ごとにまとめられていると思って見ていけば分かると思う。
インパルス応答関数の信頼区間も出すことができるのだが、これは極めて危険で、えげつない時間がかかる。
boot.1 <- bootstrap_irf(pvar.1, typeof_irf = c("OIRF"), n.ahead=10, nof_Nstar_draws=100, confidence.band = 0.95)
typeof_irfのところは、OIRFにしておくと直交化インパルス応答、GIRFにすると一般化インパルス応答関数になる。
nof_Nstar_drawsはブートストラップのサンプリングの回数。これ、{vars}パッケージで求める（パネルじゃないデータの）IRFの場合、デフォルトが確か100回になってて一瞬で終わるのだが、この{panelvar}のブートストラップは1回分やるのに私のパソコンだと数分かかる場合もあるので、とてもじゃないが気軽には回せない。
なお、内生変数を1個だけにして他を外生変数や先決変数にすればVARではなくなり、{plm}でやっているのと同じような分析になるのだが、計算の速さとかからして{plm}のほうがいいと思う。また、{panelvar}では、系列相関の検定結果が提供されないので、自分で理論を確認しながら組む必要があり、そんなことするぐらいならやはり素直に{plm}を使っておけばいいのではないだろうか。
もちろん、「パネルVAR」にしたいケースもあるとは思われ（ちゃんと読んでないが、ECBがパネルVARをつかった研究のサーベイをしてた）、その場合は{panelvar}を使うことになるけど、いろいろ調べても解説がまだ少なく、私のようなこの種の分析に慣れてない人間には使いづらそうではある。
パネル単位根検定
パネルデータも、パネルじゃない時系列の場合と同様に、単位根や共和分の問題に気をつける必要がある（といっても今のところ、どういう場合にどこまで問題になるのか自分でイマイチ頭が整理できていない）。
千木良・早川・山本（2011）では、80ページぐらいを費やして様々なタイプのパネル単位根検定・定常性検定が紹介されていて、その後にパネル共和分についても解説されている。
パネル単位根検定としてはシンプルかつ古いのがLevin, Lin & Chu(LLC)の検定で、これはADF検定をパネルに拡張したもの。これに加えて個体間の差異を認める方向で拡張したのがIm, Pesaran & Shin (IPS)の検定。これをさらに画期的な方法で拡張したと言われるのがcombination検定というやつで、Maddala & Wu (1999)とChoi (2001)によって別々に提案された。
以上のものはすべてADF検定の拡張であり、帰無仮説が単位根で対立仮説が定常という「単位根検定」である。逆に、帰無仮説が定常で対立仮説が単位根である「定常性検定」というのもあって、これは時系列分析におけるKPSS検定ってやつをHardi(2000)がパネルデータに拡張したものがある。
以上は「第１世代のパネル単位根検定」と呼ばれる手法群で、{plm}ではこれらの機能が提供されている。下記のpurtest()関数で、引数testに、LLC検定は'levinlin'、IPS検定は'ips'と設定する。combination検定については、Maddalaらの検定（これはChoiの提案したいくつかの手法のうちp検定と呼ばれるもの）は'madwu'と設定する。Choiが他にいくつか提案してる方法も実装されていて、Nが大きい（クロスセクション方向が長い）場合のp検定は'Pm'、あとはChoi(2001)を読んでないのでよくわからないがinverse normal testというのが'invnormal'、logit testというのが'logit'という設定で使える。Hardiの定常性検定はそのまんま'hardi'とする。
Hardiだけは毛色が違っていて、これは「すべての個体において定常」というのが帰無仮説なので、p値が有意なら「少なくとも1つが単位根」となる。他は、p値が有意であれば、単位根が棄却される。
> # パネル単位根検定 > purtest(object= EmplUK, + x = 'emp', + data = EmplUK, + index = c('firm','year'), # 個体、時間の順で入れる + test = 'madwu', # 'levinlin','ips','madwu','Pm','invnormal','logit'or'hadri' + exo = 'trend', # 'trend','none'or'intercept' + lags = 'AIC', # ラグはAICで適当に選んでという意味 + pmax = 6, # 考慮するラグの最大次数 + Hcons = TRUE # test="hardi"のときだけ関係ある設定 + ) Maddala-Wu Unit-Root Test (ex. var.: Individual Intercepts and Trend) data: EmplUK chisq = 597.8, df = 14, p-value < 2.2e-16 alternative hypothesis: stationarity
時系列分析のときと同じように、単位根検定にドリフト項（定数項）やトレンド項を入れるかどうかというのは大きな問題となる。上記では、exoのところを'none'にするとどっちも無し、'intercept'にするとドリフト項あり、'trend'にするとドリフト項もトレンド項もありとなる。
ところで大きな問題点として、パネル単位根検定は「すべてが単位根過程である」という帰無仮説を検定しているので、逆にいうと1つでも単位根仮定でない個体があるなら、「定常である」という対立仮説が採用されることになってしまう。これは、単位根の検出力が非常に低いことは意味していて、パネル単位根検定で単位根が棄却できるからといって、個々の個体の時系列が定常であるわけではない。しかし、千木良・早川・山本（2011）が説明してるのだが、パネルデータ分析というのはそもそも各グループ・各年に共通のパラメータを求めに行く分析なのだから、これは本質的に避けがたい問題であると考えたほうがよいようである。
以上はすべて「第一世代の単位根・定常性検定」と言われるもので、クロスセクション方向には基本的に独立性を仮定しているのだが、「第二世代」の検定手法群になるとクロスセクション方向の相関も認めるという拡張が行われる。{plm}パッケージでは第２世代のうちCross-sectionally augmented IPS testというのをcipstest()という関数で実行できる。帰無仮説は「単位根過程あり」。
> # pseriesというクラスのデータを与えないと行けないので… > Produc.p <- pdata.frame(Produc, index=c('state', 'year')) > > cipstest(Produc.p$gsp, + lags = 1, + type = 'trend', + model = 'cmg', + truncated = FALSE) Pesaran's CIPS test for unit roots data: Produc.p$gsp CIPS test = -0.70295, lag order = 1, p-value = 0.1 alternative hypothesis: Stationarity Warning message: In cipstest(Produc.p$gsp, lags = 1, type = "trend", model = "cmg", : p-value greater than printed p-value >
帰無仮説が棄却されてないので、単位根ありってことだろう。
その他
上の例で書かなかったが、固定効果が含まれているかどうかの検定をHausman検定を使って行うことができるがいまいちロジックを理解してないので後ほど確認しておきます。
千木良・早川・山本（2011）の共和分検定の章をまだ読んでないのでこれもあとで確認します。
また、（一階階差GMMの時に問題になる）初期条件の平均定常性の検定なども行うことが望ましいとされます。
パッケージの比較
現時点での、{plm}と{panelvar}の比較をまとめておきます。
- パネル単位根検定は、plmでのみ提供されている（第2世代までカバーされている）。
- Hausman検定も、plmでのみ提供されている。
- パネル「VAR」モデルの分析は、plmでは実行できずpanelvarが必要。
- 系列相関の検定が、plmでは提供されているが、panelvarでは提供されていない。
- 過剰識別制約の検定はどちらでもできる。
- モデル全体のAICを確認する手段が、panelvarでは提供されているが、plmでは提供されていない。（Stackoverflowでそれを試みているエントリはあったが詳細確認してません。）
- 計算速度は、plmは気になるほどの時間はかからないが、panelvarは遅い。
*1:内生性があるかどうかはまず理論的に考えるべきだが、結論を歪めるほどのものであるか否かは、HausmanやWooldridgeの方法による内生性検定で確認できる。 http://www.ier.hit-u.ac.jp/~kitamura/lecture/Hit/08Statsys5.pdf
*2:GLS（誤差の分散共分散行列に制約を仮定する「一般化最小二乗法」）を組み合わせた三段階最小二乗法（3SLS）というものもある。
*3:完全情最尤推定法というのもあるらしいがよく分かってない。
*4:コクラン＝オーカット法や、プレイス＝ウィンステン法。
*5:なんか種類は色々ある。
*6:対数差分を取るとそれは変化率に相当するが、経済データだと「前年比」に着目するとうまく捉えられる変数も多いので、理論的にも都合がいいとされる
*7:個別のケースごとになら、色々細かい対処法が考案されてそうな気はするが
*8:ランダム切片モデルというやつで、この場合、一種のGLSとみなすことができる。
*9:厳密に言えば、Yt-1の誤差項とYt-2が絶対に相関しないのか、つまり新しい期の値の誤差が古い期の値と相関を持たないのかというと、その変数の自己相関がうまく定式化されてなくて、古い期の値が新しい期の値に十分反映されてない場合、相関が残るとは思う。
*10:遠い過去の値はそもそも「弱操作変数」の問題を引き起こすので、千木良他（2011）では実務上は「3期まで」に限る方針がとられることが多いとしているが、どれぐらい一般的なのかは不明。
*11:ちなみに1ステップGMMはいわゆる二段階最小二乗法と同値。
*12:パネル用データフレームであるpdata.frameクラスの挙動を理解すればdplyrを使わなくても便利に計算できるのかも知れないが、確認していない。
*13:wageのほうを1:99にしたほうがいい場合もある気がするが、理論的にどういう内生性を考えてるかによるかな？ヘルプで引用されてる Blundell and Bond (1998)をあとで見て考えよう。
