I want to prove that $V$ is an unbiased estimator of the covariance matrix $$(X'X)^{-1}(X'DX)(X'X)^{-1},$$ where $D=diag(\sigma^2,...,\sigma^2)=E(ee'|X)$ in a linear model.
$$V = \frac{n}{n-k}(X'X)^{-1}\left(\sum_{i=1}^n{X_iX'_i\hat{e}^2_i}\right)(X'X)^{-1}$$
To do so, I first find the conditional expectation of $V$.
$$E[V|X]= \frac{n}{n-k}(X'X)^{-1}\left(\sum_{i=1}^n{X_iX'_iE(\hat{e}^2_i}|X)\right)(X'X)^{-1}$$
However, I am not sure how to proceed from here.
I do not think your claim is correct. You are analyzing what MacKinnon and White (1985) refer to as the heteroskedasticity-robust variance estimator $HC_1$. I will argue that, instead, the estimator they call $HC_2$ is unbiased when in fact no heteroskedasticity is present, while $HC_1$ only is in a special case. Here is more detail:
Notice that the residuals satisfy $\hat e=My$ for $M=I-H$ and hat matrix $$H=\{h_{ij}\}_{i,j=1,\ldots,n}=\{x_i'(X'X)^{-1}x_j\}_{i,j=1,\ldots,n}=X(X'X)^{-1}X'.$$ Hence, $$ \hat e_i=(1-h_{ii})y_i-\sum_{j\neq i}h_{ij}y_j, $$ so that $$\hat e_i^2=\left[(1-h_{ii})y_i-\sum_{j\neq i}h_{ij}y_j\right]^2$$ In conditional expectation and under random sampling (i.i.d.) with conditional homoskedasticity - assumptions you do not state but that I need - $$ \begin{align*} E(\hat e_i^2|X)&=Var(\hat e_i|X) \end{align*}$$ in view of $E(\hat e|X)=ME(y|X)=MX\beta=0$. Now, again by the i.i.d. assumption,
\begin{align*} Var(\hat e_i|X)&=(1-h_{ii})^2Var(y_i|X)+\sum_{j\neq i}h_{ij}^2Var(y_j|X)\\ &=(1-h_{ii})^2\sigma^2+\sum_{j\neq i}h_{ij}^2\sigma^2, \end{align*}
as covariances are zero and all conditional variances are assumed to be the same. Next, notice that $$ \sum_{j=1}^nh_{ij}^2=h_{ii},$$ as, by symmetry of $H$, \begin{align*} \sum_{j=1}^nh_{ij}^2&=\sum_{j=1}^nh_{ij}h_{ji}\\ &=\sum_{j=1}^nx_i'(X'X)^{-1}x_jx_j'(X'X)^{-1}x_i\\ &=x_i'(X'X)^{-1}\sum_{j=1}^nx_jx_j'(X'X)^{-1}x_i\\ &=x_i'(X'X)^{-1}(X'X)(X'X)^{-1}x_i\\ &=x_i'(X'X)^{-1}x_i\\ &=h_{ii}\\ \end{align*}
Hence, by multiplying out and rearranging, \begin{align*} E(\hat e_i^2|X)&=(1-2h_{ii}+h_{ii}^2)\sigma^2+\sum_{j\neq i}h_{ij}^2\sigma^2\\ &=(1-h_{ii})\sigma^2-h_{ii}\sigma^2+\sum_{j=1}^nh_{ij}^2\sigma^2\\ &=(1-h_{ii})\sigma^2, \end{align*} Hence, under homoskedasticity, unbiasedness would call for adjusting squared residuals by $1-h_{ii}$, not multiplying by $n/(n-k)$.
Only in a "balanced design" (note that $\sum_ih_{ii}=k$, so that $k/n$ is the average value of the $h_{ii}$) in which $h_{ii}=k/n$ for all $i$ would the two coincide, as we would then have $$ \frac{1}{1-h_{ii}}=\frac{1}{1-k/n}=\frac{n}{n-k} $$ So all in all, for $$ HC_2=(X'X)^{-1}\left(\sum_{i=1}^n{x_ix'_i\frac{\hat{e}^2_i}{1-h_{ii}}}\right)(X'X)^{-1} $$ we obtain \begin{align*} E(HC_2|X)&=(X'X)^{-1}\left(\sum_{i=1}^nx_ix'_i\sigma^2\right)(X'X)^{-1}\\&=\sigma^2(X'X)^{-1} \end{align*} Hence, $HC_2$'s conditional expected value equals the conditional variance of $\hat\beta$ under homoskedasticity, which is different from the unconditional variance, see here.
Here is a little simulation to illustrate. I take the regressor to be fixed in repeated samples here to avoid the distinction between conditional and unconditional variance discussed above.
I find that the bias of $HC_2$ tends to be an order of magnitude smaller, although both biases are very small for the designs considered here.
library(sandwich)
n <- 50
x <- rnorm(n, sd=3)
x <- rt(n, df=2)
x <- runif(n, -10, 10)
sigma <- .2
mc.function <- function(n){
u <- rnorm(n, sd=sigma)
y <- 2*x + u
limo <- lm(y~x-1)
return(c(vcovHC(limo, "HC1"), vcovHC(limo, "HC2")))
}
true.cond.var <- sigma^2/sum(x^2)
vcovs <- replicate(1000, mc.function(n))
(bias <- rowMeans(vcovs-true.cond.var))
abs(bias[1]) > abs(bias[2])
