Variance of a sample covariance for normal variables

Question

In normally distributed populations, sample variance $s^2$ follows chi-squared distribution and the variance of this estimator is expressed by:

$$\text{Var}(s^2) = \frac{2\sigma^4}{n-1}$$

I am wondering what is the generalisation of this result to covariances.

In particular, I am looking at vector $\mathbf{X}=[X_1, X_2]^T$ with distribution $N(\mathbf{0}, \mathbf{\Sigma})$. Given $\mathbf{Q}$ is sample covariance $$\mathbf{Q} = \frac{1}{n-1}\sum_{i=1}^n \mathbf{x}_i\mathbf{x}_i^T$$ what is the distribution of $\mathbf{Q}$ and $\text{Var}(\mathbf{Q})$?

Answer 1

The OP is interested in Var(sample covariances) in a bivariate Normal world.

• You know the solution for: Var(sample variances) (main diagonal), so ...
• All that is needed is the solution for: Var(sample covariance)

Then there is no need for any matrix notation whatsoever, and if I understand correctly, the question reduces to:

Question: Let $(X,Y)$ be bivariate Normal. Given a sample of size $n$, namely $\big( (X_1, Y_1), \dots, (X_n, Y_n) \big)$, find $\text{Var}(W)$, where the sample covariance $W$ is defined as:

$$W = \frac 1 {n-1} \sum_{i=1}^n (X_i - \bar X)(Y_i-\bar Y). $$

Answer: This problem can be solved much more generally. We will obtain a general solution to $\text{Var}(W)$ for any distribution whose moments exist, and then solve for the Normal as a special case. The modus operandi for solving such problems is to work with power sum notation, which in our bivariate world is of form $s_{r,t}$, namely:

$$s_{r,t}=\sum _{i=1}^n X_i^r Y_i^t$$

In particular, sample covariance $W$ can be written in power sum notation as:

We seek $\text{Var}(W)$. Since the variance operator is the $2^\text{nd}$ Central Moment of $W$, we can find the variance using the mathStatica (for Mathematica) package function :

where:

• $\mu _{r,s}$ denotes the product central moment:

$$\mu _{r,s}=E\left[(X-E[X]]^r (Y-E[Y])^s\right]$$

For example, $\mu_{1,1} = \text{Cov}(X,Y)$, $\mu_{2,0}= \text{Var}(X)$ and $\mu_{0,2}= \text{Var}(Y)$.

This is a general solution for any bivariate distribution whose moments exist.

Special case of bivariate Normal

In the case of a bivariate Normal with variance-covariance:

$$\Sigma = \left( \begin{array}{cc} \sigma _{11}^2 & \rho \sigma _{11} \sigma _{22} \\ \rho \sigma _{11} \sigma _{22} & \sigma _{22}^2 \\ \end{array} \right)$$

... the various central moments $\mu_{i,j}$ are:

Substituting in these values into the general solution sol yields $\text{Var}(W)$ in the bivariate Normal case as:

I am not sure if this is the same as the OP's posted solution of: $\frac{1}{n-1}(\sigma_{ij}^2 + \sigma_{ii}\sigma_{jj})$, as the notation he intended for $\sigma_{ij}$ is not tied down and appears inconsistent with that used by the OP for the diagonal cases, which the OP states correctly as: $\text{Var}(s^2) = \frac{2\sigma^4}{n-1}$

In summary, for the Normal case, the Variance of sample covariances is:

$$ \text{Var} \left( \begin{array}{cc} \frac 1 {n-1} \sum_{i=1}^n (X_i - \bar X)^2 & \frac 1 {n-1} \sum_{i=1}^n (X_i - \bar X)(Y_i-\bar Y) \\ \frac 1 {n-1} \sum_{i=1}^n (X_i - \bar X)(Y_i-\bar Y) & \frac 1 {n-1} \sum_{i=1}^n (Y_i-\bar Y)^2 \\ \end{array} \right)$$

$$= \frac{1}{n-1}\left( \begin{array}{cc} 2 \sigma _{11}^4 & \left(1+\rho ^2\right) \sigma _{11}^2 \sigma _{22}^2 \\ \left(1+\rho ^2\right) \sigma _{11}^2 \sigma _{22}^2 & 2 \sigma _{22}^4 \\ \end{array} \right)$$

Answer 2

After following the suggestion from Mark Stone I looked up Wishart distribution and estimation of covariance matrices and here's a quick summary.

For a random, normally distributed $p$ element vector with a covariance matrix $\Sigma$ the quantity: $$ \sum_{i=1}^n \mathbf{X}_i\mathbf{X}_i^T \sim W_p(\Sigma, n-1) $$ where $W_p$ is a Wishart distribution. The formula for a variance of the sample covariance is then given by:

$$ \text{Var}(Q_{ij}) = \frac{1}{n-1}(\sigma_{ij}^2 + \sigma_{ii}\sigma_{jj}), $$ where $\sigma_{ij}$ are components of covariance matrix $\Sigma$.