Given an i.i.d. sample $(x_i,y_i)$ of size $n$ from a bivariate distribution $(x,y)$, I'm trying to prove that the sample covariance $$\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})$$ is an unbiasted estimator of the true population covariance $\text{cov}(x,y)$.
I started by writing the sample covariance in matrix form as $$\frac{1}{n-1}(\mathbf{Dx})'\mathbf{Dy}=\frac{1}{n-1}\mathbf{x'D'Dy}=\frac{1}{n-1}\mathbf{x'DDy}=\frac{1}{n-1}\mathbf{x'Dy}$$ where $\mathbf{D}$ is the symmetric and idempotent "centering matrix" $$\mathbf{D}=\mathbf{I}_n-\frac{1}{n}\iota\iota'$$ that puts the data in deviation from their sample mean and $\mathbf{x}$ and $\mathbf{y}$ are the $x_i$ and $y_i$ data points collected into $n\times 1$ column vectors. Taking expectations, and using "Tr" to denote trace, I get \begin{eqnarray} E\left[\frac{1}{n-1}\mathbf{x'Dy}\right]&=&\frac{1}{n-1}E[\mathbf{x'Dy}] \nonumber\\ &=&\frac{1}{n-1}E[\text{Tr}(\mathbf{x'Dy})]\nonumber\\ &=&\frac{1}{n-1}E[\text{Tr}(\mathbf{Dyx'})]\nonumber\\ &=&\frac{1}{n-1}\text{Tr}(E[\mathbf{Dyx'}])\nonumber\\ &=&\frac{1}{n-1}\text{Tr}(\mathbf{D}E[\mathbf{yx'}])\nonumber\\ &=&\frac{1}{n-1}\text{Tr}(\mathbf{D}(\text{cov}(x,y)+\mu_x\mu_y)\mathbf{I})\nonumber\\ &=&\frac{1}{n-1}\text{Tr}(\mathbf{DI})(\text{cov}(x,y)+\mu_x\mu_y)\nonumber\\ &=&\frac{1}{n-1}\text{Tr}(\mathbf{D})(\text{cov}(x,y)+\mu_x\mu_y)\nonumber\\ &=&\frac{1}{n-1}(n-1)(\text{cov}(x,y)+\mu_x\mu_y)\nonumber\\ &=&\text{cov}(x,y)+\mu_x\mu_y\nonumber \end{eqnarray}
But this clearly can't be right, because the expectation is supposed to come out just to be $\text{cov}(x,y)$. And yet I can't spot my mistake. Can someone point out where I've made an error?
The only mistake I could anticipate is that it's wrong to claim $$E[\mathbf{yx'}]=(\text{cov}(x,y)+\mu_x\mu_y)\mathbf{I}$$ but every entry in the matrix $\mathbf{yx'}$ has the form $y_ix_j$, and since the sample is i.i.d., $E[y_ix_j]=E[yx]=\text{cov}(x,y)+\mu_x\mu_y$.
