[proofplan] We first verify that the covariance entries are well-defined [real numbers](/page/Real%20Numbers) using the finite second moment hypothesis. Symmetry follows directly from commutativity of multiplication in $\mathbb{R}$. To prove positive semidefiniteness, we evaluate the quadratic form $a^\top \Sigma a$ and identify it with the expectation of the square of the centered scalar projection $a^\top(X-\mu)$. [/proofplan] [step:Verify that the covariance entries are finite] For each $i \in \{1,\dots,p\}$, the hypothesis $\mathbb E[X_i^2] < \infty$ implies $\mathbb E[|X_i|] < \infty$, since $|t| \leq 1 + t^2$ for every $t \in \mathbb{R}$. Hence $\mu_i = \mathbb E[X_i]$ is a finite real number. Fix $i,j \in \{1,\dots,p\}$. Define the centered coordinate random variables $Y_i: \Omega \to \mathbb{R}$ and $Y_j: \Omega \to \mathbb{R}$ by \begin{align*} Y_i(\omega) &:= X_i(\omega) - \mu_i, & Y_j(\omega) &:= X_j(\omega) - \mu_j. \end{align*} Using $(u-v)^2 \leq 2u^2 + 2v^2$ for real numbers $u,v$, we have \begin{align*} \mathbb E[Y_i^2] &= \int_\Omega (X_i(\omega)-\mu_i)^2\,d\mathbb P(\omega) \\ &\leq 2\int_\Omega X_i(\omega)^2\,d\mathbb P(\omega) + 2\mu_i^2 < \infty. \end{align*} The same argument gives $\mathbb E[Y_j^2] < \infty$. Finally, by $2|uv| \leq u^2+v^2$, \begin{align*} \mathbb E[|Y_iY_j|] &= \int_\Omega |Y_i(\omega)Y_j(\omega)|\,d\mathbb P(\omega) \\ &\leq \frac{1}{2}\int_\Omega Y_i(\omega)^2\,d\mathbb P(\omega) + \frac{1}{2}\int_\Omega Y_j(\omega)^2\,d\mathbb P(\omega) < \infty. \end{align*} Therefore $\Sigma_{ij} = \mathbb E[Y_iY_j]$ is well-defined and finite. [/step] [step:Use commutativity of multiplication to prove symmetry] Let $i,j \in \{1,\dots,p\}$. With $Y_i = X_i-\mu_i$ and $Y_j = X_j-\mu_j$ as above, commutativity of multiplication in $\mathbb{R}$ gives $Y_iY_j = Y_jY_i$ pointwise on $\Omega$. Therefore \begin{align*} \Sigma_{ij} &= \mathbb E[(X_i-\mu_i)(X_j-\mu_j)] \\ &= \mathbb E[(X_j-\mu_j)(X_i-\mu_i)] \\ &= \Sigma_{ji}. \end{align*} Thus $\Sigma^\top=\Sigma$. [/step] [step:Identify the quadratic form with the variance of a scalar projection] Fix $a=(a_1,\dots,a_p) \in \mathbb{R}^p$. Define the centered scalar projection $Z_a: \Omega \to \mathbb{R}$ by \begin{align*} Z_a(\omega) := \sum_{i=1}^p a_i(X_i(\omega)-\mu_i). \end{align*} Since $Z_a$ is a finite linear combination of measurable real-valued random variables, it is measurable. Moreover, using the finite second moments of the centered coordinates, $Z_a^2$ is integrable. Indeed, \begin{align*} |Z_a(\omega)|^2 &= \left|\sum_{i=1}^p a_i(X_i(\omega)-\mu_i)\right|^2 \\ &\leq p\sum_{i=1}^p a_i^2(X_i(\omega)-\mu_i)^2, \end{align*} so \begin{align*} \mathbb E[Z_a^2] \leq p\sum_{i=1}^p a_i^2\mathbb E[(X_i-\mu_i)^2] < \infty. \end{align*} By the definition of matrix multiplication and finite additivity of expectation over finite sums, \begin{align*} a^\top \Sigma a &= \sum_{i=1}^p\sum_{j=1}^p a_i\Sigma_{ij}a_j \\ &= \sum_{i=1}^p\sum_{j=1}^p a_i a_j\mathbb E[(X_i-\mu_i)(X_j-\mu_j)] \\ &= \mathbb E\left[\sum_{i=1}^p\sum_{j=1}^p a_i a_j(X_i-\mu_i)(X_j-\mu_j)\right] \\ &= \mathbb E\left[\left(\sum_{i=1}^p a_i(X_i-\mu_i)\right)^2\right] \\ &= \mathbb E[Z_a^2]. \end{align*} Since $Z_a^2 \geq 0$ pointwise on $\Omega$, its expectation is nonnegative: \begin{align*} a^\top \Sigma a = \mathbb E[Z_a^2] \geq 0. \end{align*} Because $a \in \mathbb{R}^p$ was arbitrary, $\Sigma$ is positive semidefinite. [/step]