Proof of Affine Transformations of Gaussian Random Vectors are Gaussian

Fix a Gaussian representation $\bigl(m,(\mu_k),(a_{kj}),(Z_j)\bigr)$ of $(X_1,\dots,X_d)$, so that the events $E_k=\bigl\{X_k=\mu_k+\sum_{j=1}^{m}a_{kj}Z_j\bigr\}$ satisfy $P(E_k)=1$ for $1\le k\le d$; when $m=0$ all sums over $j$ below are empty and equal to $0$, in accordance with Gaussian Random Vectors and Jointly Gaussian Random Variables. Each $Y_i$ is a random variable by the closure of random variables under sums and scalar multiples recorded in the preliminaries of Square-Integrable Random Variables and the Mean-Square Inner Product.

Step 1: A representation for $(Y_1,\dots,Y_p)$. Let $\Omega_0=\bigcap_{k=1}^{d}E_k$. Its complement is the finite union of the complements $\Omega\setminus E_k$, each of probability $0$, so by the monotonicity and countable additivity of the probability measure $P$ (finite subadditivity being the special case with cofinitely many empty sets), $P(\Omega\setminus\Omega_0)=0$, i.e., $P(\Omega_0)=1$. For every $\omega\in\Omega_0$ and $1\le i\le p$, substituting and rearranging the finite sums,

where

Hence $P\bigl(Y_i=\nu_i+\sum_{j}b_{ij}Z_j\bigr)\ge P(\Omega_0)=1$ for each $i$, so $\bigl(m,(\nu_i),(b_{ij}),(Z_j)\bigr)$ is a Gaussian representation of $(Y_1,\dots,Y_p)$, and $(Y_1,\dots,Y_p)$ is a Gaussian random vector.

Step 2: Mean vector and covariances. Apply Claim 2 of Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector twice: to the representation of Step 1, giving $\mathbb{E}[Y_i]=\nu_i$ and $\operatorname{Cov}(Y_i,Y_l)=\sum_{j}b_{ij}b_{lj}$, and to the original representation, giving $\mathbb{E}[X_k]=\mu_k$ and $\operatorname{Cov}(X_k,X_{k'})=\sum_{j}a_{kj}a_{k'j}$. Then

and, exchanging the order of the finite sums,

Step 3: Particular cases. The linear-combination claim is the case $p=1$ with $c_1=c$ and $M_{1k}=c_k$. The subfamily claim is the case $c_i=0$ and $M_{ik}=1$ if $k=i_i$ and $M_{ik}=0$ otherwise, where $i_1<\dots<i_p$ are the selected indices, so that $Y_i=X_{i_i}$ everywhere on $\Omega$. Sums and differences are the case $p=1$ with coefficients in $\{1,-1\}$ and $c=0$. $\blacksquare$