2026-03-15
wooldridge: データセットAER: 操作変数推定 ivreg()car: 線形仮説検定 linearHypothesis()変数 \(z\) が \(x\) の操作変数になるための2条件:
操作変数の外生性 (Instrument Exogeneity) [15.4]: \[\text{Cov}(z, u) = 0\] \(z\) は誤差項と無相関（検定不可能: 経済理論で正当化）
操作変数の関連性 (Instrument Relevance) [15.5]: \[\text{Cov}(z, x) \neq 0\] \(z\) は内生変数 \(x\) と相関（検定可能: 第1段階回帰）
| IV候補 | 外生性 | 関連性 |
|---|---|---|
父親の教育年数 (fatheduc) |
ある程度成立（能力との無相関を仮定） | educと正相関（検定可能） |
兄弟の数 (sibs) |
成立しやすい（能力との無相関） | educと負相関（検定可能） |
誕生四半期 (frstqrt) |
主張される（季節と能力は無関係） | わずかな相関（問題あり） |
| IQ | 外生性を欠く（abilの代理変数） | — |
# educ を fatheduc に単回帰して関連性を確認
first_rel <- lm(educ ~ fatheduc, data = mroz_work)
# fatheduc の係数・SE・t値・p値を表示
coef(summary(first_rel))["fatheduc", ]## Estimate Std. Error t value Pr(>|t|)
## 2.694416e-01 2.858634e-02 9.425538e+00 2.764936e-19# OLS推定
res_ols151 <- lm(lwage ~ educ, data = mroz_work)
# IV推定: 書式 ivreg(被説明変数 ~ 内生変数 | 操作変数)
# fatheduc を educ の IV に指定
res_iv151 <- ivreg(lwage ~ educ | fatheduc, data = mroz_work)
# educ の係数と標準誤差を OLS・IV で比較
rbind(
OLS = coef(summary(res_ols151))["educ", 1:2],
IV = coef(summary(res_iv151))["educ", 1:2]
)## Estimate Std. Error
## OLS 0.10864866 0.01439985
## IV 0.05917348 0.03514177# wage2 データを読み込む
data(wage2)
# 関連性確認: sibs（兄弟数）が多いほど学歴は低いか
# sibs の係数（負なら関連性あり）
lm(educ ~ sibs, data = wage2)$coef["sibs"]## sibs
## -0.2279164# OLS推定
res_ols152 <- lm(lwage ~ educ, data = wage2)
# IV推定（sibs を educ の IV に指定）
res_iv152 <- ivreg(lwage ~ educ | sibs, data = wage2)
# educ の係数と標準誤差を OLS・IV で比較
rbind(
OLS = coef(summary(res_ols152))["educ", 1:2],
IV = coef(summary(res_iv152))["educ", 1:2]
)## Estimate Std. Error
## OLS 0.0598392 0.005963094
## IV 0.1224326 0.026350568\[\text{plim}\, \hat{\beta}_1^{IV} = \beta_1 + \frac{\text{Corr}(z, u)}{\text{Corr}(z, x)} \cdot \frac{\sigma_u}{\sigma_x}\]
# bwght データを読み込む
data(bwght)
# 関連性確認: cigprice（タバコ価格）が packs（喫煙量）と相関しているか
# cigprice の係数を取得（ほぼゼロなら関連性なし）
lm(packs ~ cigprice, data = bwght)$coef["cigprice"]## cigprice
## 0.0002828844# card データセットを読み込む
data(card)
# OLS推定: lwage を educ + コントロール変数に回帰
res_ols154 <- lm(lwage ~ educ + exper + expersq + black + smsa + south,
data = card)
# 2SLS推定: nearc4 が educ の操作変数
# 書式: ivreg(y ~ 内生+外生 | IV+外生)、外生変数は両側に書く
res_iv154 <- ivreg(lwage ~ educ + exper + expersq + black + smsa + south |
nearc4 + exper + expersq + black + smsa + south,
data = card)# educ の係数と標準誤差を OLS・2SLS で比較
rbind(
OLS = coef(summary(res_ols154))["educ", 1:2],
IV = coef(summary(res_iv154))["educ", 1:2]
)## Estimate Std. Error
## OLS 0.07400899 0.003505435
## IV 0.13228884 0.049233236# educ を nearc4 + コントロール変数に回帰
first_stage154 <- lm(educ ~ nearc4 + exper + expersq +
black + smsa + south,
data = card)
# nearc4 の係数・SE・t値・p値を表示
coef(summary(first_stage154))["nearc4", ]## Estimate Std. Error t value Pr(>|t|)
## 3.373208e-01 8.250044e-02 4.088715e+00 4.451508e-05# educ を両親の教育年数 + 経験変数に回帰
first_stage155 <- lm(educ ~ motheduc + fatheduc + exper + expersq,
data = mroz_work)
# F 統計量 > 10 が弱操作変数でない目安（Stock-Yogo, 2005）
linearHypothesis(first_stage155, c("motheduc = 0", "fatheduc = 0"))##
## Linear hypothesis test:
## motheduc = 0
## fatheduc = 0
##
## Model 1: restricted model
## Model 2: educ ~ motheduc + fatheduc + exper + expersq
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 425 2219.2
## 2 423 1758.6 2 460.64 55.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| 左:
内生変数（educ）＋外生変数 | 右:
除外IV（motheduc, fatheduc）＋外生変数# 2SLS推定: | 左: 内生変数+外生変数、| 右: 除外IV+外生変数
res_2sls155 <- ivreg(lwage ~ educ + exper + expersq |
motheduc + fatheduc + exper + expersq,
data = mroz_work)
# 定数・educ・exper・expersq の4係数を表示
summary(res_2sls155)$coefficients[1:4, ]## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0481003069 0.4003280776 0.1201522 0.904419479
## educ 0.0613966287 0.0314366956 1.9530242 0.051474174
## exper 0.0441703929 0.0134324755 3.2883286 0.001091838
## expersq -0.0008989696 0.0004016856 -2.2379930 0.025740027# linearHypothesis の結果を ht に保存して F 値を取り出す
ht <- linearHypothesis(first_stage155, c("motheduc = 0", "fatheduc = 0"))
# ht$F[2] が F 値（[1] は NA のため [2] を使う）、>> 10 → 弱い操作変数問題なし
cat("F統計量:", ht$F[2], "\n")## F統計量: 55.4003| 左に \(IQ\)
を含め、| 右では \(KWW\)
に置き換える## Estimate Std. Error t value Pr(>|t|)
## 1.034136e-01 1.543919e-02 6.698122e+00 3.656739e-11誘導型で \(y_2\) を全外生変数に回帰し残差 \(\hat{v}_2\) を得る
構造方程式に \(\hat{v}_2\) を追加してOLS推定
\(\hat{v}_2\) の係数が有意 \(\Rightarrow\) \(y_2\) は内生的
# Step 1: 誘導型残差を取得（educ の説明されない部分を分離）
v_hat <- residuals(lm(educ ~ motheduc + fatheduc + exper + expersq,
data = mroz_work))
# Step 2: 構造方程式に v_hat を追加（有意なら educ は内生）
endog_reg <- lm(lwage ~ educ + exper + expersq + v_hat, data = mroz_work)
# v_hat の係数・SE・t値・p値を表示
coef(summary(endog_reg))["v_hat", ]## Estimate Std. Error t value Pr(>|t|)
## 0.05816661 0.03480728 1.67110501 0.095440552SLS推定を行い残差 \(\hat{u}_1\) を得る
\(\hat{u}_1\) を全外生変数に回帰し \(R^2\) を取得
検定統計量 \(nR^2 \sim \chi^2_q\)（\(q\) = 過剰識別制約数）
# Step 3: 検定統計量 nR² ～ χ²(q), q = IV数 - 内生変数数 = 2 - 1 = 1
# n × R² を計算
nR2 <- n * summary(oir_reg)$r.squared
# χ²(1) の上側確率（p値）
pval <- pchisq(nR2, df = 1, lower.tail = FALSE)# 2SLS推定（例15.5と同じモデル）
res_hc <- ivreg(lwage ~ educ + exper + expersq |
motheduc + fatheduc + exper + expersq,
data = mroz_work)
# vcovHC(): HC1型の不均一分散頑健分散共分散行列を計算
# coeftest(): 頑健 SE を使って t 検定を再実施、educ 行のみ表示
coeftest(res_hc, vcov = vcovHC(res_hc, type = "HC1"))["educ", ]## Estimate Std. Error t value Pr(>|t|)
## 0.06139663 0.03333859 1.84160854 0.06623070# jtrain データを読み込む
data(jtrain)
# 1988 年のデータを抽出（1階差分変数が格納されている年）
j88 <- subset(jtrain, year == 1988)
# 必要な3変数のみ選択し、欠損値を除去（n = 45 になる）
j88_clean <- na.omit(j88[, c("clscrap", "chrsemp", "cgrant")])
# 2SLS: cgrant（補助金の変化）を chrsemp（訓練時間の変化）のIVとして指定
res_2sls1510 <- ivreg(clscrap ~ chrsemp | cgrant, data = j88_clean)
# chrsemp の係数・SE・t値・p値を表示
coef(summary(res_2sls1510))["chrsemp", ]## Estimate Std. Error t value Pr(>|t|)
## -0.014153164 0.007914742 -1.788202786 0.080790975| 手法 | 状況 | R関数 |
|---|---|---|
| IV (1変数) | 内生変数1個、IV1個 | ivreg(y ~ x | z) |
| 2SLS | 内生変数1個、IV複数 | ivreg(y ~ x+w | z+w) |
| 2SLS（複数内生） | 内生変数複数 | ivreg() |
| 頑健SE付き2SLS | 不均一分散 | coeftest(, vcovHC()) |