Proof of Affine Transformations of Gaussian Random Vectors are Gaussian

Step 1 ((iii)). The constant random vector $m$ equals $m$ everywhere, so it is Gaussian by the first alternative of Gaussian Random Vectors and Jointly Gaussian Random Variables (i). A standard normal vector $Z$ in $\mathbb{R}^{d}$ satisfies $Z=0+I_{d}Z$, since $I_{d}x=x$ for every $x\in\mathbb{R}^{d}$ (Algebraic Rules for Matrix Products, Matrix-Vector Products and Transposes (v)); so $(d,Z,0,I_{d})$ is a Gaussian representation of $Z$, and $Z$ is Gaussian.

Step 2 ((i)). $BX+b$ is a random vector. The constant maps $\omega\mapsto B$ and $\omega\mapsto b$ are measurable (Measurable Maps: Generators, Composition, Continuous Maps, Restrictions, and Criteria for Real- and [0,∞]-Valued Maps (iii)). So $\omega\mapsto BX(\omega)$ is measurable as a pointwise matrix-vector product of measurable maps (Sums, Products, Transposes, Dot Products and Norms of Measurable Vector- and Matrix-Valued Maps (ii)), and $\omega\mapsto BX(\omega)+b$ as a pointwise sum of measurable maps (Sums, Products, Transposes, Dot Products and Norms of Measurable Vector- and Matrix-Valued Maps (i)).

The constant case. Let $X=m$ almost surely. By Random Vectors: Coordinates, Functions, Pairs and Almost-Sure Equality (iv), applied to $BX+b$ and the constant random vector $Bm+b$, the set $\{BX+b=Bm+b\}$ is an event. It contains the event $\{X=m\}$, which has probability one, so $BX+b=Bm+b$ almost surely, and $BX+b$ is Gaussian.

The represented case. Let $(k,Z,m,A)$ be a Gaussian representation of $X$. On the event $\{X=m+AZ\}$, by the rules for matrix-vector products in Algebraic Rules for Matrix Products, Matrix-Vector Products and Transposes (v),

$$BX+b=B(m+AZ)+b=(Bm+b)+(BA)Z .$$

The random vector $(Bm+b)+(BA)Z$ is measurable as in the first paragraph, so by Random Vectors: Coordinates, Functions, Pairs and Almost-Sure Equality (iv) the set where it equals $BX+b$ is an event; it contains $\{X=m+AZ\}$, which has probability one. Hence $BX+b=(Bm+b)+(BA)Z$ almost surely, so $(k,Z,Bm+b,BA)$ is a Gaussian representation of $BX+b$, and $BX+b$ is Gaussian.

Step 3 ((ii)). Let $S\in\mathbb{R}^{p\times d}$ be the matrix whose $j$-th row is the $i_{j}$-th row of $I_{d}$, that is, $S_{jl}=(I_{d})_{i_{j}l}$. Then for every $x\in\mathbb{R}^{d}$ and $j\in[p]$,

$$(Sx)_{j}=(I_{d}x)_{i_{j}}=x_{i_{j}},$$

using $I_{d}x=x$ (Algebraic Rules for Matrix Products, Matrix-Vector Products and Transposes (v)). Hence $(X_{i_{1}},\dots,X_{i_{p}})=SX$, which is Gaussian by Step 2 with $b=0$. For $p=1$ and $i_{1}=i$ this is the random vector $(X_{i})$ in $\mathbb{R}^{1}$, whose single coordinate $X_{i}$ is a random variable (Random Vectors: Coordinates, Functions, Pairs and Almost-Sure Equality (i)); so $X_{i}$ is a Gaussian random variable (Gaussian Random Vectors and Jointly Gaussian Random Variables (iii)).