[proofplan] We compute the covariance of $X$ directly from the structural equation $X=\Lambda F+\varepsilon$. Since both latent factors and idiosyncratic errors have mean zero, $X$ is centred, so its covariance is its second moment matrix. Expanding $(\Lambda F+\varepsilon)(\Lambda F+\varepsilon)^\top$ produces two main covariance terms and two mixed terms; independence and zero means make the mixed terms vanish. Substituting $\operatorname{Cov}(F)=I_k$ and $\operatorname{Cov}(\varepsilon)=\Psi$ gives the desired formula. [/proofplan] [step:Show that the observed vector is centred] Since $F$ and $\varepsilon$ are Gaussian random vectors, all first and second moments appearing below are finite. By linearity of expectation, \begin{align*} \mathbb E[X] &= \mathbb E[\Lambda F+\varepsilon] \\ &= \Lambda \mathbb E[F]+\mathbb E[\varepsilon] \\ &= \Lambda 0+0 \\ &=0. \end{align*} Therefore \begin{align*} \operatorname{Cov}(X)=\mathbb E[XX^\top]. \end{align*} [/step] [step:Expand the second moment matrix of $X$] Using $X=\Lambda F+\varepsilon$, matrix multiplication, and linearity of expectation, \begin{align*} \operatorname{Cov}(X) &= \mathbb E[(\Lambda F+\varepsilon)(\Lambda F+\varepsilon)^\top] \\ &= \mathbb E[(\Lambda F+\varepsilon)(F^\top\Lambda^\top+\varepsilon^\top)] \\ &= \mathbb E[\Lambda FF^\top\Lambda^\top] +\mathbb E[\Lambda F\varepsilon^\top] +\mathbb E[\varepsilon F^\top\Lambda^\top] +\mathbb E[\varepsilon\varepsilon^\top] \\ &= \Lambda \mathbb E[FF^\top]\Lambda^\top +\Lambda \mathbb E[F\varepsilon^\top] +\mathbb E[\varepsilon F^\top]\Lambda^\top +\mathbb E[\varepsilon\varepsilon^\top]. \end{align*} [guided] Because $X$ is centred, the covariance computation reduces to expanding the second moment matrix $\mathbb E[XX^\top]$. Substituting the factor-model equation gives \begin{align*} \operatorname{Cov}(X) &= \mathbb E[(\Lambda F+\varepsilon)(\Lambda F+\varepsilon)^\top]. \end{align*} The transpose of $\Lambda F+\varepsilon$ is $F^\top\Lambda^\top+\varepsilon^\top$, so distributivity of matrix multiplication gives \begin{align*} (\Lambda F+\varepsilon)(\Lambda F+\varepsilon)^\top &= \Lambda FF^\top\Lambda^\top +\Lambda F\varepsilon^\top +\varepsilon F^\top\Lambda^\top +\varepsilon\varepsilon^\top. \end{align*} Taking expectations entrywise and using linearity of expectation yields \begin{align*} \operatorname{Cov}(X) &= \Lambda \mathbb E[FF^\top]\Lambda^\top +\Lambda \mathbb E[F\varepsilon^\top] +\mathbb E[\varepsilon F^\top]\Lambda^\top +\mathbb E[\varepsilon\varepsilon^\top]. \end{align*} This separates the covariance into the contribution from the common factors, the contribution from the idiosyncratic errors, and two mixed terms. [/guided] [/step] [step:Use independence and zero means to remove the mixed covariance terms] For $i\in\{1,\dots,k\}$ and $j\in\{1,\dots,d\}$, let $F_i:\Omega\to\mathbb R$ and $\varepsilon_j:\Omega\to\mathbb R$ denote the coordinate random variables of $F$ and $\varepsilon$. Since $F$ and $\varepsilon$ are independent, $F_i$ and $\varepsilon_j$ are independent. Hence \begin{align*} \mathbb E[F_i\varepsilon_j] &= \mathbb E[F_i]\mathbb E[\varepsilon_j] \\ &= 0\cdot 0 \\ &=0. \end{align*} Thus every entry of $\mathbb E[F\varepsilon^\top]\in\mathbb R^{k\times d}$ is zero, so \begin{align*} \mathbb E[F\varepsilon^\top]=0. \end{align*} Taking transposes gives \begin{align*} \mathbb E[\varepsilon F^\top] = \mathbb E[F\varepsilon^\top]^\top = 0. \end{align*} [guided] The mixed matrices measure linear dependence between the latent factors and the error coordinates. We show entry by entry that they vanish. For $i\in\{1,\dots,k\}$ and $j\in\{1,\dots,d\}$, define $F_i:\Omega\to\mathbb R$ to be the $i$-th coordinate of $F$ and $\varepsilon_j:\Omega\to\mathbb R$ to be the $j$-th coordinate of $\varepsilon$. Independence of the random vectors $F$ and $\varepsilon$ implies independence of each pair $F_i$ and $\varepsilon_j$. The $(i,j)$ entry of $\mathbb E[F\varepsilon^\top]$ is $\mathbb E[F_i\varepsilon_j]$. Since $F_i$ and $\varepsilon_j$ are independent and integrable, \begin{align*} \mathbb E[F_i\varepsilon_j] &= \mathbb E[F_i]\mathbb E[\varepsilon_j]. \end{align*} The assumptions $\mathbb E[F]=0$ and $\mathbb E[\varepsilon]=0$ mean $\mathbb E[F_i]=0$ and $\mathbb E[\varepsilon_j]=0$, so \begin{align*} \mathbb E[F_i\varepsilon_j] &=0. \end{align*} Since this holds for every entry, \begin{align*} \mathbb E[F\varepsilon^\top]=0. \end{align*} The other mixed matrix is the transpose: \begin{align*} \mathbb E[\varepsilon F^\top] = \mathbb E[F\varepsilon^\top]^\top = 0. \end{align*} Thus independence and centring eliminate both cross terms. [/guided] [/step] [step:Substitute the factor and error covariance matrices] Because $\mathbb E[F]=0$ and $\mathbb E[\varepsilon]=0$, \begin{align*} \mathbb E[FF^\top] &=\operatorname{Cov}(F) =I_k, \\ \mathbb E[\varepsilon\varepsilon^\top] &=\operatorname{Cov}(\varepsilon) =\Psi. \end{align*} Substituting these identities and the vanishing mixed terms into the expansion gives \begin{align*} \operatorname{Cov}(X) &= \Lambda I_k\Lambda^\top+\Lambda 0+0\Lambda^\top+\Psi \\ &= \Lambda\Lambda^\top+\Psi. \end{align*} This proves the covariance formula for the Gaussian factor model. [/step]