Proof of Independent Gaussian Random Variables are Jointly Gaussian

Throughout, expectations, variances, and covariances of the Gaussian random variables involved are defined and finite by Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector. Write $\mu_i=\mathbb{E}[X_i]$ and $\sigma_i^2=\operatorname{Var}(X_i)$ for $1\le i\le p$, with the variance of square-integrable random variables, and let

$$J=\{\,i\in\{1,\dots,p\}\ :\ \sigma_i^2>0\,\},$$

let $m$ be the number of elements of $J$, and enumerate $J=\{i_1,\dots,i_m\}$ with $i_1<\dots<i_m$; the case $m=0$ (empty $J$) is allowed. For $i\notin J$ we have $\sigma_i^2=0$: by the moment formulas (claim 2) of Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector applied to any Gaussian representation of the one-term tuple $(X_i)$, the variance is a sum of squares, hence nonnegative.

Step 1: Standardized variables. Suppose $m\ge1$ and fix $l\in\{1,\dots,m\}$. Let $\sigma_{i_l}$ be the positive square root of $\sigma_{i_l}^2$ and define

$$Z_l=\frac{X_{i_l}-\mu_{i_l}}{\sigma_{i_l}}.$$

By the standardization claim (claim 2) of Standardization and Cumulative Distribution Function of a Gaussian Random Variable, $Z_l$ is a standard normal random variable.

Step 2: $Z_1,\dots,Z_m$ are independent. Fix $l\in\{1,\dots,m\}$ and define $g_l:\mathbb{R}\to\mathbb{R}$ by $g_l(x)=(x-\mu_{i_l})/\sigma_{i_l}$, so that $Z_l=g_l(X_{i_l})$ pointwise on $\Omega$. For every real $a$, since $\sigma_{i_l}>0$,

$$\{x\in\mathbb{R}: g_l(x)>a\}=(\mu_{i_l}+\sigma_{i_l}a,\ \infty),$$

an open interval and hence a Borel set; by the generator criterion of Measurable Function and Real-Valued Measurable Function, $g_l$ is measurable from $\mathbb{R}$ with the Borel $\sigma$-algebra to itself. Hence for every Borel set $B$ the set $g_l^{-1}(B)$ is Borel, and

$$\{Z_l\in B\}=\{X_{i_l}\in g_l^{-1}(B)\}.$$

Now fix Borel sets $B_1,\dots,B_m$ and put $C_l=g_l^{-1}(B_l)$ for $1\le l\le m$. Since $X_1,\dots,X_p$ are independent, the events $\{X_1\in E_1\},\dots,\{X_p\in E_p\}$ are independent for the Borel choice $E_{i_l}=C_l$ ($1\le l\le m$) and $E_i=\mathbb{R}$ for $i\notin J$. Independence of events requires the product identity for every nonempty subfamily, in particular for every nonempty subfamily of the events $\{X_{i_l}\in C_l\}=\{Z_l\in B_l\}$ ($1\le l\le m$). Hence for every nonempty $S\subseteq\{1,\dots,m\}$,

$$P\Bigl(\bigcap_{l\in S}\{Z_l\in B_l\}\Bigr)=\prod_{l\in S}P(Z_l\in B_l),$$

and, the Borel sets $B_1,\dots,B_m$ being arbitrary, the random variables $Z_1,\dots,Z_m$ are independent.

Step 3: A Gaussian representation. Define real numbers $a_{ij}$ for $1\le i\le p$, $1\le j\le m$ by: $a_{ij}=\sigma_{i_l}$ if $i=i_l$ and $j=l$ for some $l\in\{1,\dots,m\}$, and $a_{ij}=0$ otherwise.

If $i=i_l\in J$, then pointwise on $\Omega$,

$$\mu_i+\sum_{j=1}^{m}a_{ij}Z_j=\mu_{i_l}+\sigma_{i_l}Z_l=X_{i_l},$$

directly from the definition of $Z_l$ in Step 1; in particular the equality holds with probability one. If $i\notin J$, then $\sigma_i^2=0$, so $P(X_i=\mu_i)=1$ by the degenerate case (claim 1) of Standardization and Cumulative Distribution Function of a Gaussian Random Variable, while the $i$-th row of $(a_{ij})$ is zero, so that

$$\mu_i+\sum_{j=1}^{m}a_{ij}Z_j=\mu_i .$$

Hence $P\bigl(X_i=\mu_i+\sum_{j=1}^{m}a_{ij}Z_j\bigr)=1$ for every $i\in\{1,\dots,p\}$, and by Steps 1 and 2 the random variables $Z_1,\dots,Z_m$ are independent standard normal. Therefore $\bigl(m,(\mu_i),(a_{ij}),(Z_j)\bigr)$ is a Gaussian representation of $(X_1,\dots,X_p)$; the case $m=0$, in which every $X_i$ is almost surely constant, is expressly allowed by Gaussian Random Vectors and Jointly Gaussian Random Variables. Thus $(X_1,\dots,X_p)$ is a Gaussian random vector.

Step 4: Distinct components are uncorrelated. By the moment formulas (claim 2) of Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector applied to the representation of Step 3, for $1\le i<k\le p$,

$$\operatorname{Cov}(X_i,X_k)=\sum_{j=1}^{m}a_{ij}\,a_{kj}.$$

Fix $j=l\in\{1,\dots,m\}$. By construction $a_{il}\ne0$ only if $i=i_l$, and $a_{kl}\ne0$ only if $k=i_l$; since $i\ne k$, at least one of the two factors vanishes, so every summand is $0$. Hence $\operatorname{Cov}(X_i,X_k)=0$. $\blacksquare$