Limiting Distribution Given Samples from Joint Distribution

Let $$(X_i, Y_i)$$, $$1 \leq i \leq n$$ be independent and identically distributed samples from a joint distribution $$F (x, y)$$. Suppose that $$E[X^4], E[Y^4] < \infty$$.

Now define $$\sigma_{XY} = E[(X - \mu_X)(Y - \mu_Y)]$$ and an estimator $$\hat{\sigma}_{XY} = \frac{1}{n} \sum \limits_{i = 1}^n (X_i - \bar{X})(Y_i - \bar{Y})$$. Find the limiting distribution for $$\sqrt{n} (\hat{\sigma}_{XY} - \sigma_{XY})$$.

I've shown that $$\hat{\sigma}_{XY}$$ converges to $$\sigma_{XY}$$ almost surely. Perhaps it's useless for this problem, but I'm having trouble translating that into a convergence in distribution statement as required. I've also had mind to try and maybe determine the limiting distribution for $$\sqrt{n} (\hat{\sigma^2}_{XY} - \sigma^2_{XY})$$ and then apply the Delta Method to transform the distribution, but to no success either. Any tips or tricks for how to solve this efficiently?

Note that $$\hat{\sigma}_{XY} = \frac{1}{n} \sum \limits_{i = 1}^n (X_i - \mu_X)(Y_i - \mu_Y) + R_n$$ where \begin{align} R_n &= (\mu_X - \bar{X}) \cdot\frac{1}{n} \sum_{i=1}^n(Y_i - \mu_Y) + (\mu_Y - \bar{Y}) \cdot\frac{1}{n} \sum_{i=1}^n(X_i - \mu_X) \\ & \;\;+ (\mu_X - \bar{X})\cdot(\mu_Y - \bar{Y}). \end{align}
It follows from CLT (for i.i.d. sequence with finite second moment) that $$\frac{1}{\sqrt{n}} \sum \limits_{i = 1}^n (X_i - \mu_X)(Y_i - \mu_Y) \stackrel{d}{\rightarrow} N(\sigma_{XY}, \,\cdots).$$ (I leave to you to write down the asymptotic variance "$$\cdots$$".)
So what remains to be shown is that $$\sqrt{n} R_n = o_p(1)$$. Take, for example, the first term in the expression of $$R_n$$, $$\sqrt{n} (\mu_X - \bar{X}) \cdot\frac{1}{n} \sum_{i=1}^n(Y_i - \mu_Y) = \underbrace{ (\mu_X - \bar{X}) }_{A} \, \cdot \, \underbrace{ \frac{1}{\sqrt{n}} \sum_{i=1}^n(Y_i - \mu_Y)}_{B}.$$
By LLN, $$A \in o_p(1)$$, and by CLT, $$B \in O_p(1)$$, so their product $$AB \in o_p(1)$$. Similarly, the other two terms in the expression of $$\sqrt{n} R_n$$ are also $$o_p(1)$$. This proves the claim.
• Wow, this is really clever. How did you know to express $\hat{\sigma}_{XY}$ in terms of the convergent term and $R_n$? The form of $R_n$ is extremely unintuitive to me. Oct 4, 2020 at 21:46
• Also, it's not clear to me what the asymptotic variance should be. It seems like I want $Var((X - \mu_X)(Y - \mu_Y))$, but I can't get around the mental idea of taking the variance of what appears to be covariance. Is this kosher just because we want the limiting distribution? Oct 4, 2020 at 21:48
• @rcmpgrc "...to express $\sigma_XY$ in terms of the convergent term and $R_n$..."?---Replacing the sample moment by the population moment, or vice versa, is the natural thing to do in most cases. Since sample moment converges to population moment when LLN holds, one can hope the remainder term is small, which in this case it is. Oct 4, 2020 at 22:12
• @rcmpgrc "...what the asymptotic variance should be...It seems like $Var\left((X - \mu_X)(Y - \mu_Y)\right)$..."---yes, $Var\left((X - \mu_X)(Y - \mu_Y)\right)$ exists in this case because the assumption that $X$ and $Y$ have finite fourth moments and the Cauchy-Schwarz inequality. So one can simply use the formula $Var(AB) = E[A^2 B^2] - E[AB]^2$, where $A = X - \mu_X$, $B = Y - \mu_Y$. (Here, $E[AB]^2 = \sigma_{XY}^2$.) Oct 4, 2020 at 22:15