Unbiased estimator of square of covariance

Question

Let $X$ and $Y$ be two correlated random variables. The covariance is:

$\mu_{1,1} = E[(X-E[X])(Y-E[Y])]$

And the unbiased estimator:

$\hat{\mu}_{1,1} = \frac{1}{n-1} \sum_{i=1}^n (x_i-\bar{x})(y_i-\bar{y})$

According to some numerical experiments I've run, $\hat{\mu}_{1,1}^2$ is not an unbiased estimator of $\mu_{1,1}^2$. So, apparently I need the "multivariate polyache" (see p.11-12 of this book), but I can't find an implementation.

Any suggestions?

Thank you!

Answer 1

I'll just go ahead and answer my own question since I've found the answer in the Mathematica documentation:

An unbiased estimator of the square of covariance is the multivariate polyache $h_{\left \{ \left \{ 1,1 \right \},\left \{ 1,1 \right \} \right \} }$

Where

$S_{r,s} = \sum_{i=1}^n x_i^r y_i^s$

MultiPolyH[P] is a Mathematica function I've written (details bellow)

---

Just to have it documented, here is a more general answer:

For $X_1,X_2,\dots,X_d$ random variables, an unbiased estimator of the product $\prod_{j=1}^m\mu_{p_{j,1},p_{j,2},...,p_{j,d}}$, where

$\mu_{p_{j,1},p_{j,2},...,p_{j,d}} = E \left[ (X_1-\bar{X}_1)^{p_{j,1}}(X_2-\bar{X}_2)^{p_{j,2}}\dots(X_d-\bar{X}_d)^{p_{j,d}} \right]$

is the multivariate polyache $h_{ P }$, where

$P = \left \{ \left \{ p_{1,1},p_{1,2},...,p_{1,d} \right \} ,\left \{ p_{2,1},p_{2,2},...,p_{2,d}\right \}, \dots \left \{ p_{m,1},p_{m,2},...,p_{m,d} \right \}\right \}$

The following function implemented in Mathematica computes $h_P$.

MultiPolyH[P_] := 
FullSimplify[
    MomentConvert[
        Product[CentralMoment[P[[i]]], {i, 1, Length[P]}], 
        "UnbiasedEstimator", "PowerSymmetricPolynomial"]
    /.{PowerSymmetricPolynomial[0] -> n}, 
    Assumptions -> Element[n, Integers]]

For example, the product of variances $\mu_{2,0}\mu_{0,2}$ can be estimated by:

Note that the Mathematica function also works for $d=1$ and $m=1$. For example the third central moment: