manual implementation of HC0-HC4 robust standard errors for reference

HC0_4_manual.R

library(car); library(sandwich)
data(auto)
# %%
fit = lm(price ~ mpg + weight, data = auto)
X = model.matrix(fit); n = nrow(X); k = ncol(X)
e = resid(fit)
A = crossprod(X)
H = X %*% solve(A) %*% t(X); h_ii = diag(H)
# cla =         (t(e) %*% e)/(n-k) %*% solve(t(X) %*% X)
CLA = as.numeric(crossprod(e)/(n-k)) * solve(A)
# HC0 =           solve(t(X) %*% X) %*% t(X) %*% diag(e)                %*% X %*% solve(t(X) %*% X)
B = wcrossprod(X, w = e^2)
HC0 = solve(A) %*% B %*% solve(A)
# HC1 = (n/(n-k))*solve(t(X) %*% X) %*% t(X) %*% diag(e^2)              %*% X %*% solve(t(X) %*% X)
HC1 = (n/(n-k)) * solve(A) %*% B %*% solve(A)
# HC2 =           solve(t(X) %*% X) %*% t(X) %*% diag(e^2/(1-h_ii))     %*% X %*% solve(t(X) %*% X)
B2 = wcrossprod(X, w = e^2/(1-h_ii))
HC2 = solve(A) %*% B2 %*% solve(A)
# HC3 =           solve(t(X) %*% X) %*% t(X) %*% diag((e^2/(1-h_ii)^2)) %*% X %*% solve(t(X) %*% X)
B3 = wcrossprod(X, w = (e^2/(1-h_ii)^2))
HC3 = solve(A) %*% B3 %*% solve(A)
# %%
semake = \(x) x |> diag() |> sqrt() |> round(2)
Map(semake, list(CLA, HC0, HC1, HC2, HC3))
# %%
Map(semake,
  list(vcov(fit), vcovHC(fit, "HC0") , vcovHC(fit, "HC1") , vcovHC(fit, "HC2"), vcovHC(fit, "HC3")))
# %%