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)

Can anyone point out my mistake? I already googled it for days and could not find a solution. Thank you!

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

1 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 * 

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

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.

  • $\begingroup$ 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? $\endgroup$
    – JJ1214
    Jan 18 at 7:31
  • $\begingroup$ @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. $\endgroup$
    – AdamO
    Jan 18 at 17:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.