[proofplan]
Each part is an "if and only if" statement, requiring both directions. For the forward direction, we derive the properties (symmetry, positive semidefiniteness, diagonal entries) directly from the definition of covariance and the non-negativity of variance. For the reverse direction, we explicitly construct a random vector with the given matrix as its covariance (or correlation) by using a multivariate Gaussian with the specified matrix as its parameter. The Gaussian construction is valid precisely because the given matrix has the required algebraic properties.
[/proofplan]
[step:Prove that every covariance matrix is symmetric and positive semidefinite]
Let $X \in \mathbb{R}^p$ be a random vector with $\mathbb{E}[X] = 0$ and $\operatorname{Cov}(X) = \Sigma$. Then $\Sigma_{ij} = \mathbb{E}[X_i X_j]$ for each $i, j \in \{1, \ldots, p\}$.
**Symmetry:** $\Sigma_{ij} = \mathbb{E}[X_i X_j] = \mathbb{E}[X_j X_i] = \Sigma_{ji}$, so $\Sigma = \Sigma^\top$.
**Positive semidefiniteness:** For any $a \in \mathbb{R}^p$,
\begin{align*}
a^\top \Sigma a = a^\top \mathbb{E}[X X^\top] a = \mathbb{E}[a^\top X X^\top a] = \mathbb{E}[(a^\top X)^2] \geq 0,
\end{align*}
where we used linearity of expectation to move $a$ inside $\mathbb{E}$ (since $a$ is deterministic), and the final inequality holds because $(a^\top X)^2 \geq 0$ pointwise.
[guided]
The forward direction extracts the two properties from the definition $\Sigma = \mathbb{E}[XX^\top]$ (using $\mathbb{E}[X] = 0$).
**Symmetry** is immediate: $(\mathbb{E}[XX^\top])_{ij} = \mathbb{E}[X_i X_j] = \mathbb{E}[X_j X_i] = (\mathbb{E}[XX^\top])_{ji}$.
**Positive semidefiniteness** follows from the non-negativity of variance. For any $a \in \mathbb{R}^p$, the scalar $a^\top X$ has variance $\operatorname{Var}(a^\top X) = a^\top \operatorname{Cov}(X) a = a^\top \Sigma a$. Since variance is non-negative, $a^\top \Sigma a \geq 0$. More directly, $a^\top \Sigma a = \mathbb{E}[(a^\top X)^2] \geq 0$.
[/guided]
[/step]
[step:Construct a random vector with any given symmetric positive semidefinite covariance]
Conversely, let $\Sigma \in \mathbb{R}^{p \times p}$ be symmetric and positive semidefinite. We must exhibit a random vector $X \in \mathbb{R}^p$ with $\mathbb{E}[X] = 0$ and $\operatorname{Cov}(X) = \Sigma$.
Since $\Sigma$ is symmetric positive semidefinite, it admits a spectral decomposition $\Sigma = V \Lambda V^\top$, where $V \in \mathbb{R}^{p \times p}$ is orthogonal and $\Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_p)$ with $\lambda_i \geq 0$ for all $i$. Define $\Sigma^{1/2} = V \Lambda^{1/2} V^\top$, where $\Lambda^{1/2} = \operatorname{diag}(\sqrt{\lambda_1}, \ldots, \sqrt{\lambda_p})$. Then $\Sigma^{1/2}$ is symmetric and satisfies $\Sigma^{1/2} \Sigma^{1/2} = \Sigma$.
Let $Z \sim \mathcal{N}(0, I_p)$ be a standard Gaussian vector in $\mathbb{R}^p$, and define $X = \Sigma^{1/2} Z$. Then $\mathbb{E}[X] = \Sigma^{1/2} \mathbb{E}[Z] = 0$ and
\begin{align*}
\operatorname{Cov}(X) = \mathbb{E}[XX^\top] = \Sigma^{1/2} \mathbb{E}[ZZ^\top] (\Sigma^{1/2})^\top = \Sigma^{1/2} I_p \Sigma^{1/2} = \Sigma.
\end{align*}
[guided]
The goal is to construct a concrete random vector with the desired covariance. The natural candidate is a linear transformation of a standard Gaussian. Why does this work?
A standard Gaussian vector $Z \sim \mathcal{N}(0, I_p)$ has identity covariance. If we set $X = AZ$ for any $p \times p$ matrix $A$, then $\operatorname{Cov}(X) = A I_p A^\top = AA^\top$. We need $AA^\top = \Sigma$, so we need a "square root" of $\Sigma$.
The spectral theorem guarantees that any symmetric positive semidefinite matrix $\Sigma$ can be written as $\Sigma = V\Lambda V^\top$ with $V$ orthogonal and $\Lambda$ diagonal with non-negative entries. The matrix square root $\Sigma^{1/2} = V\Lambda^{1/2}V^\top$ satisfies $\Sigma^{1/2} \Sigma^{1/2} = V\Lambda^{1/2}V^\top V\Lambda^{1/2}V^\top = V\Lambda V^\top = \Sigma$ (using $V^\top V = I_p$). Setting $A = \Sigma^{1/2}$, we get $X = \Sigma^{1/2} Z$ with $\operatorname{Cov}(X) = \Sigma$.
Note that $\Sigma$ need not be strictly positive definite; if some eigenvalues are zero, the corresponding components of $X$ are deterministically zero, which is permitted.
[/guided]
[/step]
[step:Prove the forward direction for correlation matrices]
Let $X \in \mathbb{R}^p$ with $\mathbb{E}[X] = 0$ and $\operatorname{Var}(X_i) > 0$ for each $i$ (otherwise the correlation matrix is not defined). The correlation matrix is
\begin{align*}
R_{ij} = \frac{\operatorname{Cov}(X_i, X_j)}{\sqrt{\operatorname{Var}(X_i)\operatorname{Var}(X_j)}} = \frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}.
\end{align*}
Equivalently, $R = D^{-1} \Sigma D^{-1}$ where $D = \operatorname{diag}(\sqrt{\Sigma_{11}}, \ldots, \sqrt{\Sigma_{pp}})$.
**Symmetry:** $R_{ij} = R_{ji}$ since $\Sigma_{ij} = \Sigma_{ji}$.
**Ones on the diagonal:** $R_{ii} = \Sigma_{ii}/\Sigma_{ii} = 1$.
**Positive semidefiniteness:** For any $a \in \mathbb{R}^p$, let $b = D^{-1}a$. Then $a^\top R a = a^\top D^{-1}\Sigma D^{-1} a = b^\top \Sigma b \geq 0$, since $\Sigma$ is positive semidefinite.
[/step]
[step:Construct a random vector with any given symmetric positive semidefinite matrix with unit diagonal as its correlation]
Conversely, let $R \in \mathbb{R}^{p \times p}$ be symmetric, positive semidefinite, and have $R_{ii} = 1$ for all $i$.
Since $R$ is symmetric positive semidefinite with $R_{ii} = 1$, it is in particular a valid covariance matrix (by part (i)). Construct $X \sim \mathcal{N}(0, R)$ using the square-root construction from the second step above. Then $\operatorname{Cov}(X) = R$, so $\operatorname{Var}(X_i) = R_{ii} = 1$ for each $i$. The correlation matrix of $X$ is therefore
\begin{align*}
\operatorname{Cor}(X_i, X_j) = \frac{\operatorname{Cov}(X_i, X_j)}{\sqrt{\operatorname{Var}(X_i)\operatorname{Var}(X_j)}} = \frac{R_{ij}}{\sqrt{1 \cdot 1}} = R_{ij}.
\end{align*}
Hence $R$ is the correlation matrix of $X$.
[guided]
The trick for the reverse direction of part (ii) is to notice that a matrix with ones on the diagonal and positive semidefiniteness is simultaneously a valid covariance matrix (by part (i)). If we use it directly as the covariance, the resulting random vector automatically has unit variances, so the correlation matrix equals the covariance matrix.
Concretely, $X \sim \mathcal{N}(0, R)$ has $\operatorname{Var}(X_i) = R_{ii} = 1$, so the standardization step in the correlation formula is trivial: $\operatorname{Cor}(X_i, X_j) = \operatorname{Cov}(X_i, X_j)/1 = R_{ij}$. This completes the characterisation: a matrix is a correlation matrix if and only if it is symmetric, positive semidefinite, and has ones on the diagonal.
[/guided]
[/step]