# Exercise

1. Replicate the simulation from the right panel of my plot above:
1. Set the seed to `4321` and generate `n = 100` uniform draws for `x`
2. Set `y` equal to `0.2 + 0.9 * x + error` where `error` is a vector of independent, mean-zero normal errors with standard deviation `sqrt(2 * x)`.
3. Replicate my plot and check that yours matches it.
2. Using `x` and `y`, replicate my regression and F-test results from above.
3. 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` and `y`. Check that your results match those of `lm_robust()`.

# Solution

## Part 1

``````library(tidyverse)
set.seed(4321)

n <- 100
x <- runif(n)
error <- rnorm(n, mean = 0, sd = sqrt(2 * x))
intercept <- 0.2
slope <- 0.9
y <- intercept + slope * x + error

tibble(x, y) |>
ggplot(aes(x, y)) +
geom_smooth(method = 'lm') +
geom_point() ``````

## Part 2

``````library(estimatr)
library(car)
library(broom)
library(modelsummary)

reg_classical <- lm_robust(y ~ x, se_type = 'classical')

reg_robust <- lm_robust(y ~ x, se_type = 'HC0')

modelsummary(list(Classical = reg_classical, Robust = reg_robust),
fmt = 2,
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

``````reg <- lm(y ~ x)
uhat <- residuals(reg)
x_demeaned <- x - mean(x)
n <- length(uhat)

# Classical
sigma_sq_hat <- sum(uhat^2) / (n - 2) # two estimated parameters
var_classical <- sigma_sq_hat / sum(x_demeaned^2)
SE_classical <- sqrt(var_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
var_HC0 <- sum(uhat^2 * x_demeaned^2) / (sum(x_demeaned^2)^2)
SE_HC0 <- sqrt(var_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 ``````