The Multivariate Central Limit Theorem For Sample Moments

Why Sample Moments Have Joint Limits

Econometric calculations often use several summaries of the same sample: a sample mean, a second raw moment, or a collection of means and proportions. These summaries are usually dependent because they are calculated from the same observations. A central limit theorem for each summary separately does not, by itself, describe their joint uncertainty. The multivariate central limit theorem does: it gives a joint normal approximation and a covariance matrix that records both variances and cross-covariances.

The key move is to treat a vector of functions of each observation as one random vector. Averaging those vectors produces the sample moments of interest. The iid assumption and a finite covariance matrix then give a joint limit at the familiar square-root-of-sample-size rate.

Definition: Let $X_1,\ldots,X_n$ be iid observations, and let $g(X_i)$ be a $k$-dimensional column vector with finite mean $\eta=E[g(X_i)]$ and finite covariance matrix $\Omega=\operatorname{Var}(g(X_i))$. Define $\bar G_n=n^{-1}\sum_{i=1}^n g(X_i)$. The multivariate central limit theorem states that $$ \sqrt{n}(\bar G_n-\eta)\xrightarrow{d}N_k(0,\Omega). $$ Here $N_k(0,\Omega)$ denotes a $k$-dimensional normal distribution with mean vector zero and covariance matrix $\Omega$. The matrix $\Omega$ may be singular.

The condition of finite covariance means that every component of $g(X_i)$ has a finite second moment. In particular, every entry of $\Omega$ exists. No independence among the components of $g(X_i)$ is required: the components are calculated from the same observation, and their dependence is represented by the off-diagonal entries of $\Omega$. Independence is required across the observations $i$.

From Vectors of Observations to Sample Moments

For a scalar observation $X_i$, define its population raw moment of order $r$ as $\mu_r=E[X_i^r]$, and its sample raw moment as $\hat{\mu}_r=n^{-1}\sum_{i=1}^n X_i^r$. To study moments of orders $1,\ldots,q$ jointly, use the observation-level vector $g(X_i)=(X_i,X_i^2,\ldots,X_i^q)'$. Its sample average is precisely $(\hat{\mu}_1,\hat{\mu}_2,\ldots,\hat{\mu}_q)'$.

The covariance between the $r$th and $s$th components of this vector is $\operatorname{Cov}(X_i^r,X_i^s)$. When the required moments exist, this is $\mu_{r+s}-\mu_r\mu_s$. Thus the limiting covariance matrix for the vector of sample raw moments has entries $\Omega_{rs}=\mu_{r+s}-\mu_r\mu_s$. This calculation identifies which population moments are needed: for moments through order $q$, a sufficient condition is $E[|X_i|^{2q}]<\infty$.

\sqrt{n} \begin{pmatrix} \hat{\mu}_1-\mu_1\\ \hat{\mu}_2-\mu_2\\ \vdots\\ \hat{\mu}_q-\mu_q \end{pmatrix} \xrightarrow{d} N_q(0,\Omega), \qquad \Omega_{rs}=\mu_{r+s}-\mu_r\mu_s.

The covariance matrix describes the joint limit, not just a list of separate variances. In particular, a positive or negative off-diagonal entry indicates that the corresponding centered sample moments tend to fluctuate together or in opposite directions. It is generally incorrect to replace $\Omega$ by a diagonal matrix unless the relevant components really are uncorrelated.

Why the Vector Theorem Holds

A short justification shows how the scalar iid central limit theorem extends to vectors. For any fixed $k$-dimensional vector $a$, take the scalar projection $a'\bar G_n$. It is the sample average of the iid scalar variables $a'g(X_i)$, whose mean is $a'\eta$ and whose variance is $a'\Omega a$. The scalar central limit theorem therefore gives a normal limit for that projection.

\sqrt{n}\,a'(\bar G_n-\eta) = \frac{1}{\sqrt{n}}\sum_{i=1}^n a'\{g(X_i)-\eta\} \xrightarrow{d}N(0,a'\Omega a).

The Cramér–Wold characterization says that a sequence of random vectors converges in distribution to a specified vector distribution if every fixed scalar projection converges to the corresponding projection of that distribution. The projections above are exactly those of $N_k(0,\Omega)$, because a projection $a'Z$ of $Z\sim N_k(0,\Omega)$ is normal with variance $a'\Omega a$. Applying this characterization establishes the stated vector limit. If $a'\Omega a=0$ for some $a$, that projection has zero variance and is identically zero after centering; the argument still applies.

Worked Example: Joint Limits for the First Two Raw Moments

Suppose, hypothetically, that each iid observation $X_i$ takes values $0,1,2$ with probabilities $0.2,0.5,0.3$, respectively. Find the limiting covariance matrix for $(\hat{\mu}_1,\hat{\mu}_2)'$, where $\hat{\mu}_1$ is the sample mean and $\hat{\mu}_2$ is the sample second raw moment.

First calculate the population moments needed for the covariance entries:

\mu_1=0(0.2)+1(0.5)+2(0.3)=1.1,\qquad \mu_2=0^2(0.2)+1^2(0.5)+2^2(0.3)=1.7, $$ $$ \mu_3=0^3(0.2)+1^3(0.5)+2^3(0.3)=2.9,\qquad \mu_4=0^4(0.2)+1^4(0.5)+2^4(0.3)=5.3.

The variance of $X_i$ is $\mu_2-\mu_1^2=1.7-1.21=0.49$. The covariance of $X_i$ and $X_i^2$ is $\mu_3-\mu_1\mu_2=2.9-(1.1)(1.7)=1.03$. The variance of $X_i^2$ is $\mu_4-\mu_2^2=5.3-2.89=2.41$. Therefore,

\sqrt{n} \begin{pmatrix} \hat{\mu}_1-1.1\\ \hat{\mu}_2-1.7 \end{pmatrix} \xrightarrow{d} N_2\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix} 0.49 & 1.03\\ 1.03 & 2.41 \end{pmatrix} \right).

All moments exist because $X_i$ has finite support, so the theorem's moment condition holds; the observations are assumed iid as specified. For $n=100$, the joint normal approximation assigns covariance matrix $\Omega/n$ to the unscaled sample-moment deviations. The approximate standard errors are $\sqrt{0.49/100}=0.0700$ for $\hat{\mu}_1$ and $\sqrt{2.41/100}=0.155$ for $\hat{\mu}_2$. These are large-sample approximations, not exact finite-sample statements. The positive covariance indicates that samples with a relatively high first raw moment tend also to have a relatively high second raw moment.

Joint Limits for Means and Proportions

The same construction applies when the components are not powers. For example, a vector may include a numeric outcome and an indicator for whether that outcome exceeds a threshold. The sample average of the indicator is a sample proportion; the covariance matrix then accounts for its dependence with the sample mean. This is useful whenever multiple summaries come from the same observations.

Worked Example: Two Sample Proportions From Paired Outcomes

Suppose, hypothetically, each iid unit has two binary outcomes $X_i,Y_i$, with joint probabilities $P(0,0)=0.4$, $P(1,0)=0.2$, $P(0,1)=0.1$, and $P(1,1)=0.3$. Find the joint limiting distribution of the two sample proportions $(\bar X,\bar Y)'$, where $\bar X=n^{-1}\sum X_i$ and $\bar Y=n^{-1}\sum Y_i$.

The means are $E[X_i]=P(X_i=1)=0.2+0.3=0.5$ and $E[Y_i]=P(Y_i=1)=0.1+0.3=0.4$. Since a binary variable equals its square, $\operatorname{Var}(X_i)=0.5-(0.5)^2=0.25$ and $\operatorname{Var}(Y_i)=0.4-(0.4)^2=0.24$. Also, $E[X_iY_i]=P(X_i=1,Y_i=1)=0.3$, so their covariance is $0.3-(0.5)(0.4)=0.10$.

\sqrt{n} \begin{pmatrix} \bar X-0.5\\ \bar Y-0.4 \end{pmatrix} \xrightarrow{d} N_2\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix} 0.25 & 0.10\\ 0.10 & 0.24 \end{pmatrix} \right).

The observations are iid pairs, and the components are bounded, so the finite-covariance condition is satisfied. For $n=400$, the approximate standard errors of the two proportions are $\sqrt{0.25/400}=0.0250$ and $\sqrt{0.24/400}=0.0245$; their approximate covariance is $0.10/400=0.000250$. Treating these proportions as independent would discard that covariance and misdescribe their joint sampling variation.

More Than Two Moments and Singular Limits

A valid multivariate normal limit need not have a positive-definite covariance matrix. If one component of the observation vector is a linear combination of others, then the same relation holds for the sample averages, and the limiting covariance matrix is singular. This does not invalidate the multivariate central limit theorem: it means the normal limit is concentrated on a lower-dimensional subspace.

Worked Example: A Singular Limit for Three Sample Moments

Suppose, hypothetically, $X_i$ takes values $-1,0,1$ with probabilities $0.25,0.5,0.25$, independently across observations. Consider the vector of sample moments $(\hat{\mu}_1,\hat{\mu}_2,\hat{\mu}_3)'$. By symmetry, $\mu_1=\mu_3=0$, while $\mu_2=0.5$. Also, $\mu_4=0.5$, $\mu_5=0$, and $\mu_6=0.5$.

Using $\Omega_{rs}=\mu_{r+s}-\mu_r\mu_s$, the diagonal entries are $\operatorname{Var}(X_i)=0.5$, $\operatorname{Var}(X_i^2)=0.5-(0.5)^2=0.25$, and $\operatorname{Var}(X_i^3)=0.5$. The first–second and second–third covariances are zero by symmetry. The first–third covariance is $E[X_i^4]-E[X_i]E[X_i^3]=0.5$. Thus,

\sqrt{n} \begin{pmatrix} \hat{\mu}_1\\ \hat{\mu}_2-0.5\\ \hat{\mu}_3 \end{pmatrix} \xrightarrow{d} N_3\left( \begin{pmatrix}0\\0\\0\end{pmatrix}, \begin{pmatrix} 0.5 & 0 & 0.5\\ 0 & 0.25 & 0\\ 0.5 & 0 & 0.5 \end{pmatrix} \right).

Here $X_i^3=X_i$ for every possible observation, so $\hat{\mu}_3=\hat{\mu}_1$ in every sample. The first and third components of the centered vector are therefore identical, which explains the singular covariance matrix. The finite support ensures all required moments exist, and iid sampling supplies the other condition.

Common Pitfalls and Interpretation

Ignoring dependence within an observation. Components such as $X_i$ and $X_i^2$ are generally dependent. Calculate their covariance rather than assuming it is zero.
Checking only the variance of $X_i$. A joint limit for moments through order $q$ requires finite second moments of every component, including $X_i^q$. A convenient sufficient condition is $E[|X_i|^{2q}]<\infty$.
Using the wrong scaling. The theorem concerns $\sqrt n(\bar G_n-\eta)$. Consequently, the unscaled average has approximate covariance $\Omega/n$, not $\Omega$.
Rejecting a singular covariance matrix. Singularity can reflect exact relationships among components. It does not make the joint limit invalid.
Confusing an approximation with an exact distribution. The multivariate normal statement is a limit as $n$ grows. It does not claim that the finite-sample vector is normally distributed.

Key takeaway: Choose the observation-level vector whose averages give the sample moments of interest, verify that its components have finite second moments, and calculate the covariance matrix from that vector. The multivariate central limit theorem then provides a joint normal limit, including the dependence among the sample moments.

Check Your Understanding

Use the vector formulation and covariance calculations to answer the following.

For a sample raw moment vector through order $q$, what finite-moment condition is sufficient for the multivariate central limit theorem?
Why can the limiting covariance between $\hat{\mu}_1$ and $\hat{\mu}_2$ be nonzero even though both are calculated from the same iid sample?
If the limiting covariance matrix for $\sqrt{n}(\bar G_n-\eta)$ is $\Omega$, what is the corresponding large-sample covariance approximation for $\bar G_n$?
For paired binary outcomes, what population quantity determines the off-diagonal covariance between their sample proportions?
What does a singular covariance matrix imply about the joint normal limit, and does it invalidate the central limit theorem?