# Distribution of Variance of draws from Multivariate Normal Distribution

Let the vector $$\boldsymbol x$$ be a draw of $$n$$ values from a multivariate normal distribution with zero mean. $$\boldsymbol x \sim \mathcal{N}(\boldsymbol 0, \Sigma)$$ It may be assumed that $$\Sigma$$ is standardised.

Let the scalar $$y$$ be the variance of $$\boldsymbol x$$: $$y = var(\boldsymbol x)$$ Then, what is the distribution of $$y$$, given $$n$$ and $$\Sigma$$? $$y \sim ?(n, \Sigma)$$

Context: an autoregressive process (AR1) generates time series where values close in time are more correlated than values further apart. Generating $$n$$ such values is equivalent to sampling from a multivariate normal distribution with a particular covariance structure, namely $$\Sigma_{ij} = \rho^{|t_i - t_j|}$$. What is the variance $$y$$ of the values of these time series, and how is that distributed when multiple time series are generated?

• Just to check, $x$ and $\mu$ are vectors and $\Sigma$ is a matrix, but you want the scalar $y$. This is going to be messy, especially if the elements of $\mu$ are not all equal, or if $\Sigma$ is not exchangeable Jan 18, 2022 at 18:18
• Correct. I will edit to make vectors bold, for clarity. Also, it may be assumed that $\mu=\overrightarrow{0}$. Jan 18, 2022 at 18:43
• The answer to the question you asked will be truly messy: it's a generalization of a non-central Chi distribution where the components can have different variances. It gets a little less messy if you don't take the square root: that would be the variance. If you also assume $\Sigma$ is constant diagonal, you will have a multiple of the familiar Chi-square distribution. If you don't, you will have have a Gamma mixture.
– whuber
Jan 18, 2022 at 19:14
• If $\Sigma$ is diagonal (and standardised), then all $x_i$ are i.i.d., right? That would indeed make things easier. Likewise, a perfectly correlated $\Sigma$ also yields a straight forward answer. I am curious about the area in between. A function which (efficiently) simulates the answer is also fine, or some close approximation. I have changed s.d. to var per your suggestion. Jan 18, 2022 at 20:04

You mentioned simulation in one of your comments. This is easily done and the R code below does so for the particular case where the vector has $$10$$ elements, $$\mu=0$$ and $$\Sigma$$ is the identity matrix, but this can easily be changed for other $$\mu$$ and $$\Sigma$$. Note that the variance formula used is the $$\frac1{n-1}$$ version.

As whuber commented, in this particular case you get a chi-square distribution with $$10-1=9$$ degrees of freedom, with values multiplied by $$\frac19$$, so this distribution of variances has mean $$1$$ and variance $$\frac29$$.

library(MASS)
library(matrixStats)
set.seed(2022)
mu <- rep(0,10)
Sigma <- diag(10)
sims <- rowVars(mvrnorm(10^5, mu, Sigma))
mean(sims)
# 0.9987011
var(sims)
# 0.2225065


To illustrate that this particular case does have that distribution, here are the cdf and density from the simulation (black) compared to the theoretical curves (red). The overlap is so good that it is difficult to distinguish the curves

plot.ecdf(sims) plot(density(sims)) • That's great, thanks. But what about non-trivial $\Sigma$? I understand the direct simulation (black line) is possible in that case, and your code easily allows for it, but is the resulting distribution (red line) obtainable via some more efficient method? For example, a converging sum over terms or some iterative method, if not an analytical expression? Perhaps it is simply not possible to be more efficient than pure simulation, that's a fine answer too. Jan 19, 2022 at 9:52