Wishart Fluctuations and the Covariance Matrix CLT

Why Covariance Matrices Fluctuate Together

A sample covariance matrix is not a collection of unrelated estimates. Its entries are computed from the same observations, so their sampling fluctuations can be correlated. Under multivariate normal sampling, the centered sample covariance has an exact Wishart representation, and its large-sample fluctuations have a Gaussian limit. Knowing the full covariance of those fluctuations is essential whenever several entries are considered together.

Throughout, \(p\) is fixed as \(n\) grows, and \(X_1,\ldots,X_n\) are independent and identically distributed \(p\)-dimensional random vectors with distribution \(N_p(\mu,\Sigma)\), where \(\Sigma\) is positive definite. Define the sample mean and unbiased sample covariance by

\bar X=\frac{1}{n}\sum_{r=1}^n X_r, \qquad S=\frac{1}{n-1}\sum_{r=1}^n (X_r-\bar X)(X_r-\bar X)^\mathsf{T}.

We use column-wise vectorization: \(\operatorname{vec}(A)\) stacks the first column of \(A\), then the second, and so on. Thus the coordinate of \(\operatorname{vec}(S)\) corresponding to \(S_{ij}\) is indexed by \(i+(j-1)p\).

The Exact Wishart Representation

Definition: If \(Z_1,\ldots,Z_\nu\) are independent \(N_p(0,\Sigma)\) vectors, the random matrix \(\sum_{r=1}^{\nu} Z_rZ_r^\mathsf{T}\) has a Wishart distribution with scale matrix \(\Sigma\) and \(\nu\) degrees of freedom, written \(W_p(\Sigma,\nu)\).

For normal observations, centering by the sample mean changes the degrees of freedom from \(n\) to \(n-1\), but preserves the Wishart form.

Result: Under the assumptions above, \((n-1)S\sim W_p(\Sigma,n-1)\). Moreover, \(S\) is independent of \(\bar X\).

To see why, arrange the observations as the rows of an \(n\) by \(p\) matrix. Choose an orthogonal \(n\) by \(n\) matrix whose first column is \(n^{-1/2}\mathbf 1\), where \(\mathbf 1\) is the vector of ones. Multiplying the data matrix by the transpose of this orthogonal matrix transforms the first row into \(\sqrt n\,\bar X^\mathsf{T}\). The remaining \(n-1\) rows are orthogonal contrasts: their coefficients sum to zero, so the common mean \(\mu\) contributes nothing to them.

The transformed rows are jointly normal. Their cross-covariances are zero because the rows of the orthogonal matrix are orthogonal; jointly normal vectors with zero cross-covariance are independent. Consequently, the first transformed row is independent of the remaining rows, and each remaining row, transposed, has distribution \(N_p(0,\Sigma)\). Orthogonal transformation preserves the sum of centered outer products, so that sum is the sum of the outer products of these \(n-1\) independent normal vectors. This proves both the Wishart representation and the independence assertion.

The Covariance Matrix CLT

Write \(K_{p,p}\) for the commutation matrix, the \(p^2\) by \(p^2\) permutation matrix satisfying \(K_{p,p}\operatorname{vec}(A)=\operatorname{vec}(A^\mathsf{T})\) for every \(p\) by \(p\) matrix \(A\).

Theorem: Suppose \(X_1,\ldots,X_n\) are independent \(N_p(\mu,\Sigma)\) vectors, with fixed \(p\) and positive definite \(\Sigma\), and define \(S\) using divisor \(n-1\). Then

\sqrt n\,\operatorname{vec}(S-\Sigma) \xrightarrow{d} N_{p^2}\!\left(0,\,(I_{p^2}+K_{p,p})(\Sigma\otimes\Sigma)\right).

Equivalently, for any \(1\leq i,j,k,l\leq p\), the limiting covariance between the fluctuations of entries \(S_{ij}\) and \(S_{kl}\) is

\lim_{n\to\infty} n\,\operatorname{Cov}(S_{ij},S_{kl}) = \Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}.

The covariance formula follows from the fourth moments of a centered multivariate normal vector. For completeness, if \(Z\sim N_p(0,\Sigma)\), then

E[Z_iZ_jZ_kZ_l] = \Sigma_{ij}\Sigma_{kl} +\Sigma_{ik}\Sigma_{jl} +\Sigma_{il}\Sigma_{jk}.

One way to obtain this identity is to expand the moment-generating function \(E[\exp(t^\mathsf{T}Z)]=\exp(t^\mathsf{T}\Sigma t/2)\). The coefficient of \(t_it_jt_kt_l\), with repeated indices interpreted by differentiation, comes from pairing the four factors in the three possible ways; each pairing contributes the product of the corresponding two covariances. Subtracting \(E[Z_iZ_j]E[Z_kZ_l]=\Sigma_{ij}\Sigma_{kl}\) therefore gives

\operatorname{Cov}(Z_iZ_j,Z_kZ_l) = \Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}.

This identity verifies the covariance used in the theorem and also explains the entries in its matrix form. Here is a direct proof of the limit.

Use the Wishart sum.
By the exact representation, \((n-1)S=\sum_{r=1}^{n-1}Z_rZ_r^\mathsf{T}\), where the \(Z_r\) are independent \(N_p(0,\Sigma)\) vectors. Set \(m=n-1\) and \(W_r=\operatorname{vec}(Z_rZ_r^\mathsf{T}-\Sigma)\). Then \(E[W_r]=0\), and the fourth-moment calculation gives the covariance of \(W_r\).

Apply the multivariate central limit theorem.
The vectors \(W_r\) are independent and identically distributed with finite second moments. The multivariate iid central limit theorem therefore gives \(m^{-1/2}\sum_{r=1}^{m}W_r\xrightarrow{d}N_{p^2}(0,\Gamma)\), where \(\Gamma\) has entry \(\Gamma_{(ij),(kl)}=\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}\).

Account for the divisor.
Since \(S-\Sigma=m^{-1}\sum_{r=1}^{m}(Z_rZ_r^\mathsf{T}-\Sigma)\), multiplying its vectorization by \(\sqrt n\) multiplies the central-limit expression by \(\sqrt{n/m}\). As \(m=n-1\), \(\sqrt{n/m}\to1\), so this factor does not change the limit.

Identify the matrix covariance.
The \(((i,j),(k,l))\) entry of \((I_{p^2}+K_{p,p})(\Sigma\otimes\Sigma)\) is \(\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}\). Thus \(\Gamma=(I_{p^2}+K_{p,p})(\Sigma\otimes\Sigma)\), which proves the stated limit. \(\square\)

Worked Examples

Worked Example: The One-Dimensional Limit

Let \(p=1\), so \(X_r\sim N(\mu,\sigma^2)\), and \(S\) is the usual unbiased sample variance. The entry covariance formula has only one entry: set \(i=j=k=l=1\). Then

\lim_{n\to\infty}n\,\operatorname{Var}(S) = \sigma^2\sigma^2+\sigma^2\sigma^2 = 2\sigma^4.

Thus \(\sqrt n(S-\sigma^2)\xrightarrow{d}N(0,2\sigma^4)\). This agrees with the exact normal-theory variance \(\operatorname{Var}(S)=2\sigma^4/(n-1)\): multiplying by \(n\) gives \(2\sigma^4n/(n-1)\), which tends to \(2\sigma^4\).

Worked Example: All Entry Covariances in a Two-Variable Model

Suppose the population covariance matrix is

\Sigma= \begin{pmatrix} 4&1\\ 1&9 \end{pmatrix}.

This is positive definite because its leading diagonal entry is \(4>0\) and its determinant is \(4(9)-1^2=35>0\). For the distinct entries, use the order \((S_{11},S_{12},S_{22})\). The limiting variance of \(S_{11}\), on the \(\sqrt n\) scale, is \(2\Sigma_{11}^2=2(4^2)=32\). The limiting variance of \(S_{22}\) is \(2\Sigma_{22}^2=2(9^2)=162\). For the off-diagonal entry it is

\Sigma_{11}\Sigma_{22}+\Sigma_{12}^2 =4(9)+1^2=37.

The cross-covariance between \(S_{11}\) and \(S_{12}\) is \(2\Sigma_{11}\Sigma_{12}=2(4)(1)=8\). Between \(S_{22}\) and \(S_{12}\), it is \(2\Sigma_{22}\Sigma_{12}=2(9)(1)=18\). Finally, between \(S_{11}\) and \(S_{22}\), it is \(2\Sigma_{12}^2=2\). Therefore,

\sqrt n \begin{pmatrix} S_{11}-4\\ S_{12}-1\\ S_{22}-9 \end{pmatrix} \xrightarrow{d} N_3\!\left( \begin{pmatrix}0\\0\\0\end{pmatrix}, \begin{pmatrix} 32&8&2\\ 8&37&18\\ 2&18&162 \end{pmatrix} \right).

The nonzero off-diagonal terms show why separate standard errors do not describe the joint fluctuation of the covariance estimates.

Worked Example: A Linear Summary of the Covariance Matrix

Consider the sample covariance \(S_{12}\) when the two variables are uncorrelated, with

\Sigma= \begin{pmatrix} 2&0\\ 0&5 \end{pmatrix}.

Taking \(i=1,j=2,k=1,l=2\) gives

\lim_{n\to\infty}n\,\operatorname{Var}(S_{12}) = \Sigma_{11}\Sigma_{22}+\Sigma_{12}\Sigma_{21} = 2(5)+0=10.

Hence \(\sqrt n S_{12}\xrightarrow{d}N(0,10)\), since the population covariance is zero. For comparison, the limiting variances of \(S_{11}\) and \(S_{22}\) are \(2(2^2)=8\) and \(2(5^2)=50\), respectively. Their limiting covariances with \(S_{12}\) are zero: for example,

\lim_{n\to\infty}n\,\operatorname{Cov}(S_{11},S_{12}) =2\Sigma_{11}\Sigma_{12}=0.

The joint Gaussian limit therefore has independent distinct entries in this diagonal-\(\Sigma\) case, because its covariance matrix is diagonal. Zero covariance implies independence here because the limiting vector is jointly Gaussian.

Reading the Limit Correctly

The vectorized limit is generally a degenerate Gaussian distribution in \(p^2\) coordinates. This is not a defect: \(S\) is symmetric for every sample, so \(S_{ij}=S_{ji}\) exactly. The vector \(\operatorname{vec}(S-\Sigma)\) therefore lies in a lower-dimensional subspace, and its covariance matrix cannot be full rank. One may instead retain only the \(p(p+1)/2\) distinct entries, but their covariance still follows from the same entrywise formula.

The limit concerns fixed \(p\) and \(n\to\infty\). It does not by itself establish an approximation when the number of variables grows with the sample size. Also, the multiplier is \(\sqrt n\), although the exact Wishart degrees of freedom are \(n-1\); the proof shows why the two normalizations have the same limit.

Pitfall: treating the entries as independent. The incorrect reasoning is that each covariance entry is estimated from the same data but can nevertheless be assigned an independent Gaussian error. The correction is to calculate \(\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}\) for every pair of entries; these terms need not vanish.
Pitfall: using the wrong variance for an off-diagonal entry. The variance of \(S_{ij}\), for \(i\ne j\), is not generally \(2\Sigma_{ij}^2/n\). Substituting \(k=i,l=j\) gives the correct limiting variance \(\Sigma_{ii}\Sigma_{jj}+\Sigma_{ij}^2\).
Pitfall: overlooking symmetry in vector form. A covariance matrix for all \(p^2\) vectorized entries is singular in the symmetric-matrix setting. Singularity reflects duplicated entries, not failure of the central limit theorem.

Key takeaway: Under fixed-dimensional multivariate normal sampling, the sample covariance has an exact Wishart representation, and its vectorized centered error is asymptotically Gaussian. The covariance between entries \(S_{ij}\) and \(S_{kl}\) on the \(\sqrt n\) scale is \(\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}\).

Check Your Understanding

Use the entrywise covariance formula and the assumptions stated above to answer these questions.

For \(p=1\), what is the limiting variance of \(\sqrt n(S-\sigma^2)\), and how does it compare with the exact variance of \(S\)?
For general \(i\) and \(j\), derive the limiting variance of \(\sqrt n(S_{ij}-\Sigma_{ij})\).
Why is the Gaussian limit for \(\operatorname{vec}(S-\Sigma)\) generally degenerate when it is written using all \(p^2\) entries?
If \(\Sigma\) is diagonal, what are the limiting covariances between distinct entries \(S_{ii}\), \(S_{jj}\), and \(S_{ij}\), for \(i\ne j\)?
Why does the exact Wishart representation have \(n-1\), rather than \(n\), degrees of freedom?