library(tidyverse)

1 Data

Heiss (2020) Using R for Introductory Econometrics, Section 15.1 に依拠して，Mroz (ECMA 1987) の既婚女性に関する賃金データ mroz を利用する．

このデータは Wooldridge の教科書で使用されており，wooldridge パッケージに収められている．

主な変数

wage: estimated wage rate (per hour), 賃金率
- inlf = 0 (in labor force) の個人は NA
educ: years of schooling, 教育年数
exper: actual labor market experience
motheduc: mother’s years of schooling, 母親の教育年数
fatheduc: father’s years of schooling, 父親の教育年数

# install.packages("wooldridge")
data(mroz, package = "wooldridge")  # load data 
nrow(mroz)  # sample size

## [1] 753

hist(mroz$wage)

summary(mroz$wage)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.1282  2.2626  3.4819  4.1777  4.9708 25.0000     325

推定には fixest パッケージを用いる． fixest はパネルデータ分析だけではなく操作変数推定も可能．使い方はパッケージの公式Vignette：Fast Fixed-Effects Estimation: Short Introduction を参照．

# install.packages("fixest")
library(fixest)

注：以前は AER::ivreg が主流だった．

2 IV estimation: One endogenous, one instrument

\[ \mbox{First stage:} \quad x = \gamma_0 + \gamma_1 z + v \] \[ \mbox{Second stage:} \quad y = \beta_0 + \beta_1 x + u \]

IV 推定量

\[ \hat{\beta}_{IV} = \frac{cov(z, y)}{cov(z, x)} \]

mroz データを用いて教育の収益率を推定する．

アウトカム log(wage) を educ に回帰したとき，educ の回帰係数 $\beta$ が半弾力性（educ が1単位増加したとき賃金率が $100 \times \beta$ % 増加する）として解釈できる
しかし educ は内生変数なので，fatheduc で instrument する

\[ \mbox{First stage:} \quad (\mbox{Education}) = \gamma_0 + \gamma_1 (\mbox{Father's education}) + v \] \[ \mbox{Second stage:} \quad \log(\mbox{Wage}) = \beta_0 + \beta_1 (\mbox{Education}) + u \]

2.1 Estimation by hand

mroz_with_complete_wage <- mroz %>% 
  filter(!is.na(wage))
y <- log(mroz_with_complete_wage$wage)
x <- mroz_with_complete_wage$educ
z <- mroz_with_complete_wage$fatheduc
cov(z, y) / cov(z, x)

## [1] 0.05917348

IV推定量は教育年数が1年増加することで賃金率が 6% 上昇することを示す．ただし，この計算からは統計的に有意かどうかは分からない．

3 2SLS (Two-stage least squares): One endogenous, one instrument

3.1 Estimation by hand

3.1.1 First stage

内生変数 x を操作変数 z に回帰し fitted value $\hat{x}$ を計算する．

この例では first stage の説明変数は z のみなので，F 値は summary(lmオブジェクト) で確認できる
Fitted value は lm オブジェクトに $fitted.values として格納されている
lm オブジェクトを構成する要素は str(lmオブジェクト) で確認できる

first_stage <- lm(x ~ z)
summary(first_stage)

## 
## Call:
## lm(formula = x ~ z)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4704 -1.1231 -0.1231  0.9546  5.9546 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.23705    0.27594  37.099   <2e-16 ***
## z            0.26944    0.02859   9.426   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.081 on 426 degrees of freedom
## Multiple R-squared:  0.1726, Adjusted R-squared:  0.1706 
## F-statistic: 88.84 on 1 and 426 DF,  p-value: < 2.2e-16

x_hat <- first_stage$fitted.values
summary(x_hat)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.24   12.12   12.12   12.66   13.47   14.82

計算された回帰係数から計算して確認する．

\[ \hat{x} = \hat{\gamma}_0 + \hat{\gamma}_1 z \]

first_stage$coefficients

## (Intercept)           z 
##  10.2370514   0.2694416

x_hat2 <- first_stage$coefficients[1] + first_stage$coefficients[2] * z  # fitted value 
summary(x_hat2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.24   12.12   12.12   12.66   13.47   14.82

3.1.2 Second stage

アウトカム y を first stage で計算された $\hat{x}$ に回帰する．

lm(y ~ x_hat)

## 
## Call:
## lm(formula = y ~ x_hat)
## 
## Coefficients:
## (Intercept)        x_hat  
##     0.44110      0.05917

summary 関数を使えば x_hat 変数の標準誤差が出力されるが，これは不正確であるため，以下のセクションのように専用の関数を使用するべき．

3.2 `fixest::feols` 関数で推定

fixest::feols 関数では アウトカム ~ 外生変数 | 内生変数 ~ 操作変数 の形式で変数を指定する．

今回は外生変数がないため以下のように 1 と書く．

iv1 <- feols(log(wage) ~ 1 | educ ~ fatheduc, data = mroz)
iv1

## TSLS estimation, Dep. Var.: log(wage), Endo.: educ, Instr.: fatheduc
## Second stage: Dep. Var.: log(wage)
## Observations: 428 
## Standard-errors: IID 
##             Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 0.441103   0.446102 0.988795 0.323325    
## fit_educ    0.059173   0.035142 1.683850 0.092943 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.687777   Adj. R2: 0.09131
## F-test (1st stage), educ: stat = 88.8    , p < 2.2e-16 , on 1 and 426 DoF.
##               Wu-Hausman: stat =  2.47035, p = 0.116756, on 1 and 425 DoF.

First stage は以下で表示できる．

F 統計量もレポートされる． F 値が10を超えているので，内生変数と操作変数の関連性条件は満たされているだろう．

summary(iv1, stage = 1)

## TSLS estimation, Dep. Var.: educ, Endo.: educ, Instr.: fatheduc
## First stage: Dep. Var.: educ
## Observations: 428 
## Standard-errors: IID 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 10.237051   0.275936 37.09933 < 2.2e-16 ***
## fatheduc     0.269442   0.028586  9.42554 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 2.07643   Adj. R2: 0.170617
## F-test (1st stage): stat = 88.8, p < 2.2e-16, on 1 and 426 DoF.

feols 関数では vcov = "hetero" 引数を指定することで不均一分散に頑健な標準誤差を計算できる．

以下のように操作変数を使用しない場合を含めて比較する．頑健な標準誤差を用いると，教育年数変数は有意ではない．

ols1 <- feols(log(wage) ~ educ, data = mroz)
iv2 <- feols(log(wage) ~ 1 | educ ~ fatheduc, data = mroz, vcov = "hetero")
library(modelsummary)
modelsummary(list(OLS = ols1, summary(iv1, stage = 1:2), summary(iv2, stage = 1:2)), 
             gof_omit = "AIC|BIC|RMSE", stars = TRUE)

	OLS	iv: First stage: educ	iv: Second stage	iv: First stage: educ	iv: Second stage
(Intercept)	-0.185	10.237***	0.441	10.237***	0.441
	(0.185)	(0.276)	(0.446)	(0.272)	(0.465)
educ	0.109***
	(0.014)
fatheduc		0.269***		0.269***
		(0.029)		(0.029)
fit_educ			0.059+		0.059
			(0.035)		(0.037)
Num.Obs.	428	428	428	428	428
R2	0.118	0.173	0.093	0.173	0.093
R2 Adj.	0.116	0.171	0.091	0.171	0.091
Std.Errors	IID	IID	IID	Heteroskedasticity-robust	Heteroskedasticity-robust
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

注：「OLS」とはモデルをデータに当てはめる原理の名称であり，「IV」の counterpart を表すラベルとして使うことは本来は適切ではない．しかし，実際には（少なくとも経済学の論文では）操作変数を使っていないナイーブなモデルという意味で「OLS」というラベルが便宜的に用いられている．

外生変数として exper などを追加．

iv3 <- feols(log(wage) ~ exper + I(exper^2) | educ ~ fatheduc, data = mroz, vcov = "hetero")
iv3

## TSLS estimation, Dep. Var.: log(wage), Endo.: educ, Instr.: fatheduc
## Second stage: Dep. Var.: log(wage)
## Observations: 428 
## Standard-errors: Heteroskedasticity-robust 
##              Estimate Std. Error   t value  Pr(>|t|)    
## (Intercept) -0.061117   0.458134 -0.133404 0.8939372    
## fit_educ     0.070226   0.035939  1.954043 0.0513532 .  
## exper        0.043672   0.015566  2.805513 0.0052549 ** 
## I(exper^2)  -0.000882   0.000431 -2.045618 0.0414105 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.668704   Adj. R2: 0.136959
## F-test (1st stage), educ: stat = 87.7    , p < 2.2e-16 , on 1 and 424 DoF.
##               Wu-Hausman: stat =  1.43731, p = 0.231246, on 1 and 423 DoF.

4 2SLS (Two-stage least squares): One endogenous, multiple instrument

4.1 `fixest::feols` 関数で推定

操作変数が複数ある場合，fixest::feols 関数では アウトカム ~ 外生変数 | 内生変数 ~ 操作変数1 + 操作変数2 の形式で指定する．

iv4 <- feols(log(wage) ~ exper + I(exper^2) | educ ~ motheduc + fatheduc, data = mroz, vcov = "hetero")

比較する．

library(modelsummary)
modelsummary(list(summary(iv3, stage = 1:2), summary(iv4, stage = 1:2)), 
             gof_omit = "AIC|BIC|RMSE", stars = TRUE)

	iv: First stage: educ	iv: Second stage	iv: First stage: educ	iv: Second stage
(Intercept)	9.887***	-0.061	9.103***	0.048
	(0.376)	(0.458)	(0.424)	(0.430)
fatheduc	0.271***		0.190***
	(0.029)		(0.032)
exper	0.047	0.044**	0.045	0.044**
	(0.043)	(0.016)	(0.042)	(0.016)
I(exper^2)	-0.001	-0.001*	-0.001	-0.001*
	(0.001)	(0.000)	(0.001)	(0.000)
fit_educ		0.070+		0.061+
		(0.036)		(0.033)
motheduc			0.158***
			(0.035)
Num.Obs.	428	428	428	428
R2	0.176	0.143	0.211	0.136
R2 Adj.	0.170	0.137	0.204	0.130
Std.Errors	Heteroskedasticity-robust	Heteroskedasticity-robust	Heteroskedasticity-robust	Heteroskedasticity-robust
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

複数の操作変数がある場合は過剰識別検定を行うことができる．

fitstat(iv4, "sargan")

## Sargan: stat = 0.378071, p = 0.538637, on 1 DoF.

p = 0.5386 であるため，「すべての操作変数は外生」という帰無仮説を棄却することはできない．言い換えれば，操作変数のいずれかが invalid という証拠はこの検定からは得られていない．

ただし，操作変数の外生性を直接的に検証しているわけではない点に要注意．

5 Simulation: OLS vs. IV

中級者向け．

除外変数バイアスがある場合に操作変数を利用しないと推定値にバイアスが生じることをモンテカルロシミュレーションによって確認する．

5.1 Basic

「Step 1. データ生成 → Step 2. 推定 → Step 3. 推定値を記録」を100回繰り返す
記録する推定値は，OLS の内生変数の係数，IV の内生変数の係数，内生変数と操作変数の相関係数，First stage の F 値
- 推定値は coefs オブジェクト（matrix 型）に記録．Iteration のたびに推定値の行を追加．
- First stage の F 値は fitstat(iv_iter, "ivf") で出力される．$ivf1$stat で F 値を抽出
Data generating process は以下の通り

\[ z_i = \begin{cases} 1 \quad \mbox{if } \ \eta_i > 0 \\ 0 \quad \mbox{if } \ \eta_i \le 0 \end{cases}, \quad \eta_i \sim N(0, 1) \]

\[ x_i = \begin{cases} 1 \quad \mbox{if } \ 0.5 (\mbox{omitted var.})_i + 5 z_i + \xi_i > 0 \\ 0 \quad \mbox{if } \ 0.5 (\mbox{omitted var.})_i + 5 z_i + \xi_i \le 0 \end{cases}, \quad (\mbox{omitted var.})_i \sim N(0, 1), \quad \xi_i \sim N(0, 1) \]

\[ y_i = 5 x_i + 3 (\mbox{omitted var.})_i + \zeta_i, \quad \zeta_i \sim N(0, 1) \]

n <- 1000
coefs <- NULL  # 推定値を格納するオブジェクト
for (iter in 1:100) { 
  # Step 1. Generate data 
  set.seed(iter)  # 再現性のため
  sample_iter <- tibble::tibble(
    omitted = rnorm(n),  # omitted variable 
    inst = 1*(rnorm(n) > 0),  # instrument, same as [rbinom(n, 1, .5)] 
    endog = 1*(0.5 * omitted + 5 * inst + rnorm(n) > 0),  # endogenous variable 
    outcome = 5 * endog + 3 * omitted + rnorm(n)  # true coeff. = 5
    )
  # Step 2. Estimate 
  ols_iter <- feols(outcome ~ endog, data = sample_iter)
  iv_iter <- feols(outcome ~ 1 | endog ~ inst, data = sample_iter)
  # Step 3. Keep track of estimated values 
  coefs <- rbind(coefs, c(ols = as.numeric(ols_iter$coefficients[2]), 
                          iv = as.numeric(iv_iter$coefficients[2]),
                          cor_xz = cor(sample_iter$inst, sample_iter$endog),
                          iv_f = fitstat(iv_iter, "ivf")$ivf1$stat))
}
summary(coefs[, "cor_xz"])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5321  0.5617  0.5750  0.5761  0.5903  0.6211

summary(coefs[, "iv_f"])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   394.1   460.0   492.8   498.3   533.7   626.8

ggplot() + 
  geom_histogram(aes(x = coefs[, "ols"], linetype = "OLS"), alpha = 0.5) + 
  geom_histogram(aes(x = coefs[, "iv"], linetype = "IV"), alpha = 0.3, fill = "blue") + 
  geom_vline(xintercept = 5, color = "red") +  # true value 
  xlab("Coefficient")

操作変数を使わないと系統的なバイアスが生じる．操作変数を使うと平均的には概ね正しく推定される．

5.2 With defier

上の例ではサンプルの 100% が complier であることを仮定している（単調性の仮定を満たす）．

たとえば，サンプルの 30% が defier の場合は次の通り．

\[ x_{i: \ 30 \% \ of \ the \ sample} = \begin{cases} 1 \quad \mbox{if } \ 0.5 (\mbox{omitted var.})_i + 5 z_i + \xi_i < 0 \\ 0 \quad \mbox{if } \ 0.5 (\mbox{omitted var.})_i + 5 z_i + \xi_i \ge 0 \end{cases} \]

n <- 1000
coefs <- NULL
for (iter in 1:100) {
  set.seed(iter)
  sample_iter <- tibble::tibble(
    id = 1:n, 
    omitted = rnorm(n), 
    inst = 1*(rnorm(n) > 0), 
    endog = 1*(0.5 * omitted + 5 * inst + rnorm(n) > 0) 
    ) %>% 
    mutate(
      endog = ifelse(id %in% 1:(n*.3), 1*(0.5 * omitted + 5 * inst + rnorm(n) < 0), endog),  # defier (30%) 
      outcome = 5 * endog + 3 * omitted + rnorm(n)
    )
  ols_iter <- feols(outcome ~ endog, data = sample_iter)
  iv_iter <- feols(outcome ~ 1 | endog ~ inst, data = sample_iter)
  coefs <- rbind(coefs, c(ols = as.numeric(ols_iter$coefficients[2]), 
                          iv = as.numeric(iv_iter$coefficients[2]),
                          cor_xz = cor(sample_iter$inst, sample_iter$endog),
                          iv_f = fitstat(iv_iter, "ivf")$ivf1$stat))
}
summary(coefs[, "cor_xz"])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1176  0.1855  0.2009  0.2030  0.2235  0.2936

summary(coefs[, "iv_f"])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.99   35.55   41.97   43.84   52.49   94.15

ggplot() + 
  geom_histogram(aes(x = coefs[, "ols"], linetype = "OLS"), alpha = 0.5) + 
  geom_histogram(aes(x = coefs[, "iv"], linetype = "IV"), alpha = 0.3, fill = "blue") + 
  geom_vline(xintercept = 5, color = "red") +  # true value 
  xlab("Coefficient")

F 値が 10 を超えていても推定値にはかなりのばらつきがある．ましてや weak instrument だったら，まともな推定を期待するのは無理があるだろう．

Data Analysis Using Statistical Packages: Instrumental Variables

Sho Kuroda / 黒田翔

2024年4月 (Last update: 2024年5月)

1 Data

2 IV estimation: One endogenous, one instrument

2.1 Estimation by hand

3 2SLS (Two-stage least squares): One endogenous, one instrument

3.1 Estimation by hand

3.1.1 First stage

3.1.2 Second stage

3.2 `fixest::feols` 関数で推定

4 2SLS (Two-stage least squares): One endogenous, multiple instrument

4.1 `fixest::feols` 関数で推定

5 Simulation: OLS vs. IV

5.1 Basic

5.2 With defier

Data Analysis Using Statistical Packages: Instrumental Variables

Sho Kuroda / 黒田翔

2024年4月 (Last update: 2024年5月)

1 Data

2 IV estimation: One endogenous, one instrument

2.1 Estimation by hand

3 2SLS (Two-stage least squares): One endogenous, one instrument

3.1 Estimation by hand

3.1.1 First stage

3.1.2 Second stage

3.2 fixest::feols 関数で推定

4 2SLS (Two-stage least squares): One endogenous, multiple instrument

4.1 fixest::feols 関数で推定

5 Simulation: OLS vs. IV

5.1 Basic

5.2 With defier

3.2 `fixest::feols` 関数で推定

4.1 `fixest::feols` 関数で推定