必要なパッケージの読み込み

wooldridgeパッケージの読み込み

library(wooldridge)

plmパッケージの読み込み（パネルデータ分析に使用）

library(plm)

14-1 固定効果推定

固定効果変換 (Fixed Effects Transformation)

モデル: \(y_{it} = \beta_1 x_{it} + a_i + u_{it}, \quad t = 1, 2, \dots, T\) [14.1]
個体 \(i\) に関して時間平均をとると： \[\bar{y}_i = \beta_1 \bar{x}_i + a_i + \bar{u}_i \hfill [14.2]\]
[14.1] から [14.2] を引く（固定効果変換 / within変換）： \[\ddot{y}_{it} = \beta_1 \ddot{x}_{it} + \ddot{u}_{it}, \quad t = 1, 2, \dots, T \hfill [14.3]\]
\(\ddot{y}_{it} = y_{it} - \bar{y}_i\)：時間平均除去データ (time-demeaned data)
\(a_i\) が消去される \(\Rightarrow\) OLSで \(\beta_1\) を一致推定できる

固定効果推定量の特性

[14.3] を Pooled OLS で推定したものが 固定効果推定量（within推定量）
時間不変の変数（出生地、性別など）は時間平均除去によって消去され、推定不可能
厳密外生性仮定（\(\text{E}(u_{it}|\mathbf{X}_i, a_i) = 0\)）のもとで不偏・一致
\(a_i\) と \(x_{it}\) が任意に相関していても構わない（第13章のFDと同じ強み）
自由度: \(df = NT - N - k = N(T-1) - k\)（時間平均除去により個体ごとに自由度を1失う）

複数の説明変数への拡張

一般的な観測不能効果モデル [14.4]: \[y_{it} = \beta_1 x_{it1} + \beta_2 x_{it2} + \cdots + \beta_k x_{itk} + a_i + u_{it}\]
各変数を時間平均除去 \(\Rightarrow\) [14.5]: \[\ddot{y}_{it} = \beta_1 \ddot{x}_{it1} + \beta_2 \ddot{x}_{it2} + \cdots + \beta_k \ddot{x}_{itk} + \ddot{u}_{it}\]
これをPooled OLSで推定する

例14.1: 職業訓練が不良品率に与える影響 (JTRAIN)

データ: 54社 × 3年（1987〜1989年）
従属変数: \(\log(scrap_{it})\)（不良品率の対数）
説明変数: 年ダミー \(d88\), \(d89\), 補助金 \(grant\), ラグ補助金 \(grant_{-1}\)
\(df = 54(3-1) - 4 = 104\)

data(jtrain)
pdata_jt <- pdata.frame(jtrain, index = c("fcode", "year"))
res_fe141 <- plm(lscrap ~ d88 + d89 + grant + grant_1,
                 data = pdata_jt, model = "within",
                 na.action = na.omit)
summary(res_fe141)$coefficients

##            Estimate Std. Error    t-value   Pr(>|t|)
## d88     -0.08021567  0.1094751 -0.7327297 0.46537158
## d89     -0.24720279  0.1332183 -1.8556220 0.06633916
## grant   -0.25231487  0.1506290 -1.6750751 0.09692392
## grant_1 -0.42158951  0.2102000 -2.0056593 0.04748974

例14.1: 結果の解釈

\(\hat{\beta}_{grant_{-1}} \approx -0.422\)：前年に補助金を受けると不良品率が約34%低下 [\(\exp(-0.422) - 1 \approx -0.344\)]、5%水準で有意
\(\hat{\beta}_{grant} \approx -0.252\)：当年補助金は10%水準で有意（即時効果は小さい）
\(d89\) の係数が負 \(\Rightarrow\) 補助金とは無関係に1989年の生産性は上昇 \(\Rightarrow\) 年ダミーを含めることの重要性

例14.2: 教育収益率の時間変化 (WAGEPAN)

データ: 545人 × 8年（1980〜1987年）、WAGEPAN
固定効果では educ（時間不変）は推定不可 \(\Rightarrow\) educ × 年ダミー の交差項を使用

data(wagepan)
pdata_wp <- pdata.frame(wagepan, index = c("nr", "year"))
res_fe142 <- plm(lwage ~ married + union + exper + expersq +
                   d81 + d82 + d83 + d84 + d85 + d86 + d87 +
                   I(educ * d81) + I(educ * d82) + I(educ * d83) +
                   I(educ * d84) + I(educ * d85) + I(educ * d86) +
                   I(educ * d87),
                 data = pdata_wp, model = "within")
# 交差項の係数のみ表示
coef(res_fe142)[grep("educ", names(coef(res_fe142)))]

## I(educ * d81) I(educ * d82) I(educ * d83) I(educ * d84) I(educ * d85) 
##   0.004990619   0.001650951  -0.002662147  -0.009825689  -0.009214474 
## I(educ * d86) I(educ * d87) 
##  -0.012138164  -0.015789160

例14.2: 解釈

\(d87 \cdot educ\) の係数 \(\approx 0.030\)（\(t \approx 2.48\)）：最大
1987年の教育収益率は基準年（1980年）より約3%ポイント高い
交差項7本の結合有意性検定: \(p値 \approx 0.28\)（まとめて有意でない）
教育収益率は1980〜1987年を通じてほぼ一定という結論

14-1a ダミー変数回帰

ダミー変数回帰 (Dummy Variable Regression)

\(a_i\) を各個体のパラメータとして推定する伝統的な方法
\(N\) 個の切片ダミー変数を説明変数に含めてOLS推定 \(\Rightarrow\) ダミー変数回帰
重要な性質: ダミー変数回帰は固定効果推定量と全く同じ \(\hat{\beta}_j\) を与える
固定効果推定後、各個体の固定効果の推定値 [14.6]: \[\hat{a}_i = \bar{y}_i - \hat{\beta}_1 \bar{x}_{i1} - \cdots - \hat{\beta}_k \bar{x}_{ik}\]
\(N\) が大きい場合、ダミー変数回帰は非実用的（変数が多すぎる）

14-1b 固定効果 vs 1階差分

固定効果法 (FE) と 1階差分法 (FD) の比較

特徴	FE (固定効果)	FD (1階差分)
\(T=2\) のとき	同一	同一
\(T \geq 3\) のとき	異なる推定量	異なる推定量
\(u_{it}\) が系列無相関	効率的	非効率
\(u_{it}\) がランダムウォーク	非効率	効率的
測定誤差に対する感度	高い	高い

\(T = 2\): FEとFDは同一 \(\Rightarrow\) どちらでも良い
\(T \geq 3\): \(u_{it}\) の系列相関構造で選択
- 系列無相関 \(\Rightarrow\) FEが優れる
- ランダムウォーク的正の系列相関 \(\Rightarrow\) FDが優れる

14-1c 不均衡パネル

不均衡パネルの固定効果推定

不均衡パネル (unbalanced panel): 一部の個体で欠損期間がある
個体 \(i\) の観測期間数を \(T_i\) とすると、総観測数 \(= T_1 + T_2 + \cdots + T_N\)
時間平均除去は個体ごとに \(T_i\) 期間の平均を使う
固定効果推定は不均衡パネルでも同様に適用可能
注意: 欠損が自己選択的（離脱が誤差項と相関）な場合はバイアスが生じる
- ただしFEは \(a_i\) と選択指標の相関を許容するため、FDより頑健な場合がある

例14.3: 不均衡パネルでの職業訓練効果 (JTRAIN)

sales（売上の対数）と employ（従業員数の対数）を追加
3社が完全に脱落、5観測が追加欠損 \(\Rightarrow\) \(n = 148\)

res_fe143 <- plm(lscrap ~ d88 + d89 + grant + grant_1 +
                   lsales + lemploy,
                 data = pdata_jt, model = "within",
                 na.action = na.omit)
summary(res_fe143)$coefficients[c("grant", "grant_1"), ]

##           Estimate Std. Error   t-value   Pr(>|t|)
## grant   -0.2967542  0.1570861 -1.889119 0.06206040
## grant_1 -0.5355783  0.2242060 -2.388778 0.01896951

\(\hat{\beta}_{grant} \approx -0.297\)、\(\hat{\beta}_{grant_{-1}} \approx -0.536\)：基本結果と整合的

14-2 変量効果モデル

変量効果モデルの設定

同じ不観測効果モデルから出発 [14.7]: \[y_{it} = \beta_0 + \beta_1 x_{it1} + \cdots + \beta_k x_{itk} + a_i + u_{it}\]
変量効果 (RE) の仮定 [14.8]: \[\text{Cov}(x_{itj}, a_i) = 0, \quad t = 1, \dots, T;\ j = 1, \dots, k\]
\(a_i\) がすべての説明変数と無相関と仮定 \(\Rightarrow\) \(a_i\) を誤差項に含める
合成誤差項: \(v_{it} = a_i + u_{it}\)
\(a_i\) が各期間の誤差に入るため、\(v_{it}\) は時間を通じて正の系列相関を持つ

GLS変換（準時間平均除去）

合成誤差の系列相関を除去するため一般化最小二乗法 (GLS) を適用
変換パラメータ [14.10]: \[\theta = 1 - \left[\sigma_u^2 / (\sigma_u^2 + T\sigma_a^2)\right]^{1/2}\]
変換後の方程式 [14.11]: \[y_{it} - \theta\bar{y}_i = \beta_0(1-\theta) + \beta_1(x_{it1} - \theta\bar{x}_{i1}) + \cdots + (v_{it} - \theta\bar{v}_i)\]
quasi-demeaned（準時間平均除去）データに基づくPooled OLS \(\Rightarrow\) 変量効果推定量
\(\theta = 0\): Pooled OLS と同一
\(\theta = 1\): 固定効果推定量と同一

例14.4: 賃金方程式 (WAGEPAN)

3手法（Pooled OLS, RE, FE）の比較

# Pooled OLS（参考）
res_pols <- plm(lwage ~ educ + black + hisp + exper + expersq +
                  married + union +
                  d81 + d82 + d83 + d84 + d85 + d86 + d87,
                data = pdata_wp, model = "pooling")
# 変量効果 (RE)
res_re <- plm(lwage ~ educ + black + hisp + exper + expersq +
                married + union +
                d81 + d82 + d83 + d84 + d85 + d86 + d87,
              data = pdata_wp, model = "random")
# 固定効果 (FE)
res_fe <- plm(lwage ~ expersq + married + union +
                d81 + d82 + d83 + d84 + d85 + d86 + d87,
              data = pdata_wp, model = "within")

例14.4: 結果の比較

# 時間変動変数の係数比較
vars <- c("expersq", "married", "union")
rbind(
  POLS = coef(res_pols)[vars],
  RE   = coef(res_re)[vars],
  FE   = coef(res_fe)[vars]
)

##           expersq    married      union
## POLS -0.002411703 0.10825295 0.18246128
## RE   -0.004723943 0.06398602 0.10613443
## FE   -0.005185498 0.04668036 0.08000186

例14.4: 解釈

\(\hat{\theta} \approx 0.643\)：REはPooled OLSとFEの中間
married の係数: POLS(0.108) \(>\) RE(0.064) \(>\) FE(0.047)
- 能力の高い人が結婚しやすい \(\Rightarrow\) POLS/REは上方バイアス
union の係数: POLS(0.182) \(>\) RE(0.106) \(>\) FE(0.080)
- 組合効果の一部は個人固有特性に起因
FEでは educ, black, hispan は時間不変のため推定不可

14-2a 変量効果 vs Pooled OLS

変量効果とPooled OLSの選択

Breusch-Pagan (1980) 検定: \(H_0: \sigma_a^2 = 0\)（個体効果が存在しない）
\(\sigma_a^2 = 0\) なら観測不能効果なし \(\Rightarrow\) Pooled OLSが適切
ただし、この検定はあくまで \(\sigma_a^2 > 0\) を確認するもの
- \(\sigma_a^2 > 0\) でも \(a_i\) が \(x_{it}\) と相関していれば RE/POLS ともに不一致
実務上の含意: REはPooled OLSより好まれることが多い（特に誤差構造がRE仮定に近い場合）
- 系列相関の除去により効率的
- \(a_i\) の一部を誤差項から除去 \(\Rightarrow\) バイアスを軽減

14-2b 変量効果 vs 固定効果

FE vs RE の選択基準

FE: \(a_i\) と \(x_{it}\) が相関していてもよい \(\Rightarrow\) より頑健
RE: \(a_i\) と \(x_{it}\) が無相関という仮定が必要 \(\Rightarrow\) 効率的（標準誤差が小さい）
RE を選ぶ場合: 主要な説明変数が時間不変（例: educ）\(\Rightarrow\) FEでは推定不可
Hausman (1978) 検定: RE仮定 [14.8] が成立しているか検定

\[H_0: \text{Cov}(x_{itj}, a_i) = 0 \quad \forall j\]

Hausman検定の実施

# Hausman検定（同じ変数セットでFEとREを推定）
res_re_h <- plm(lwage ~ expersq + married + union +
                  d81 + d82 + d83 + d84 + d85 + d86 + d87,
                data = pdata_wp, model = "random")
res_fe_h <- plm(lwage ~ expersq + married + union +
                  d81 + d82 + d83 + d84 + d85 + d86 + d87,
                data = pdata_wp, model = "within")
phtest(res_fe_h, res_re_h)

## 
##  Hausman Test
## 
## data:  lwage ~ expersq + married + union + d81 + d82 + d83 + d84 + d85 +  ...
## chisq = 37.01, df = 10, p-value = 5.637e-05
## alternative hypothesis: one model is inconsistent

\(p値\) が小さい \(\Rightarrow\) RE仮定を棄却 \(\Rightarrow\) FEを採用

14-3 相関変量効果アプローチ

CREアプローチの基本アイデア

\(a_i\) が観測変数と相関することを明示的にモデル化
\(a_i\) と \(\{x_{it}\}\) の相関を、時間平均 \(\bar{x}_i\) との線形関係で表す [14.12]: \[a_i = \alpha + \gamma \bar{x}_i + r_i, \quad \text{Cov}(\bar{x}_i, r_i) = 0\]
これを元のモデルに代入 [14.14]: \[y_{it} = \alpha + \beta x_{it} + \gamma \bar{x}_i + r_i + u_{it}\]
\(\bar{x}_i\) を追加してRE推定 \(\Rightarrow\) 相関変量効果 (CRE) 推定量

CREとFEの等価性

重要な結果 [14.16]: CRE推定量はFE推定量と同一 \[\hat{\beta}_{CRE} = \hat{\beta}_{FE}\]
CREアプローチの利点:
1. 時間不変変数の係数を推定できる（\(\delta_j\) として）
2. RE vs FE 検定の形式的な枠組みを提供 (\(H_0: \gamma = 0\))
3. 不均衡パネルへの拡張が容易

CREの推定: WAGEPAN

# 時間変動変数の個体平均を計算
wagepan$mean_expersq <- ave(wagepan$expersq, wagepan$nr)
wagepan$mean_married <- ave(wagepan$married, wagepan$nr)
wagepan$mean_union   <- ave(wagepan$union, wagepan$nr)

pdata_wp2 <- pdata.frame(wagepan, index = c("nr", "year"))

res_cre <- plm(lwage ~ expersq + married + union +
                 mean_expersq + mean_married + mean_union +
                 d81 + d82 + d83 + d84 + d85 + d86 + d87,
               data = pdata_wp2, model = "random")
# 時間平均変数の係数（γの推定値）
coef(res_cre)[c("mean_expersq", "mean_married", "mean_union")]

## mean_expersq mean_married   mean_union 
##  0.003212666  0.161764098  0.161241547

CRE: RE vs FE の検定

\(H_0: \gamma_1 = \gamma_2 = \cdots = \gamma_k = 0\)（REで十分）[14.19]
時間平均変数 \(\bar{x}_{i1}, \dots, \bar{x}_{ik}\) が全て統計的に有意でなければ RE を採用

# 時間平均変数の結合有意性検定（Wald検定）
library(car)
linearHypothesis(res_cre,
                 c("mean_expersq = 0",
                   "mean_married = 0",
                   "mean_union = 0"))

## 
## Linear hypothesis test:
## mean_expersq = 0
## mean_married = 0
## mean_union = 0
## 
## Model 1: restricted model
## Model 2: lwage ~ expersq + married + union + mean_expersq + mean_married + 
##     mean_union + d81 + d82 + d83 + d84 + d85 + d86 + d87
## 
##   Res.Df Df  Chisq Pr(>Chisq)    
## 1   4349                         
## 2   4346  3 35.992  7.516e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

14-4 パネルデータを用いた政策分析

一般的な政策分析フレームワーク [14.22]

\[y_{it} = \eta_1 + \alpha_2 d2_t + \cdots + \alpha_T dT_t + \beta w_{it} + \mathbf{x}_{it}\boldsymbol{\psi} + a_i + u_{it}\]

\(w_{it}\): 介入指標（ \(w_{it} = 0\) または \(1\)）
\(\beta\): 政策効果（関心パラメータ）
\(\mathbf{x}_{it}\): その他の時変コントロール変数
FEまたはFDで推定（\(a_i\) と \(w_{it}\) の相関を許容）
DiD（第13章）は \(T = 2\) の特殊ケース

動的効果の推定 [14.24]

\(t-1\)、\(t-2\) 期までのラグ効果を含む動的モデル: \[y_{it} = \eta_1 + \sum_{t=2}^{T}\alpha_t dt_t + \beta_0 w_{it} + \beta_1 w_{i,t-1} + \beta_2 w_{i,t-2} + \mathbf{x}_{it}\boldsymbol{\psi} + a_i + u_{it}\]
長期効果: \(\beta_0 + \beta_1 + \beta_2\)

因果分析の注意点

\(w_{it}\) が過去の誤差項に反応する場合（フィードバック）\(\Rightarrow\) 厳密外生性違反
検定方法: \(w_{i,t+1}\)（次期の政策）を追加し、有意かどうかを確認 [14.25]
- 有意 \(\Rightarrow\) フィードバックの存在、FE/FDの結果に注意が必要
偽造検定 (Falsification Test): 例えば、今年の介入が去年の成果を予測すべきでない

異質トレンドモデル (Heterogeneous Trend Model) [14.26]

個体固有の時間トレンド \(g_i t\) が存在する場合: \[y_{it} = \eta_1 + \sum_{t=2}^T \alpha_t dt_t + \beta w_{it} + \mathbf{x}_{it}\boldsymbol{\psi} + a_i + g_i t + u_{it}\]
1階差分で \(a_i\) を除去 \(\Rightarrow\) \(g_i\) が残る [14.27]: \[\Delta y_{it} = \sum_{t=2}^T \alpha_t \Delta dt_t + \beta \Delta w_{it} + \Delta\mathbf{x}_{it}\boldsymbol{\psi} + g_i + \Delta u_{it}\]
FEを再度適用して \(g_i\) を除去 \(\Rightarrow\) 最低 \(T \geq 3\) が必要

14-5 他のデータ構造への応用

マッチドペアとクラスターサンプル

パネルデータ手法は時間以外の構造にも応用可能
マッチドペア (matched pairs): 例) 兄弟・双子ペアの差分
- 家族固有効果 \(a_f\) を消去
- 例: 十代妊娠が将来所得に与える影響（姉妹ペア差分）
クラスターサンプル (cluster sample): 個体がクラスター（学校、企業）から抽出
- クラスターが「個体」、クラスター内観測が「時点」に対応
- クラスター効果を固定効果または差分で除去

クラスターサンプルにおける注意点

政策変数 \(w_{it}\) がクラスター内で変動しない場合（例: 教師の質の効果） \(\Rightarrow\) 固定効果は使えない \(\Rightarrow\) 変量効果（または Pooled OLS + クラスター標準誤差）
クラスターを無視してOLSを使う場合 \(\Rightarrow\) 標準誤差が過小（過信） \(\Rightarrow\) クラスター頑健標準誤差 (cluster-robust SE) を使う
Pooled OLSでクラスター頑健標準誤差を使うことで対処可能

クラスター頑健標準誤差の適用場面

クラスタリングが正当化されるのは真のクラスターサンプルの場合
- 無作為サンプルをグループに事後分割しても正当化されない
政策変数がグループレベルで変化する場合（例: 最低賃金法）
- 処置の割り当て不確実性を考慮するためクラスタリングが適切

第14章のまとめ

手法	\(a_i\) と \(x_{it}\) の相関	時間不変変数	効率性
Pooled OLS	不可	推定可	低
固定効果 (FE)	許容	推定不可	中
1階差分 (FD)	許容	推定不可	中
変量効果 (RE)	不可（仮定）	推定可	高
相関変量効果 (CRE)	許容	推定可	中

まとめ（続き）

FE vs FD: \(u_{it}\) が系列無相関 \(\Rightarrow\) FE優位；ランダムウォーク \(\Rightarrow\) FD優位
RE vs FE: Hausman検定または CRE アプローチで選択
CRE: FE の点推定を維持しつつ、時間不変変数の係数と RE vs FE の検定を提供
政策分析: 処置指標を適切に定義すれば DiD はFEの特殊ケース；動的効果・フィードバック検定も重要
他のデータ構造: 兄弟ペア、クラスターサンプルにも同様の手法が適用可能

第14章: 発展的なパネルデータ法

Jeffrey Wooldridge (2016). Introductory Econometrics: A Modern Approach Seventh Edition. Cengage Learning.

準備

必要なパッケージの読み込み

14-1 固定効果推定

固定効果変換 (Fixed Effects Transformation)

固定効果推定量の特性

複数の説明変数への拡張

例14.1: 職業訓練が不良品率に与える影響 (JTRAIN)

例14.1: 結果の解釈

例14.2: 教育収益率の時間変化 (WAGEPAN)

例14.2: 解釈

14-1a ダミー変数回帰

ダミー変数回帰 (Dummy Variable Regression)

14-1b 固定効果 vs 1階差分

固定効果法 (FE) と 1階差分法 (FD) の比較

14-1c 不均衡パネル

不均衡パネルの固定効果推定

例14.3: 不均衡パネルでの職業訓練効果 (JTRAIN)

14-2 変量効果モデル

変量効果モデルの設定

GLS変換（準時間平均除去）

例14.4: 賃金方程式 (WAGEPAN)

例14.4: 結果の比較

例14.4: 解釈

14-2a 変量効果 vs Pooled OLS

変量効果とPooled OLSの選択

14-2b 変量効果 vs 固定効果

FE vs RE の選択基準

Hausman検定の実施

14-3 相関変量効果アプローチ

CREアプローチの基本アイデア

CREとFEの等価性

CREの推定: WAGEPAN

CRE: RE vs FE の検定

14-4 パネルデータを用いた政策分析

一般的な政策分析フレームワーク [14.22]

動的効果の推定 [14.24]

因果分析の注意点

異質トレンドモデル (Heterogeneous Trend Model) [14.26]

14-5 他のデータ構造への応用

マッチドペアとクラスターサンプル

クラスターサンプルにおける注意点

クラスター頑健標準誤差の適用場面

まとめ

第14章のまとめ

まとめ（続き）

Jeffrey Wooldridge (2016).
Introductory Econometrics: A Modern Approach
Seventh Edition. Cengage Learning.