固定効果モデルの基礎的な理論と，Rでの実行方法の基本について学んだのでメモ。
理論
パネルデータ
- パネルデータ(panel data)：複数の主体を複数の期間にわたって観察したデータ
- 同じ人を追跡調査したデータとか
- 国×年度の公的統計データとか
- \(k\)個の説明変数\(X\)と被説明変数\(Y\)の主体\(i\)時間\(t\)のパネルデータ：
\[ (X_{1it}, X_{2it},…,X_{kit}, Y_{it}), \hspace{1em} i = 1,…,n \text{ and } t = 1,…,T \]
- balanced panel：各主体と各時間の観測値がすべて揃っているパネルデータ
- unbalanced panel：欠損のあるパネルデータ
固定効果モデル
固定効果モデルはパネルデータの分析で使われる分析手法で，固有効果（個体差とか）のうち推定の邪魔になるもの（固定効果）を除去して分析を行うものです。
固有効果，固定効果，変量効果
- 固有効果\(F_i\)：時間によって変わらない，主体に固有の要素
- 例えば\(Y_{it} = \beta_0 + \beta_1 X_{it} + u_{it}\)というモデルを作るとき，固有効果がある場合は誤差項\(u_{it}\)に含まれ，\(u_{it} = F_i + \varepsilon_{it}\)と分解できる。
- 「固有効果」という語は一般的でないかも。山本（2015）しか使っていないかも。
- 固定効果（fixed effect）：説明変数\(X_{it}\)と独立でない（相関がある）固有効果\(F_i\)
- 問題：固定効果があるとき，誤差項と説明変数が独立でなくなり（\(E(u_{it}|X_{it})\neq 0\)），そのままOLSすると欠落変数バイアスにより誤った推定をしてしまう
- \(F_i\)の効果も\(X_{it}\)の効果も混ぜ合わせて推定してしまう（奥井2015）
- 一致性（consistency, サンプルサイズを大きくすれば，推定値がある値に収束する性質）が得られなくなる（山本2015）
- 対処：固定効果を除去するようなモデルにする（固定効果モデルにする） → 今回の本題
- 問題：固定効果があるとき，誤差項と説明変数が独立でなくなり（\(E(u_{it}|X_{it})\neq 0\)），そのままOLSすると欠落変数バイアスにより誤った推定をしてしまう
- 変量効果（random effect）：説明変数\(X_{it}\)と独立の（相関がない）固有効果\(F_i\)
one-way固定効果モデル
主体の固定効果モデル（within model）の推定方法
\[ Y_{it} = \beta_1 X_{it} +\alpha_i + u_{it} \]
主体の固定効果（entity fixed effects）\(\alpha_i\)の除去については，以下のような推定方法があります。
1. 差分モデル（first difference model）のOLS推定(\(T=2\)のときのみ)
\[ (Y_{i,t+1} - Y_{it}) = (\beta_0 - \beta_0) + \beta_1 (X_{i,t+1}-X_{it}) + (F_i - F_i) + (\varepsilon_{i,t+1} - \varepsilon_{it}) \]
記号を置き換えると \[ \Delta Y_{it} = \beta_1 \Delta X_{it}+ \Delta \varepsilon_{it} \]
- 推定方法：
- 説明変数，被説明変数それぞれ\(t+1\)期から\(t\)期を引く
- 上の式をOLS推定する
2. “\(n-1\)個のダミー説明変数”を用いたOLS推定
- 最小二乗ダミー変数推定（Least Squares Dummy Variables (LSDV) 推定）とも呼ばれるもの
\[ Y_{it} = \beta_0 + \beta_1 X_{it} + \gamma_2 D2_i + \cdots + \gamma_n Dn_i + u_{it} \\ \text{where } D2_i = \begin{cases} 1 & \text{for } i = 2\\ 0 & \text{otherwise} \end{cases} \text{, etc.} \]
3. ”主体の平均除去（Entity-demeaned）”を用いたOLS推定
\[ \begin{align} \tilde{Y}_{it} &= \beta_1 \tilde{X}_{it} + \tilde{u}_{it} \\ \text{where } \tilde{Y}_{it} &= Y_{it} - \bar{Y}_i, \hspace{1em} \bar{Y}_i = \frac{1}{T} \sum^T_{t=1} Y_{it}\\ \tilde{X}_{it} &= X_{it} - \bar{X}_i, \hspace{1em} \bar{X}_i = \frac{1}{T} \sum^T_{t=1} X_{it}\\ \tilde{u}_{it} &= u_{it}- \bar{u}_i, \hspace{1em} \bar{u}_i = \frac{1}{T} \sum_{t=1}^T u_{it} \end{align} \]
時間の固定効果モデル（between）
\[ Y_{it} = \beta_1 X_{it} +\alpha_t + u_{it} \]
時間の固定効果（time fixed effects）\(\alpha_t\)を除きたい場合も主体の場合とほぼ同様です。
- ”\(T-1\)個のダミー説明変数”を用いたOLS推定
- ”time-demeaned”を用いたOLS推定
two-way固定効果モデル
two-way固定効果モデルあるいは主体と時間の固定効果モデル（entity and time fixed effects model） \[ Y_{it} =\beta_1 X_{it} +\alpha_i + \alpha_t + u_{it} \] によって主体の固定効果\(\alpha_i\)と時間の固定効果\(\alpha_t\)の両方を除去したい場合，それぞれの推定方法の組み合わせになります。
- 差分モデルのOLS推定（T=2）
- ”entity demeaning”と”\(T-1\)個のダミー変数”を用いたOLS推定
- ”time demeaning”と“\(n-1\)個のダミー変数”を用いたOLS推定
- ”entity & time demeaning”を用いたOLS推定
- 説明変数と被説明変数について，主体と時間両方の平均を引いてOLS推定
固定効果モデルの仮定
説明変数\(X\)が１つ，主体の固定効果のみのモデルの場合について表記すると \[ Y_{it} = β_1 X_{it} + α_i + u_{it}, i = 1,…,n, t = 1,…, T \]
- 仮定１： \(E(u_{it}|X_{i1},…,X_{iT},α_i) = 0\)
- 仮定２： \((X_{i1},…,X_{iT},u_{i1},…,u_{iT}), i =1,…,n\) は独立かつ同一の同時分布に従う
- 仮定３： 大きな外れ値は存在しない：\((X_{it},u_{it})\)はゼロでない有限の4次のモーメントを持つ
- 仮定４： 完全な多重共線性は存在しない（Xが複数ある場合）
- 仮定５：各主体についての誤差項は，説明変数の下で，自己相関なし：\(\text{cov}(u_{it},u_{is}|X_{i1},…,X_{iT}, \alpha_i) =0, \ t \neq s\)
- この仮定が満たされない場合は頑健標準誤差を使う
Rによる実践
library(plm)
data("Gasoline", package="plm")
head(Gasoline)
> head(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}のGasolineデータセットを使います。
変数の意味は以下の通り：
- country：18ヶ国の国名
- year：年（1960～1978）
- lgaspcar：自動車1台あたりのガソリン消費量の自然対数
- lincomep：一人あたり実質所得の自然対数
- lrpmg：実質ガソリン価格の自然対数
- lcarpcap：一人あたりの自動車の量（ストック）の自然対数
ガソリン価格の変化が１台あたりのガソリン消費量に与えた影響（エコカー化に与えた影響）を考える，という状況を想定してみます。
lm()による推定
ダミー変数推定（LSDV推定）
factor型にして{makedummies}でダミー変数を作るか，変数をfactor型にしてそのままlmに入れます
library(tidyverse) library(makedummies) # ダミー変数を作成 Gasoline_dummies <- Gasoline %>% mutate(year = as.factor(year)) %>% select(country, year) %>% makedummies() # 作成したダミー変数を結合 data <- Gasoline %>% select(-country, -year) %>% bind_cols(Gasoline_dummies) # 推定 lm_LSDV <- lm(lgaspcar ~ ., data = data) summary(lm_LSDV)
（コード文中の有意水準「***」がはてなmarkdownの太字指定とみなされてレイアウトが崩れるため，summaryの返り値を画像で貼り付けています）
と出ました。
{lmtest}でロバスト標準誤差を計算してみます。
library(lmtest) coeftest(lm_LSDV, vcov = vcovHC(lm_LSDV, type = "HC1"))
平均除去（demean）による推定
# 主体固定効果の除去 entity_demeaned <- Gasoline %>% group_by(country) %>% mutate(lgaspcar = lgaspcar - mean(lgaspcar), lincomep = lincomep - mean(lincomep), lrpmg = lrpmg - mean(lrpmg), lcarpcap = lcarpcap - mean(lcarpcap)) %>% ungroup() # 時間固定効果の除去 twoways_demeaned <- entity_demeaned %>% group_by(year) %>% mutate(lgaspcar = lgaspcar - mean(lgaspcar), lincomep = lincomep - mean(lincomep), lrpmg = lrpmg - mean(lrpmg), lcarpcap = lcarpcap - mean(lcarpcap)) %>% ungroup() # 推定 lm_demeaned <- lm(lgaspcar ~ -1 + lrpmg + lincomep + lcarpcap, data = twoways_demeaned) # クラスターロバスト標準誤差 library(plm) coeftest(lm_demeaned, vcov = plm::vcovHC(lm_demeaned, method = "arellano", type = "HC1", cluster = c("group", "time")))
となりました。理論通り，ダミー変数を使っての推定と同じ結果です。
plm()による推定
まずpdata.frame()でデータをplm用の形式にします。
重要な引数は以下の２つ：
index：個体，時間の列名を文字列ベクトルで指定するdrop.index：pdata.frame変換時にindex列を除くか否か
library(plm) Gasoline_pdf <- pdata.frame(Gasoline, index=c("country","year"), drop.index=TRUE) head(Gasoline_pdf)
> Gasoline_pdf <- pdata.frame(Gasoline, 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()で推定します。
重要な引数は次の２つです：
method：モデルを指定しますfd：差分モデル（first-difference model）within：主体固定効果モデルbetween：時間固定効果モデル
effect：推定する固定効果を指定しますindividual：主体固定効果time：時間固定効果twoways：主体・時間の両方
plm_result <- plm(lgaspcar ~ lrpmg + lincomep + lcarpcap, data = Gasoline_pdf, method = "within", effect = "twoways") summary(plm_result)
coeftest(plm_result, vcov = plm::vcovHC(plm_result, method = "arellano", type = "HC1", cluster = c("group", "time")))
と出ました。lmの推定量と一致しました。
しかし，平均除去によるOLS推定を行ったときの結果と標準誤差の値が違いますね。なんでだろう…ここの解明は今後の課題です。
まとめ
- 固定効果モデルは主体（個体）や時間（年度とか）に固有の効果が推定の邪魔になるときに，その固有の効果（固定効果）を除去する手法です
- ダミー変数を作って固有の効果を吸収させたりとかします
- Rで実行するときは{plm}パッケージを使うと楽です。
参考文献
- Angrist, J. D., & Pischke, J. S. (2013)『「ほとんど無害」な計量経済学』
- James H. Stock & Mark W. Watson (2016) 『入門計量経済学』
- Panel Data Econometrics in R: The plm Package - R/CRAN - R Project
- 太田浩司. (2013). パネル・データ分析におけるクラスター頑健手法の使用について.
- 奥井亮. (2015). 固定効果と変量効果. 日本労働研究雑誌, (657), 6-9.
- 山本勲（2015）『実証分析のための計量経済学』，中央経済社。
