# Manually calculate robust standard errors of 2SLS regression

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!

• Without having digested your code fully, doesn't X represent the regressors when it should be the instruments in the 2nd-to-last line? Jan 18 at 9:05

The systemfit package reverse depends the sandwich package which says the following about their vcovHC program:

Long & Ervin (2000) conduct a simulation study of HC estimators (HC0 to HC3) in the linear regression model, recommending to use HC3 which is thus the default in vcovHC.

Your by-hand approach is applying the HC1 type which will be different. This is really a question about replicating robust standard errors on a generic model fit. If you inspect the functions in sandwich, the type= argument is passed to meatHC as a switch with the following

 HC1 = {
omega <- function(residuals, diaghat, df) residuals^2 *
length(residuals)/df

HC3 = {
omega <- function(residuals, diaghat, df) residuals^2/(1 -
diaghat)^2


Where diaghat <- try(hatvalues(x), silent = TRUE) is the vector of hat values from the influence function. See the referenced "Belsley, D. A., Kuh, E. and Welsch, R. E. (1980). Regression Diagnostics": if you want more detail on the calculation. R's internals use a C-subroutine that doesn't include explicit calculation instructions.

• The question is not about which standard errors to use but how to replicate them. Can you show how to manually calculate the HC3 type standard errors? When I apply coeftest using vcov type HC3 or HC1 after estimating the model via systemfit I get different standard errors compared to when I use the iv_robust or ivreg function with vcov type HC3 or HC1. How can I replicate the standard errors I get when using iv_robust or ivreg? Jan 18 at 7:31
• @JJ1214 It seemed apropos to point out why the default is different than you have calculated, and perhaps why that would be the case. Anyway, I edited the answer as much as probably is necessary. Jan 18 at 17:59