Heteroscedasticity consistent (HC) standard error analysis and interaction effects in an OLS

I have made a model with several variables, and 8 of them interact with a dummy to find interaction effects. These are added stepwise, resulting in three models. Now, through a Breusch-Pagan test I have determined that my model contains heteroscedasticity.

I have determined that I will use HC4 (Cribari-Neto, 2004) for the HC standard error analysis, but am unsure how to continue. Thus, my question is:

Do I do the coeftest of the HC standard errors on the base model, or also on the models containing interaction effects?

Intuitively I would say that heteroscedasticity concerns the variables themselves and not necessarily their interaction effect, which means I would only do it on the base model. However, I have not been able to find any literature on this topic or situation.

Kind regards, DB

• There is in general no reason to suppose that heteroskedasticity would only relate to the variables. Consider the classical example of increased variability of earnings ($y$) given increased schooling ($X$). You may interact $X$ with a gender dummy, and it seems quite conceivable that the differential return to schooling for men vs women would also exhibit heteroskedasticity. Jul 29 '19 at 11:37

The heteroscedasticity-robust standard error take into account the whole regressor matrix $$X$$, in the case of HC4 this is the following "sandwich" formula
$$HC4 = (X'X)^{-1} X' \text{diag} \Big[ \frac{e_i^2}{(1-h_{ii})^{\delta_i}}\Big] X(X'X)^{-1}$$ with $$\delta_i = \text{min} \Big\{4,\frac{nh_{ii}}{1-p}\Big\}$$
Since an interaction term is always also a column of $$X$$, a heteroscedasticity-robust will take into account also the variance component of the interaction term. Everything else would make no sense, since a model without interaction term is a completely different model with different standard error for the coefficients but also different (higher) residual standard error.