TRUE_BETA <- 0.08Practice session: OLS, endogeneity and consistency
Applied Economics Issues, Master 2 DADEE
Duration: 30 minutes for Parts 0–3, plus about 20 minutes for the extension (Part 4).
Objective
The question investigated in this notebook is the following: does OLS recover the true value of an estimate, and what happens when the sample grows?
Note
In this notebook, we simulate the data, so we know the true return to one year of schooling:
In real life, we never do.
How to work. Double click on the RProj file in the folder that contains the present HTML file. Run the chunks one by one. Fill in the blanks (___) and answer the questions in the text below each block.
The data (one csv per case and sample size, e.g. data/ovb_n1000.csv):
| Variable | Meaning |
|---|---|
educ |
years of schooling, as observed by the econometrician |
lnwage |
log wage |
The four cases are called exog, ovb, irrelevant and measerr. You do not know yet what is behind them.
| Time | Part |
|---|---|
| 0–3 min | 0. Look at the data |
| 3–10 min | 1. OLS by hand, exog case |
| 10–17 min | 2. The four cases at n = 1 000, and what is behind them |
| 17–27 min | 3. Convergence: what happens when n grows? |
| 27–30 min | Wrap-up |
| 30–50 min | 4. Extension: an instrument to the rescue |
We will need the following packages:
library(tidyverse)
library(here)Let us generate the data, in the data/ folder, using write_all_datasets() from the "dgp.R" script (do not look at this script file yet):
source("dgp.R")
purrr::map(
.x = c(100, 1000, 10000, 100000),
.f = ~ write_all_datasets(ns = .x, dir = "data", seed_base = 2031)
)[[1]]
[1] "exog_n100.csv" "exog_n1000.csv" "exog_n10000.csv"
[4] "exog_n100000.csv" "irrelevant_n100.csv" "irrelevant_n1000.csv"
[7] "irrelevant_n10000.csv" "irrelevant_n100000.csv" "measerr_n100.csv"
[10] "measerr_n1000.csv" "measerr_n10000.csv" "measerr_n100000.csv"
[13] "ovb_n100.csv" "ovb_n1000.csv" "ovb_n10000.csv"
[16] "ovb_n100000.csv"
[[2]]
[1] "exog_n100.csv" "exog_n1000.csv" "exog_n10000.csv"
[4] "exog_n100000.csv" "irrelevant_n100.csv" "irrelevant_n1000.csv"
[7] "irrelevant_n10000.csv" "irrelevant_n100000.csv" "measerr_n100.csv"
[10] "measerr_n1000.csv" "measerr_n10000.csv" "measerr_n100000.csv"
[13] "ovb_n100.csv" "ovb_n1000.csv" "ovb_n10000.csv"
[16] "ovb_n100000.csv"
[[3]]
[1] "exog_n100.csv" "exog_n1000.csv" "exog_n10000.csv"
[4] "exog_n100000.csv" "irrelevant_n100.csv" "irrelevant_n1000.csv"
[7] "irrelevant_n10000.csv" "irrelevant_n100000.csv" "measerr_n100.csv"
[10] "measerr_n1000.csv" "measerr_n10000.csv" "measerr_n100000.csv"
[13] "ovb_n100.csv" "ovb_n1000.csv" "ovb_n10000.csv"
[16] "ovb_n100000.csv"
[[4]]
[1] "exog_n100.csv" "exog_n1000.csv" "exog_n10000.csv"
[4] "exog_n100000.csv" "irrelevant_n100.csv" "irrelevant_n1000.csv"
[7] "irrelevant_n10000.csv" "irrelevant_n100000.csv" "measerr_n100.csv"
[10] "measerr_n1000.csv" "measerr_n10000.csv" "measerr_n100000.csv"
[13] "ovb_n100.csv" "ovb_n1000.csv" "ovb_n10000.csv"
[16] "ovb_n100000.csv" Part 0 (3 min): look at the data
files <- list.files("data", pattern = "\\.csv$")
files [1] "exog_n100.csv" "exog_n1000.csv" "exog_n10000.csv"
[4] "exog_n100000.csv" "irrelevant_n100.csv" "irrelevant_n1000.csv"
[7] "irrelevant_n10000.csv" "irrelevant_n100000.csv" "measerr_n100.csv"
[10] "measerr_n1000.csv" "measerr_n10000.csv" "measerr_n100000.csv"
[13] "ovb_n100.csv" "ovb_n1000.csv" "ovb_n10000.csv"
[16] "ovb_n100000.csv" d <- read_csv(here("data", "exog_n1000.csv"))
head(d)# A tibble: 6 × 5
educ lnwage z u_oracle educ_true_oracle
<dbl> <dbl> <dbl> <dbl> <dbl>
1 10.9 2.48 0 0.368 10.9
2 11.2 2.41 1 0.958 11.2
3 9.16 1.58 0 -0.635 9.16
4 12.1 2.47 0 0.267 12.1
5 10.5 2.47 0 0.298 10.5
6 12.5 2.78 1 -0.613 12.5 summary(d) educ lnwage z u_oracle
Min. : 7.614 Min. :1.385 Min. :0.000 Min. :-2.82059
1st Qu.:11.437 1st Qu.:2.267 1st Qu.:0.000 1st Qu.:-0.65502
Median :12.536 Median :2.506 Median :1.000 Median :-0.02244
Mean :12.494 Mean :2.498 Mean :0.507 Mean :-0.02667
3rd Qu.:13.583 3rd Qu.:2.725 3rd Qu.:1.000 3rd Qu.: 0.61335
Max. :18.387 Max. :3.601 Max. :1.000 Max. : 3.03141
educ_true_oracle
Min. : 7.614
1st Qu.:11.437
Median :12.536
Mean :12.494
3rd Qu.:13.583
Max. :18.387
Q0
How many observations does d have? Which variables would you use in a wage regression?
Your answer:
Part 1 (7 min): OLS in the exog case, with lm() and by hand
1.1 Estimate \(\ln w_i = \alpha + \beta\,\text{educ}_i + \varepsilon_i\) with lm().
fit <- lm(lnwage ~ educ, data = d)
summary(fit)
Call:
lm(formula = lnwage ~ educ, data = d)
Residuals:
Min 1Q Median 3Q Max
-1.09971 -0.22130 0.00121 0.20929 1.17042
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.50903 0.08038 18.77 <2e-16 ***
educ 0.07916 0.00638 12.41 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.326 on 998 degrees of freedom
Multiple R-squared: 0.1336, Adjusted R-squared: 0.1327
F-statistic: 153.9 on 1 and 998 DF, p-value: < 2.2e-16confint(fit) 2.5 % 97.5 %
(Intercept) 1.35129362 1.66676842
educ 0.06663561 0.091676541.2 Recompute the coefficients with matrix algebra, \(\hat{\boldsymbol\beta}=(\mathbf X'\mathbf X)^{-1}\mathbf X'\mathbf y\).
X <- cbind(1, d$educ) # n x 2: a constant and educ
y <- d$lnwage # n x 1
b_hat <- ___ # TODO
b_hat
all.equal(as.numeric(b_hat), as.numeric(coef(fit))) # should be TRUEError in parse(text = input): <text>:4:11: unexpected input
3:
4: b_hat <- __
^1.3 Compare with the truth.
beta_hat <- b_hat[2]Error:
! object 'b_hat' not foundbeta_hat - TRUE_BETA # the "sampling error"Error:
! object 'beta_hat' not found
Q1
- Is 0.08 inside the 95% confidence interval? (b) Is the estimate exactly 0.08? Why (not)?
Your answer:
Part 2 (7 min): the other cases, n = 1 000
2.1 Below is a helper function that reads a file and returns the OLS estimate on educ with its 95% CI. Complete the two lines marked TODO (b is the estimate, s its standard error).
ols_educ <- function(file) {
d <- read_csv(here("data", file))
fit <- lm(lnwage ~ educ, data = d)
cf <- summary(fit)$coefficients
b <- cf["educ", "Estimate"]
s <- cf["educ", "Std. Error"]
tibble(
file = file,
n = nrow(d),
beta_hat = b,
se = s,
lo = ___, # TODO: lower bound of the 95% CI (use b and s)
hi = ___ # TODO: upper bound of the 95% CI
)
}Error in parse(text = input): <text>:12:17: unexpected input
11: se = s,
12: lo = __
^We consider the four simulated datasets with 1,000 observations, load them and regress wages on education.
res_1000 <- map(
.x = c("exog_n1000.csv", "ovb_n1000.csv",
"irrelevant_n1000.csv", "measerr_n1000.csv"),
.f = ols_educ
) |>
bind_rows()Error:
! object 'ols_educ' not foundThen, we add a binary variable covers_truth indicating whether the true effect is within the bounds of the 95% CI.
res_1000 <- res_1000 |>
mutate(
covers_truth = lo <= TRUE_BETA & TRUE_BETA <= hi
)Error:
! object 'res_1000' not foundres_1000Error:
! object 'res_1000' not found
Q2a
In which cases is the estimate far from 0.08? Is it bad luck (sampling error) or something systematic? Can you tell from n = 1 000 only?
Your answer:
2.2 Lifting the veil. Here is what is behind each case. All the random terms below are drawn independently of each other.
- \(u\sim\mathcal N(0,1)\) is ability (unobserved in real life);
- \(z\sim\text{Bernoulli}(1/2)\) is a scholarship offer, a fair coin toss (it shifts schooling in every case; it is used as an instrument in the bonus);
- \(e_x\sim\mathcal N(0,1.5^2)\) is everything else that drives schooling;
- \(v\sim\mathcal N(0,0.3^2)\) is a wage shock;
- \(m\sim\mathcal N(0,1)\) is measurement noise (used only in
measerr).
Schooling and wages are generated as \[ \begin{aligned} \text{educ}_{\text{true}} & = 12 + a\,u + z + e_x\\ \ln w & = 1.5 + \beta\,\text{educ}_{\text{true}} + \gamma\,u + v,\qquad \beta=0.08. \end{aligned} \]
The observed value of schooling is \[\text{educ} = \text{educ}_{\text{true}} + \text{err}\times m,\] where \(\text{err}=1\) only in the measerr case, and \(0\) otherwise.
| Case | \(a\) (ability \(\to\) schooling) | \(\gamma\) (ability \(\to\) wage) | \(\text{err}\) | description |
|---|---|---|---|---|
exog |
0 | 0.10 | 0 | ability affects wages, but is unrelated to schooling |
ovb |
1 | 0.10 | 0 | ability affects both: omitted variable bias |
irrelevant |
1 | 0 | 0 | ability is related to schooling, but does not affect wages |
measerr |
0 | 0 | 1 | schooling is observed with noise |
Hence, this gives the four following data generating processes:
| Case | \(\text{educ}_{\text{true}}\) | Wage equation | Observed educ |
Regression error \(\varepsilon\) |
|---|---|---|---|---|
exog |
\(12+z+e_x\) | \(1.5+0.08\,\text{educ}_{\text{true}}+0.10\,u+v\) | \(\text{educ}_{\text{true}}\) | \(0.10\,u+v\) |
ovb |
\(12+u+z+e_x\) | \(1.5+0.08\,\text{educ}_{\text{true}}+0.10\,u+v\) | \(\text{educ}_{\text{true}}\) | \(0.10\,u+v\) |
irrelevant |
\(12+u+z+e_x\) | \(1.5+0.08\,\text{educ}_{\text{true}}+v\) | \(\text{educ}_{\text{true}}\) | \(v\) |
measerr |
\(12+z+e_x\) | \(1.5+0.08\,\text{educ}_{\text{true}}+v\) | \(\text{educ}_{\text{true}}+m\) | \(v-0.08\,m\) |
Variance of schooling. First, note that \(\text{Var}(u)=\text{Var}(m)=1\), \(\text{Var}(z)=p(1-p)=\frac{1}{2} \times \frac{1}{2}=0.25\), and \(\text{Var}(e_x)=1.5^2=2.25\). Because all terms are independent, variances add:
exog: \(\text{Var}(\text{educ})=0.25+2.25=2.5\)ovbandirrelevant: \(\text{Var}(\text{educ})=1+0.25+2.25=3.5\)measerr: \(\text{Var}(\text{educ}_{\text{true}})=0.25+2.25=2.5\) and \(\text{Var}(\text{educ})=2.5+1=3.5\)
What OLS should converge to (read this carefully, no need to derive it yourself).
In every case, we regress \(\ln w\) on a constant and \(\text{educ}\): \[\ln w_i = \alpha + \beta\,\text{educ}_i + \varepsilon_i,\qquad i=1,\dots,n,\] and estimate \((\alpha,\beta)\) by OLS. The true value of the slope in the data-generating process is \(\beta=0.08\), and \(\hat\beta\) is our estimate of it. The error term \(\varepsilon_i\) of this regression is the last column of the table above.
The probability limit of the slope is \[\text{plim}\,\hat\beta=0.08+\frac{\text{Cov}(\text{educ},\varepsilon)}{\text{Var}(\text{educ})}.\]
pred_plim <- c(
exog = TRUE_BETA,
ovb = TRUE_BETA + 0.10 * 1 / (1 + 0.25 + 2.25),
irrelevant = TRUE_BETA,
measerr = TRUE_BETA * (0.25 + 2.25) / (0.25 + 2.25 + 1)
)
pred_plim exog ovb irrelevant measerr
0.08000000 0.10857143 0.08000000 0.05714286 Part 3 (10 min): convergence, what happens when n grows?
3.1 Estimate the OLS coefficient for each case and each sample size.
The data folder contains, for each data generating process, 4 different samples with increasing number of observations: \(100\), \(1,000\), \(10,000\), \(100,000\).
all_cases <- expand_grid(
case = c("exog", "ovb", "irrelevant", "measerr"),
n = c(100, 1000, 10000, 100000)
) |>
mutate(
file = sprintf("%s_n%d.csv", case, n)
)
results <- map(
.x = all_cases$file,
.f = ols_educ
) |>
bind_rows() |>
left_join(all_cases, by = c("file", "n")) |>
relocate(case, .before = file) |>
select(-file)Error:
! object 'ols_educ' not foundresultsError:
! object 'results' not found3.2 Let us plot the estimate and its 95% CI against \(n\). Solid line: truth. Dashed red line: the predicted plim from 2.2.
ggplot(
data = results
) +
geom_hline(yintercept = TRUE_BETA, colour = "black") +
geom_hline(
data = tibble(pred_plim = pred_plim, case = names(pred_plim)),
mapping = aes(yintercept = pred_plim),
colour = "red", linetype = "dashed"
) +
geom_pointrange(
mapping = aes(x = n, y = beta_hat, ymin = lo, ymax = hi)
) +
facet_wrap(facets = vars(case)) +
scale_x_log10() +
labs(
x = "Sample size, log scale",
y = "OLS estimate of beta",
caption = "Black line: true value. Dashed red line: p-lim."
)Error:
! object 'results' not found
Q3a–c
- What happens to the standard error when \(n\) grows?
- In the
exogcase, where does \(\hat\beta\) converge? And in the other cases? - In the
ovbcase, is 0.08 inside the 95% CI for n = 100? For n = 100 000? What does this tell you about “collecting more data”?
Your answer:
3.3 How many standard errors away from the truth are we?
results <- results |>
mutate(t_vs_truth = ___) # TODO
results |> select(case, n, t_vs_truth)Error in parse(text = input): <text>:2:24: unexpected input
1: results <- results |>
2: mutate(t_vs_truth = __
^
Q3d
How does this ratio behave with \(n\) in the exog case? In the ovb case?
Your answer:
Wrap-up (3 min)
Discussion
- Why does a larger sample not fix endogeneity? What does a narrow CI tell you, and what does it not tell you?
- In real data, we have no
u_oraclecolumn and no true \(\beta\). How could we detect (or get around) the problem?
Answer:
Part 4 (extension, about 20 min): an instrument to the rescue
In ovb and measerr, OLS is inconsistent: educ is correlated with the error term. To recover \(\beta=0.08\) we need variation in schooling that is unrelated to the error. In our data, we have exactly that: z, a scholarship offer randomly assigned by a coin toss. It is a valid instrument for educ because
- relevance: it shifts schooling (in the DGP, by 1 year on average);
- exogeneity: it is randomly assigned, hence independent of ability \(u\), of the wage shock \(v\) and of the noise \(m\);
- exclusion: it does not enter the wage equation directly. It affects wages only through schooling.
With \(\mathbf Z=[\mathbf 1\;\; z]\) (an \(n\times 2\) matrix, in the same way as \(\mathbf X=[\mathbf 1\;\;\text{educ}]\)), the IV estimator is \[\hat{\boldsymbol\beta}_{IV}=(\mathbf Z'\mathbf X)^{-1}\mathbf Z'\mathbf y .\]
4.1 First stage and reduced form. We work with the ovb case and \(n=10{,}000\).
d_iv <- read_csv(here("data", "ovb_n10000.csv"), show_col_types = FALSE)
first_stage <- lm(educ ~ z, data = d_iv) # schooling on the instrument
reduced_form <- lm(lnwage ~ z, data = d_iv) # wage on the instrument
summary(first_stage)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.9961953 0.02564575 467.76545 0.000000e+00
z 0.9971754 0.03623597 27.51894 9.202734e-161summary(first_stage)$fstatistic[["value"]] # first-stage F statistic[1] 757.292coef(reduced_form)["z"] z
0.07643758
Q4a
What does the first-stage coefficient on z tell you? Is the instrument strong (rule of thumb: \(F>10\))? What does the reduced-form coefficient measure?
Your answer:
4.2 The Wald ratio. The offer raises schooling by \(\hat\pi_1\) years (first stage) and raises wages by \(\hat\rho\) (reduced form). Each additional year of schooling therefore raises wages by \(\hat\rho/\hat\pi_1\).
wald <- ___ # TODO: reduced form / first stage
waldError in parse(text = input): <text>:1:10: unexpected input
1: wald <- __
^4.3 IV by hand. Compute \((\mathbf Z'\mathbf X)^{-1}\mathbf Z'\mathbf y\) and compare it with OLS, with the Wald ratio and with the truth.
X <- cbind(1, d_iv$educ)
Z <- cbind(1, d_iv$z)
y <- d_iv$lnwage
b_iv <- ___ # TODO: (Z'X)^{-1} Z'y
b_ols <- solve(t(X) %*% X) %*% t(X) %*% y
tibble(
estimator = c("OLS", "IV", "Wald ratio", "Truth"),
beta = as.numeric(c(b_ols[2], b_iv[2], wald, TRUE_BETA))
)Error in parse(text = input): <text>:5:11: unexpected input
4:
5: b_iv <- __
^
Q4b
Compare OLS, IV and the truth. Why do IV and the Wald ratio give the same number?
Your answer:
4.4 Standard errors. For a just-identified model with homoskedastic errors, \[\widehat{\text{Var}}(\hat{\boldsymbol\beta}_{IV})=\hat\sigma^2\,(\mathbf Z'\mathbf X)^{-1}\mathbf Z'\mathbf Z\,(\mathbf X'\mathbf Z)^{-1},\qquad \hat\sigma^2=\frac1n\sum_i \hat\varepsilon_i^2,\] where the residuals \(\hat\varepsilon_i=y_i-\mathbf x_i'\hat{\boldsymbol\beta}_{IV}\) are computed with the observed educ, not with fitted values from the first stage.
n <- nrow(X)
res <- y - X %*% b_iv
sigma2 <- sum(res^2) / n
ZXinv <- solve(t(Z) %*% X)
V_iv <- ___ # TODO
se_iv <- sqrt(diag(V_iv))
tibble(
estimator = c("OLS", "IV"),
estimate = c(b_ols[2], b_iv[2]),
se = c(summary(lm(lnwage ~ educ, data = d_iv))$coefficients["educ", "Std. Error"],
se_iv[2])
) |>
mutate(lo = estimate - 1.96 * se, hi = estimate + 1.96 * se)Error in parse(text = input): <text>:6:11: unexpected input
5:
6: V_iv <- __
^
Q4c
Does the 95% CI of each estimator contain 0.08? Is IV more or less precise than OLS, and why?
Your answer:
4.5 IV at different sample sizes. The helper below returns the IV estimate, its 95% CI and the first-stage \(F\) for one file. We apply it to ovb and measerr, for every sample size.
iv_educ <- function(file) {
d <- read_csv(here("data", file), show_col_types = FALSE)
X <- cbind(1, d$educ)
Z <- cbind(1, d$z)
y <- d$lnwage
n <- nrow(X)
ZXinv <- solve(t(Z) %*% X)
b <- ZXinv %*% t(Z) %*% y
res <- y - X %*% b
sigma2 <- sum(res^2) / n
V <- sigma2 * ZXinv %*% t(Z) %*% Z %*% t(ZXinv)
tibble(
file = file,
n = n,
beta_iv = b[2],
se_iv = sqrt(V[2, 2]),
lo = beta_iv - 1.96 * se_iv,
hi = beta_iv + 1.96 * se_iv,
F_first = summary(lm(educ ~ z, data = d))$fstatistic[["value"]]
)
}
iv_cases <- all_cases |> filter(case %in% c("ovb", "measerr"))
iv_results <- map(iv_cases$file, iv_educ) |>
bind_rows() |>
left_join(iv_cases, by = c("file", "n")) |>
select(case, n, beta_iv, se_iv, lo, hi, F_first)
iv_results# A tibble: 8 × 7
case n beta_iv se_iv lo hi F_first
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ovb 100 0.167 0.224 -0.272 0.606 0.627
2 ovb 1000 0.0922 0.0178 0.0574 0.127 87.8
3 ovb 10000 0.0767 0.00639 0.0641 0.0892 757.
4 ovb 100000 0.0798 0.00202 0.0758 0.0837 7488.
5 measerr 100 0.222 0.101 0.0248 0.419 5.69
6 measerr 1000 0.0721 0.0177 0.0373 0.107 82.9
7 measerr 10000 0.0827 0.00639 0.0702 0.0952 744.
8 measerr 100000 0.0772 0.00196 0.0734 0.0811 7675. Now compare OLS and IV, as \(n\) grows:
compare <- bind_rows(
results |>
filter(case %in% c("ovb", "measerr")) |>
transmute(case, n, estimator = "OLS", estimate = beta_hat, lo, hi),
iv_results |>
transmute(case, n, estimator = "IV", estimate = beta_iv, lo, hi)
)Error:
! object 'results' not foundggplot(data = compare, mapping = aes(x = n, y = estimate, ymin = lo, ymax = hi, colour = estimator)) +
geom_hline(yintercept = TRUE_BETA) +
geom_pointrange(position = position_dodge(width = 0.1)) +
facet_wrap(facets = vars(case)) +
scale_x_log10() +
coord_cartesian(ylim = c(-0.1, 0.4)) +
labs(
x = "Sample size, log scale",
y = "Estimate of beta",
colour = NULL,
caption = "Black line: true value. Confidence intervals are clipped at the bounds of the y-axis."
)Error:
! object 'compare' not found
Q4d
- What do you see for \(n=100\) and \(n=1{,}000\)? Look at the first-stage \(F\).
- What happens for large \(n\)?
- Compare with OLS.
Your answer:
Discussion
Which of the three assumptions (relevance, exogeneity, exclusion) can we check with the data we have? Which can we not, and what do we rely on instead?
Your answer: