Consistency of heteroskedasticity-robust standard errors Is the following statement true or not?
When applying OLS with model $y=a+bx+u$,
the heteroskedasticity-robust standard errors are consistent because $\hat u_i^2$ (the squared OLS residual) is a consistent estimator of $E(u_i^2|x_i)$ for each $i$.
 A: This looks like the standard assertion that one can find in any econometrics textbook that deals with the subject, but it really isn't.  The question asks about individual consistency of every OLS residual separately:
$$\hat u^2_i \xrightarrow{p}??\;\; E(u^2_i\mid x)$$
Certainly not, since 
$$\hat u^2_i = (y_i - \hat a -\hat b x_i)^2 \xrightarrow{p} (y_i - a -bx_i)^2 = u^2_i$$
given the usual assumptions of the linear regression model. What it does hold and is true, is that
$$\frac 1n\sum_{i=1}^n \hat u^2_i\mathbf x_i\mathbf x_i' \xrightarrow{p} E(u^2_i\mathbf x_i\mathbf x_i')$$
...under the additional assumption that the fourth moments of the regressors exist and are finite (an assumption that is not needed in the homoskedasticity benchmark case for the other OLS properties).
For an alternative expression, defining the diagonal matrix $\mathbf {\hat U} = \text{diag}(\hat u^2_1,...,\hat u^2_n)$, and denoting $\mathbf U$ the corresponding matrix containing the error terms, we have equivalently
$$\frac 1n \mathbf X\mathbf {\hat U}\mathbf X' \xrightarrow{p} E(\mathbf X\mathbf { U}\mathbf X')$$
The reason why the additional assumption is needed, is easily shown for the case of one regressor (and no constant term), where  $\mathbf x_i\mathbf x_i' =x_i^2$, but it of course generalizes:
$$\hat u^2_i = (y_i  -\hat b x_i)^2 = [u_i - (\hat b-b)x_i]^2  = u_i^2 - 2(\hat b-b)u_ix_i + (\hat b-b)^2ux_i^2$$
Therefore
$$\frac 1n\sum_{i=1}^n \hat u^2_ix_i^2 = \frac 1n\sum_{i=1}^n  u^2_ix_i^2 - 2(\hat b-b)\frac 1n\sum_{i=1}^n  u_ix_i^3 +(\hat b-b)^2\frac 1n\sum_{i=1}^n  u^2_ix_i^4$$
\hat b-b
The fourth sample moment of the regressor has appeared. By assuming that the population counterpart exists and is finite, implying that the same holds for the third moment also, we guarantee that the second and the third sums of the right-hand side converge to something finite, and then the elements themselves to zero due to consistency of the OLS estimator ($\hat b\rightarrow b$), leaving us with
$$\text{plim} \frac 1n\sum_{i=1}^n \hat u^2_ix_i^2 = \text{plim} \frac 1n\sum_{i=1}^n  u^2_ix_i^2 =E(u_i^2x_i^2)$$
