Second moment of draws from a multivariate normal covariance matrix

Question

Suppose I have an $N\times N$ covariance matrix that describes a multivariate normal joint distribution.

Now take 100,000 draws of the covariance matrix. I measure the variance of values for each draw -- call this the "cross sectional dispersion". If I average the result of cross-sectional dispersion across draws what is the analytical expression for the mean and variance of this dispersion?

Here is some code if you'd like to toy with this or test your idea:

# create 3x3 covariance matrix
m1 <- c( .1 , .15 , 0 )
S1 <- matrix( c( .1 , .05 , .02 , .05 , .1 , .03 , .02 , .03 , .1 ), nrow = 3 )

# draw 10,000 times from covariance matrix
set.seed(1000)   
library(MASS)
draws <- MASS::mvrnorm(10000,mu=m1,Sigma=S1)

# measure the "cross-sectional" dispersion at each row
rowDispersion <- apply( draws , 1 , var )

# measure the mean cross-sectional dispersion
mean( rowDispersion )

# measure the variance of cross-sectional dispersion
var( rowDispersion )

Answer 1

I don't think you need to create the term "cross-sectional dispersion". From your code, it seems that you are seeking the first two moments of the sample variance when the sample values are not iid draws from a normal distribution but rather are correlated draws (a vector drawn from a multivariate normal).

You have $(X_1, ..., X_N) \sim MVN$ with known mean vector and covariance matrix. You calculate the sample variance, $S^2 = \sum (X_i-\bar{X})^2/(N-1)$ and want to know the expected value and variance of $S^2$.

You can write $S^2$ as a function of $\sum X_i^2$ and $(\sum X_i)^2$; the latter can be written as $\sum_{i,j} X_i X_j$.

The expected value will be relatively simple to obtain using the means, variances, and covariances.

The variance will be a bit tedious and will require you to look at the 3rd and 4th moments of the normal distribution (see the table in the Moments section of the Wikipedia page on the normal distribution), but shouldn't be too deep.

---

Let me spell out how to calculate the expected value:

$$\mathbb{E}(S^2) = \mathbb{E}\left[\frac{1}{N-1} \sum (X_i - \bar{X})^2\right] = \frac{1}{N-1}\left\{ \mathbb{E}\left[\sum X_i^2\right] - \frac{1}{N}\mathbb{E}\left[\left(\sum X_i\right)^2\right]\right\}$$

And note that $\mathbb{E}\left[\sum X_i^2\right] = \sum \mathbb{E}\left(X_i^2\right) = \sum (\mu_i^2 + \sigma_i^2)$

And then $\mathbb{E}\left[\left( \sum X_i \right)^2\right] = \mathbb{E}\left[ \sum_{i,j} X_i X_j \right] = \sum_{i,j} \mathbb{E}\left( X_i X_j \right) = \sum_{i,j} (\mu_i \mu_j + \sigma_{ij}) = \left(\sum_i \mu_i\right)^2 + \sum_{i,j} \sigma_{i,j}$

Thus we obtain $$\mathbb{E}(S^2) = \frac{1}{N-1} \left\{ \sum_i \; \mu_i^2 + \sum_i \; \sigma_i^2 - \frac{1}{N}\left(\sum_i \; \mu_i\right)^2 - \frac{1}{N} \sum_{i,j} \; \sigma_{ij} \right\}$$

Here $\sigma_i^2$ is the $i$th diagonal element of the covariance matrix and $\sigma_{ij}$ is the $(i,j)$th element, so $\sigma^2_i = \sigma_{ii}$.

---

For your example code, I calculate an expected value of 0.0725:

(sum(m1^2) + sum(diag(S1)) - sum(m1)^2/3 - sum(S1)/3)/2

Note that you could also write this as:

n <- nrow(S1) 
var(m1) + sum(diag(S1) - mean(S1)) / (n-1)