library(tidyverse)
set.seed(4321)
<- 100
n <- runif(n)
x <- rnorm(n, mean = 0, sd = sqrt(2 * x))
error <- 0.2
intercept <- 0.9
slope <- intercept + slope * x + error
y
tibble(x, y) |>
ggplot(aes(x, y)) +
geom_smooth(method = 'lm') +
geom_point()
Robust Standard Errors - Solutions
Exercise
- Replicate the simulation from the right panel of my plot above:
- Set the seed to
4321
and generaten = 100
uniform draws forx
- Set
y
equal to0.2 + 0.9 * x + error
whereerror
is a vector of independent, mean-zero normal errors with standard deviationsqrt(2 * x)
. - Replicate my plot and check that yours matches it.
- Set the seed to
- Using
x
andy
, replicate my regression and F-test results from above. - Use the formulas from earlier in this lecture to compute the “classical” and “HC0” standard errors for the regression slope “by hand” based on
x
andy
. Check that your results match those oflm_robust()
.
Solution
Part 1
Part 2
library(estimatr)
library(car)
library(broom)
library(modelsummary)
<- lm_robust(y ~ x, se_type = 'classical')
reg_classical
<- lm_robust(y ~ x, se_type = 'HC0')
reg_robust
modelsummary(list(Classical = reg_classical, Robust = reg_robust),
fmt = 2,
gof_omit = 'R2 Adj.|AIC|BIC')
Classical | Robust | |
---|---|---|
(Intercept) | 0.34 | 0.34 |
(0.22) | (0.17) | |
x | 0.78 | 0.78 |
(0.38) | (0.40) | |
Num.Obs. | 100 | 100 |
R2 | 0.041 | 0.041 |
RMSE | 1.03 | 1.03 |
linearHypothesis(reg_classical, 'x = 0') |> tidy()
# A tibble: 1 × 8
term null.value estimate std.error statistic p.value df.residual df
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 x 0 0.777 0.379 4.21 0.0402 98 1
linearHypothesis(reg_robust, 'x = 0') |> tidy()
# A tibble: 1 × 8
term null.value estimate std.error statistic p.value df.residual df
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 x 0 0.777 0.403 3.72 0.0538 98 1
Part 3
<- lm(y ~ x)
reg <- residuals(reg)
uhat <- x - mean(x)
x_demeaned <- length(uhat)
n
# Classical
<- sum(uhat^2) / (n - 2) # two estimated parameters
sigma_sq_hat <- sigma_sq_hat / sum(x_demeaned^2)
var_classical <- sqrt(var_classical)
SE_classical
c(lm_robust = tidy(reg_classical) |>
filter(term == 'x') |>
pull(std.error),
by_hand = SE_classical)
lm_robust by_hand
0.378506 0.378506
# HC0
<- sum(uhat^2 * x_demeaned^2) / (sum(x_demeaned^2)^2)
var_HC0 <- sqrt(var_HC0)
SE_HC0
c(lm_robust = tidy(reg_robust) |>
filter(term == 'x') |>
pull(std.error), by_hand = SE_HC0)
lm_robust by_hand
0.4027311 0.4027311