Suponemos variables de media cero para simplificar. Entonces, como la muestra es i.i.d, es ergódica, y las medias muestrales convergen a los valores teóricos esperados. También debemos suponer que todos los momentos, individuales y conjuntos que aparecerán, existen y son finitos. Dados éstos, tenemos
$$S_n \equiv \hat \sigma_{xy} = \frac 1n \sum_{i=1}^n x_iy_i \to_p E(XY) = \sigma_{xy}$$
También,
$$E\left(\frac 1n \sum_{i=1}^n x_iy_i\right) = E(XY)$$
para muestras finitas.
Consideremos ahora la varianza de esta suma, escribiendo para mayor claridad $x_iy_i = w_i$ ,
$$\text{Var}\left(\frac 1n \sum_{i=1}^n w_i\right) = \frac 1 {n^2}\cdot \left[E\left(\sum w_i\right)^2 + n(n-1)\sum E(w_i)E(w_j) - \left (E\left[\sum w_i\right]\right)^2\right]$$
$$=\frac 1 {n^2} \cdot\left [ nE(w_i^2) + n(n-1)\sigma_{xy}^2- n^2 \sigma_{xy}^2 \right]= \frac 1 {n} \cdot\left [ E(x_i^2y_i^2) - \sigma_{xy}^2 \right]$$
Entonces tenemos que la expresión normalizada
$$\frac{S_n - E(S_n)}{\sqrt{\text{Var}(S_n)}} = \frac {\sqrt{n}(\hat \sigma_{xy} - \sigma_{xy})}{\sqrt{E(x_i^2y_i^2) - \sigma_{xy}^2}} \to_p N(0,1)$$
debido a la i.i.d. / ergodicidad de la muestra, y a los momentos existentes y finitos. De ello se deduce que
$$\sqrt{n}(\hat \sigma_{xy} - \sigma_{xy}) \approx N\left(0, E(x_i^2y_i^2) - \sigma_{xy}^2\right) $$