はじめに

パネルデータ分析について興味があり少し調べてみました。

パネルデータ分析と動学的パネルデータについての基本的な部分を簡単にまとめました。

こちらの書籍も目を通してみました。パネルデータを扱うための基本的な注意点などがまとまってり、あまり数理的な部分には深入りしていないです。

パネルデータの調査と分析・入門

発売日: 2016/11/30
メディア: 単行本

パネルデータ

パネルデータとは複数の観測個体を複数の時点に渡って観測したデータ。

例えば、複数の人物に数年ごとに生活状況や労働環境について調査したデータや、野生動物をロギングし接触ごとに測定したデータなどがあげられる。

$\displaystyle{ y_{it} , i=1, \cdots, N, t=1,\cdots, T }$

$i$ ：観測個体 (個人、企業、国など) $t$ ：時点 (年、４半期など)

一方で、複数年のクロス・セクションデータを調査度年の区別することなくプールして使うデータはプールド・データと呼ばれている。

パネルデータ分析

パネルデータを用いた回帰分析

固定効果モデル

観測されていない個体特性を考慮しないと、観測されている変数の効果を大きく誤って推定してしまう場合がある。

そのため、次のようなモデルを検討する

$\displaystyle{ y_{it} = \beta x_{it} + \mu_i + \epsilon_{it} }$

$\mu_i$ は、時刻に依存じない個体特性を表現するダミー変数であり、LSDVと呼ばれている。

なお、観測時刻が少なく、個体数が多い場合にはLSDVの推定は難しくなる。

また、誤差項と説明変数の相関があるので最小二乗推定では、適切な推定ができない。

そのため、個人内の平均値を用いる次のモデルを考える。

$\displaystyle{ \bar y_{i} = \beta \bar x_{t} + \mu_i + \bar \epsilon_{t} }$

上記式から上記の式を引くと次のような式となり、個体特性を推定式から消すことができる。

$\displaystyle{ y_{it} - \bar y_{i} = \beta ( x_{it} - \bar x_{t} ) + (\epsilon_{it} - \bar \epsilon_{t}) }$

これにより観測期間中に変化する個人特性の影響を除くことで、OLSによる $\beta$ の推定を行うことができる。

このような推定はwithin(個体内)推定/固定効果(fixed effect)推定と呼ばれる。

変量効果モデル

個人特性を孤立変数とは独立にランダムに発生しているもとして扱う。

具体的には次のように個人特性を撹乱項の一つの要素として扱るモデルを推定することになる。

$\displaystyle{ u_{it} = \mu_i+\epsilon_{it} }$

$\displaystyle{ \bar y_{i} = \alpha + \beta \bar x_{t} + \mu_{it} }$

このモデルが推定できる条件は、撹乱項と独立変数が相関していないことである。

しかし、勉強時間を独立変数とする時に、勉強意欲などの個人特性と相関してしまうことが多い。

特に個体数が多く観測時点が少ないパネルデータではこの条件を満たすことが難しい。

また、変量効果モデルはOLSによって $\beta$ の推定を行えるが、 $t \neq s$ の時でも，同一の対象内では[tex: u { it } ]と[tex: u { is }]は相関しているので，通常のOLSの標準誤差では適切でない。

そのため、誤差項間の相関に頑健な標準誤差であるクラスターロバスト標準誤差を使用する必要がある。

クラスターロバスト標準誤差は不均一分散と自己相関を考慮した標準誤差であり、ある同じクラスター(同じグループ)内では相関があり，異なるクラスターの誤差項では無相関となる。

なお、固定効果推定でもクラスターロバスト標準誤差を使用するべきであるとされている。固定効果変換しても相関が完全に消えるわけではない。

固定効果モデルと変量効果モデルのどちらを選択するか

検定を用いて機械的に選択する方法は次のようなものがある。

F検定：固定効果モデルかプールドOLSのどちらが適切か選択固定効果モデルにおける固定特性ダミーがゼロか否かを検定し、個体特性ダミーをモデルに加える意味があるかを検証する。
Breusch-Paganのラグランジュ乗数検定：変量効果モデルかプールドOLSのどちらが適切か選択変量効果モデルにおける個体特性の分散が0か否かを検定し、個体特製をモデルに加える意味があるのかを検証する。
ハウスマン検定：固定効果モデルか変量効果モデルのどちらが適切かを選択個体特性が独立変数と相関しているかを検定する。具体的には、両方のモデルにおける時間変化する変数の推定された係数を比較し、差があるかを検定することで、両モデルに差があるかを見る。帰無仮説が棄却された場合は固定効果モデルを、採択された場合はどちらのモデルでも良いということになる。

観測サンプルの脱落について

パネルデータでは、間隔を開けて複数回観測するために、途中から観測ができなくなってしまう=未回答のサンプルが発生してしまうことが多い。このような未回答を単位未回答と呼ぶ。

また、回答が回収できたとしても未回答の部分があることも考えられる。このような未回答を項目未回答と呼ぶ。

これらの観測の脱落に対する対処法は主に、代入を用いる方法とウェイトを用いた補正があげられる。

代入法

いくつかの欠測値への代入法方が存在する。

回帰分析を用いて代入する方法
データセット内の背景情報が類似したデータを用いて補完するHot Deck法
確率的に発生させた複数の値を代入し総合的に評価する多重代入法

ウェイトを用いた補正

今までに得られている回答などを用いて、回答が得られるサンプルと未回答のサンプルを識別するロジスティック回帰やプロビット回帰など作成する。

得られたモデルによる予測確率(傾向スコア)の逆数を各サンプルかけることで、補正を行う。

この時、作成する識別モデルにおける共変量の選択は慎重に行う必要がある。

動学的パネルデータ分析

動学的パネルデータ分析は、パネルデータ分析で用いるモデルに自己回帰などの時系列分析の要素を加える分析方法である。

パネルAR(1)モデル

最も代表的なパネル AR(1) モデルは次のようなモデルに表現される。

$\displaystyle{ y_{it} = \alpha y_{i,t-1} + \mu_i + \epsilon_{it} }$

説明変数を含むモデルは次のように表現される。

$\displaystyle{ y_{it} = \alpha y_{i,t-1} + \beta x_{i,t} + \mu_i + \epsilon_{it} }$

一階差分モデル

個別特性 $\mu_i$ を取り除くために、時刻が１ズレてる情報と差分を取る。

$\displaystyle{ y_{it} - y_{i,t-1} = \alpha (y_{i,t-1} - y_{i,t-2} ) + ( \epsilon_{it} -\epsilon_{i,t-1} ) }$

このようなモデルを一階差分モデルと呼ぶ。

このモデルに対してOLSを適用することで、 $\alpha$ を推定する。

$\displaystyle{ \hat \alpha_{FD} = \frac{\sum_{i=1}^N \sum_{t=3}^T (y_{i,t-1} - y_{i,t-2} ) (y_{it} - y_{i,t-1} )}{\sum_{i=1}^N \sum_{t=3}^T (y_{i,t-1} - y_{i,t-2} )^2} }$

この推定量を一階差分（FD）推定量と呼ぶ。

この推定量は、 $N$ → $\inf$ で $T$ を固定したときには一致性をもたない。

操作変数法

パネルAR(1)モデルは次のように変形できる。

$\displaystyle{ y_{it} = \alpha y_{i,t-1} + \beta x_{i,t} + \mu_i + \epsilon_{it} }$

$\displaystyle{ y_{it} = \alpha ( \alpha y_{i,t-2} + \beta x_{i,t-1} + \mu_i + \epsilon_{i,t-1} ) + \beta x_{i,t} + \mu_i + \epsilon_{it} }$

この時、 $\alpha$ を通じて $y_{i,t-2}$ と $\mu_i$ が相関してしまう。

そのため、固定効果モデルを用いる。

$\displaystyle{ y_{it} - y_{i,t-1} = \alpha (y_{i,t-1} - y_{i,t-2} ) + ( \epsilon_{it} -\epsilon_{i,t-1} ) }$

推定される係数 $\hat \beta$ の確率極限は次のようになり、推定にバイアスがある。

$\displaystyle{ p lim( \hat \beta) = \beta + \frac{cor( (y_{i,t-1} - y_{i,t-2} ),(u_{it}-u_{i,t-1})) }{var(y_{i,t-1} - y_{i,t-2} )} }$

そのため、操作変数法を用いる。

2期前と3期前の $y$ の差分は、操作変数としての条件を満たしている。

2期前と3期前の $y$ の差分を操作変数として用いた操作変数法により推定される係数は次にようになる。

$\displaystyle{ \beta_{IV} = \frac{\sum_{i=1}^N \sum_{t=3}^T (y_{i,t-2} - y_{i,t-3} )(y_{i,t} - y_{i,t-1} )} { \sum_{i=1}^N \sum_{t=3}^T (y_{i,t-2} - y_{i,t-3} )(y_{i,t-1} - y_{i,t-2} ) } }$

他の説明変数を用いると推定値の分散が大きくなってしまうことが知られている。

操作変数を増やすことによって、分散を小さくすることができる。

2期前と3期前の $y$ の差分だけではなく、k期前の $y$ の差分も操作変数にすることができる。

しかし、パラメータが増えることになり、解析的な係数の推定が困難になる。

そのため、推定には一般化モーメント法(GMM)を用いる。GMMについては理解が追いついてないので割愛。

Rでやってみる

パネルデータ分析をRで行うplmパッケージを用いて行う。

今回用いるデータは、ガソリンに関する統計と所得に関する統計のデータを用いる。

カラムの意味は以下の通り。

　country：18ヶ国の国名

　year：年（1960～1978）

　lgaspcar：自動車1台あたりのガソリン消費量の自然対数

　lincomep：一人あたり実質所得の自然対数

　lrpmg：実質ガソリン価格の自然対数

　lcarpcap：一人あたりの自動車の量（ストック）の自然対数

library(tidyverse)
library(plm)

data(Gasoline)

データの中身を簡単に確認

  country year lgaspcar  lincomep      lrpmg  lcarpcap
1 AUSTRIA 1960 4.173244 -6.474277 -0.3345476 -9.766840
2 AUSTRIA 1961 4.100989 -6.426006 -0.3513276 -9.608622
3 AUSTRIA 1962 4.073177 -6.407308 -0.3795177 -9.457257
4 AUSTRIA 1963 4.059509 -6.370679 -0.4142514 -9.343155
5 AUSTRIA 1964 4.037689 -6.322247 -0.4453354 -9.237739
6 AUSTRIA 1965 4.033983 -6.294668 -0.4970607 -9.123903

データをplm()で扱えるようにpdata.frame型に変換する。 index引数で個体識別と時刻を表すカラムを指定する。

Gasoline_pdf <- Gasoline %>% 
  pdata.frame(index=c("country","year"), drop.index=TRUE)

> head(Gasoline_pdf)
             lgaspcar  lincomep      lrpmg  lcarpcap
AUSTRIA-1960 4.173244 -6.474277 -0.3345476 -9.766840
AUSTRIA-1961 4.100989 -6.426006 -0.3513276 -9.608622
AUSTRIA-1962 4.073177 -6.407308 -0.3795177 -9.457257
AUSTRIA-1963 4.059509 -6.370679 -0.4142514 -9.343155
AUSTRIA-1964 4.037689 -6.322247 -0.4453354 -9.237739
AUSTRIA-1965 4.033983 -6.294668 -0.4970607 -9.123903

このデータを用いて、所得やガソリン価格がガソリン消費量に影響を与えると仮定するモデルを生成する。

モデルの推定にはplm()関数を用いて行います。

重要な引数は次の２つです。

model引数でモデルを指定します。

　fd：差分モデル

　within：主体固定効果モデル

　between：時間固定効果モデル

　pooling：プールドOLS

effect引数で推定する固定効果を指定します

　individual：主体固定効果

　time：時間固定効果

　twoways：主体・時間の両方

今回は、主体固定効果モデル

plm_result <- plm(lgaspcar ~ lrpmg + lincomep + lcarpcap, data = Gasoline_pdf,
                                model = "within", effect = "twoways")

結果を見てみます。はじめにデータです。

> summary(plm_result)
Twoways effects Within Model

Call:
plm(formula = lgaspcar ~ lrpmg + lincomep + lcarpcap, data = Gasoline_pdf, 
    effect = "twoways", model = "within")

Balanced Panel: n = 18, T = 19, N = 342

Residuals:
       Min.     1st Qu.      Median     3rd Qu.        Max. 
-0.41920085 -0.03886111  0.00018502  0.04199566  0.23067839 

Coefficients:
          Estimate Std. Error  t-value  Pr(>|t|)    
lrpmg    -0.192850   0.042860  -4.4995 9.718e-06 ***
lincomep  0.051369   0.091386   0.5621    0.5745    
lcarpcap -0.593448   0.027669 -21.4479 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    10.644
Residual Sum of Squares: 1.997
R-Squared:      0.81239
Adj. R-Squared: 0.78886
F-statistic: 437.338 on 3 and 303 DF, p-value: < 2.22e-16
>

クラスターロバスト標準誤差とそれに基づく検定結果を確認する。

lmtest::coeftest(plm_result,
                 vcov = vcovHC(plm_result,
                               method = "arellano",
                               type = "HC1", 
                               cluster = c("group", "time")))

ガソリン価格や所得がガソリン消費量に与える影響は統計的に有意ではないと考えられる。

t test of coefficients:

          Estimate Std. Error t value  Pr(>|t|)    
lrpmg    -0.192850   0.124139 -1.5535    0.1213    
lincomep  0.051369   0.231977  0.2214    0.8249    
lcarpcap -0.593448   0.080677 -7.3559 1.787e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rで動学的パネルデータ分析

Rで動学的パネルデータ分析を行うには, pgmm関数を用います。

用いたデータはplmパッケージに含まれる、EmplUKというイギリスにおける企業の雇用活動に関するデータを用います。

カラムの意味は以下の通り。

　firm：各企業のインデックス

　year：年（1976～1984）

　sector：分野での経済活動量

　emp：雇用

　wage：実質賃金

　capital：総資本

　output：総生産高

> data("EmplUK")
> head(EmplUK)
  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

こちらのデータを用いて、次のような今年のある企業の雇用数は「賃金・資本量・業種自体への需要・経済全体でのトレンド」に影響を受けると仮定したモデルを作成します。

$\displaystyle{ log(emp) = \alpha_1 log(emp_{i,t-1})+\alpha_2 log(emp_{i,t-2}) + \beta_1 * capital_{i,t}+ \beta_2 * output_{i,t} + \\ \beta_3 * output_{i,t-1} + year_t + \lambda_i+u_{it} }$

目的変数を対数変換したのちにモデルの推定を行う。なお、これを実行する際にはモデル定義の中にあるlag()がdplyrのlag()と競合してモデル推定でエラーになってしまう。そのため、dplyrパッケージを読み込んでいない状態にする必要がある。

# 変数の対数変換
UK_panel_df <- dplyr::mutate(EmplUK, emp = log(emp), capital = log(capital),output=log(output))
  
# モデル式の定義
dynam_model <- emp ~ lag(emp, 1:2) +    # 自己相関項（ラグ1〜2）
                     lag(wage, 0:1) +   # ラグ付き説明変数
                     capital +          # ラグ無し説明変数
                     lag(output, 0:1)|  # ラグ付き説明変数
                     lag(emp, 2:99)     # 操作変数（最大限使う）

# モデル推定
est_gmm <- pgmm(formula = dynam_model,
                   data = UK_panel_df,
                 effect = "twoways",
                  model = "twosteps")

結果を見てみます。

> summary(est_gmm, robust = FALSE)
Twoways effects Two steps model

Call:
pgmm(formula = dynam_model, data = UK_panel_df, effect = "twoways", 
    model = "twosteps")

Unbalanced Panel: n = 140, T = 7-9, N = 1031

Number of Observations Used: 611

Residuals:
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.7579094 -0.0244755  0.0000000  0.0000899  0.0308044  0.6617895 

Coefficients:
                    Estimate Std. Error z-value  Pr(>|z|)    
lag(emp, 1:2)1     0.4174437  0.0797793  5.2325 1.672e-07 ***
lag(emp, 1:2)2    -0.0434140  0.0250933 -1.7301 0.0836119 .  
lag(wage, 0:1)0   -0.0210767  0.0022284 -9.4582 < 2.2e-16 ***
lag(wage, 0:1)1    0.0083009  0.0032005  2.5937 0.0094963 ** 
capital            0.2950395  0.0401729  7.3442 2.069e-13 ***
lag(output, 0:1)0  0.5956860  0.1065380  5.5913 2.254e-08 ***
lag(output, 0:1)1 -0.4113409  0.1160983 -3.5430 0.0003955 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sargan test: chisq(25) = 28.47986 (p-value = 0.28624)
Autocorrelation test (1): normal = -2.039692 (p-value = 0.041381)
Autocorrelation test (2): normal = -0.1806517 (p-value = 0.85664)
Wald test for coefficients: chisq(7) = 290.7633 (p-value = < 2.22e-16)
Wald test for time dummies: chisq(6) = 28.72852 (p-value = 6.8472e-05)