# Heteroskedasticly Consistent Estimators for Var-Cov Matrix, Large Sample OLS Regression

I have a cross-sectional data sample of nearly 40,000 observations and tests for heteroskedasticity fail to reject the assumption of homoskedasticity. However, it seems common practice to report heteroskedastic errors as well. From what I've been researching, it seems like there are multiple ways to estimate heteroskedasticly-consistent standard errors (e.g., Sandwich Estimators, etc.) Does anyone know the most appropriate heteroskedastically consistent estimator for a sample size such as mine?

Addition: I don't know which software you are using, but in R you can use the sandwich package, and in STATA you can just add ,robust after your regression, i.e. reg y x, robust.