# Variance of a sample covariance for normal variables

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})$?

• The Wishart distribution is your friend stat.pitt.edu/sungkyu/course/2221Fall13/lec2.pdf . Jun 24, 2017 at 17:47
• Yes, thanks. I found it the moment I formulated my question, but I am still not able to figure out the variance estimate for off-diagonal components. Jun 24, 2017 at 17:55
• I think Proposition 8.3. (II) of projecteuclid.org/download/pdf_1/euclid.lnms/1196285114 is what you want. Jun 24, 2017 at 18:49
• Even though the parent distribution has zero means, are you certain that it is correct or appropriate (depending on defn) to remove the sample means from the defn of sample (co)variance, when your intention is to calculate moments of the sample (co)variance. Are you imposing that the sample means are zero? Or are they being calculated? Jun 24, 2017 at 21:23
• I think you're right, @wolfies. Since you providendan answer to a more general case with means and sample means, shall I rewrite the question slightly? Jun 25, 2017 at 21:56

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)$$

• Thanks, my solution, as you say, didn't have consistent notation as I dropped squares when defining $\Sigma$. Your answer takes a more diligent approach to this notation but I think the final formulas are equivalent. Thanks also for mathStatica examples too. Jun 25, 2017 at 21:52

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$.

• I don't think so. See my comment on the question. Jun 24, 2017 at 18:50
• Yes, I've just done a numerical experiment and there's still something wrong. Just reading it now. Jun 24, 2017 at 18:54
• @MarkL.Stone, I think my answer is actually correct. I found a mistake in my numerical experiment and it agrees exactly with the theory now. I was taking twice as many samples as I thought. This is, of course, not a proof, but an indication. Jun 24, 2017 at 19:36
• Your formula for $Q_{ij}$ doesn't agree with your $Var(s^2)$ above in one dimension. Jun 24, 2017 at 20:16
• I see the problem. I used inconsistently the notation for $\sigma$. In 1D variance is $\sigma^2$ but in 2D it's $\sigma$ Let me think how to fix it. Jun 24, 2017 at 20:39