Manually calculate robust standard errors of 2SLS regression

Question

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!

Answer 1

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.