I used the "systemfit" function in R to estimate a 2SLS model as it allows to specify the first and second stage separately which is important for my estimation. I need robust standard errors but I am having problems replicating the robust standard errors obtained when estimating the model using "iv_robust" or "ivreg" in combination with "coeftest" and vcov type "HC1".
I tried to calculate the robust standard errors manually as follows:
firststagehfsystem <- H ~ HL1
secondstageivsystem <- Index ~ H
system = list(firststagehfsystem, secondstageivsystem)
inst <- ~ HL1
ivhfsystem = systemfit(system, method = "2SLS", inst = inst, data = DataControls)
coef <- ivhfsystem$eq[[2]]$coefficients
vcov <- ivhfsystem$eq[[2]]$coefCov
se <- sqrt(diag(vcov))
residuals <- ivhfsystem$eq[[2]]$residuals
n <- nrow(DataControls)
k <- length(coef)
X <- model.matrix(ivhfsystem$eq[[2]])
hc1 <- (n/(n-k))*ivhfsystem$eq[[2]]$residuals^2
vce_hc1 <- solve(t(X) %*% X) %*% (t(X) %*% (hc1*diag(n)) %*% X) %*% solve(t(X) %*% X)
sqrt(diag(vce_hc1))
Can anyone point out my mistake? I already googled it for days and could not find a solution. Thank you!
X
represent the regressors when it should be the instruments in the 2nd-to-last line? $\endgroup$