[proofplan] We prove nonnegative definiteness by testing the spectral matrix against an arbitrary vector $(a,b) \in \mathbb{C}^2$. The linear combination $Z_t = aX_t+bY_t$ is a scalar weakly stationary process, and its spectral density is exactly the quadratic form associated with $F(\omega)$. We then prove directly that scalar spectral densities coming from absolutely summable covariance functions are nonnegative, using the nonnegativity of variances of finite Fourier sums. Hermitian symmetry follows from the covariance identity $\gamma_{VU}(h)=\overline{\gamma_{UV}(-h)}$, and the determinant inequality is the $2 \times 2$ consequence of nonnegative definiteness. [/proofplan] [step:Verify the Hermitian symmetry of the spectral density matrix] For $U,V \in \{X,Y\}$ and $h \in \mathbb{Z}$, the definition of cross-covariance gives \begin{align*} \gamma_{VU}(h) &= \mathbb{E}\big[(V_{t+h}-\mu_V)\overline{(U_t-\mu_U)}\big] \\ &= \overline{\mathbb{E}\big[(U_t-\mu_U)\overline{(V_{t+h}-\mu_V)}\big]} \\ &= \overline{\gamma_{UV}(-h)}. \end{align*} Therefore, for every $\omega \in \mathbb{R}$, \begin{align*} f_{VU}(\omega) &= \frac{1}{2\pi}\sum_{h \in \mathbb{Z}}\gamma_{VU}(h)e^{-ih\omega} \\ &= \frac{1}{2\pi}\sum_{h \in \mathbb{Z}}\overline{\gamma_{UV}(-h)}e^{-ih\omega} \\ &= \overline{\frac{1}{2\pi}\sum_{k \in \mathbb{Z}}\gamma_{UV}(k)e^{-ik\omega}} \\ &= \overline{f_{UV}(\omega)}. \end{align*} In particular $f_X(\omega)$ and $f_Y(\omega)$ are real-valued, and $f_{YX}(\omega)=\overline{f_{XY}(\omega)}$. Hence $F(\omega)$ is Hermitian. [/step] [step:Express every quadratic form of $F(\omega)$ as a scalar spectral density] Fix $a,b \in \mathbb{C}$, and define the scalar process $Z: \mathbb{Z} \to L^2(\Omega,\mathcal{F},\mathbb{P};\mathbb{C})$ by \begin{align*} Z_t := aX_t+bY_t. \end{align*} Its mean is $\mu_Z := a\mu_X+b\mu_Y$. Since $(X_t,Y_t)$ is jointly weakly stationary, $Z$ is weakly stationary, and its covariance function $\gamma_Z:\mathbb{Z}\to\mathbb{C}$ is \begin{align*} \gamma_Z(h) &:= \mathbb{E}\big[(Z_{t+h}-\mu_Z)\overline{(Z_t-\mu_Z)}\big] \\ &= |a|^2\gamma_{XX}(h)+a\overline{b}\,\gamma_{XY}(h) +b\overline{a}\,\gamma_{YX}(h)+|b|^2\gamma_{YY}(h). \end{align*} The absolute summability of the four covariance functions implies \begin{align*} \sum_{h \in \mathbb{Z}}|\gamma_Z(h)| < \infty. \end{align*} Define the scalar spectral density $f_Z:\mathbb{R}\to\mathbb{R}$ by \begin{align*} f_Z(\omega) := \frac{1}{2\pi}\sum_{h \in \mathbb{Z}}\gamma_Z(h)e^{-ih\omega}. \end{align*} Substituting the expression for $\gamma_Z$ and using absolute convergence to distribute the sum gives \begin{align*} f_Z(\omega) &= |a|^2 f_X(\omega)+a\overline{b}\,f_{XY}(\omega) +b\overline{a}\,f_{YX}(\omega)+|b|^2 f_Y(\omega) \\ &= \begin{pmatrix} a & b \end{pmatrix} F(\omega) \begin{pmatrix} \overline{a}\\ \overline{b} \end{pmatrix}. \end{align*} [/step] [step:Prove that the scalar spectral density is nonnegative] Let $Z$ be the scalar weakly stationary process from the previous step. For $n \in \mathbb{N}$ and $\omega \in \mathbb{R}$, define the finite Fourier sum $S_n(\omega):\Omega\to\mathbb{C}$ by \begin{align*} S_n(\omega) := \sum_{t=1}^{n}(Z_t-\mu_Z)e^{-it\omega}. \end{align*} Since $S_n(\omega)$ is square-integrable, $\mathbb{E}[|S_n(\omega)|^2]\geq 0$. Expanding the square and using weak stationarity, \begin{align*} \mathbb{E}[|S_n(\omega)|^2] &= \sum_{s=1}^{n}\sum_{t=1}^{n} \mathbb{E}\big[(Z_t-\mu_Z)\overline{(Z_s-\mu_Z)}\big]e^{-it\omega}e^{is\omega} \\ &= \sum_{s=1}^{n}\sum_{t=1}^{n} \gamma_Z(t-s)e^{-i(t-s)\omega} \\ &= \sum_{|h|<n}(n-|h|)\gamma_Z(h)e^{-ih\omega}. \end{align*} Dividing by $2\pi n$ gives \begin{align*} 0 \leq \frac{1}{2\pi n}\mathbb{E}[|S_n(\omega)|^2] = \frac{1}{2\pi}\sum_{|h|<n}\left(1-\frac{|h|}{n}\right)\gamma_Z(h)e^{-ih\omega}. \end{align*} Because $\sum_{h \in \mathbb{Z}}|\gamma_Z(h)|<\infty$ and \begin{align*} \left|\left(1-\frac{|h|}{n}\right)\mathbb{1}_{\{|h|<n\}}\gamma_Z(h)e^{-ih\omega}\right| \leq |\gamma_Z(h)|, \end{align*} the dominated convergence theorem for counting measure on $\mathbb{Z}$ gives \begin{align*} \lim_{n\to\infty} \frac{1}{2\pi}\sum_{|h|<n}\left(1-\frac{|h|}{n}\right)\gamma_Z(h)e^{-ih\omega} = \frac{1}{2\pi}\sum_{h \in \mathbb{Z}}\gamma_Z(h)e^{-ih\omega} = f_Z(\omega). \end{align*} The limit of nonnegative real numbers is nonnegative, so $f_Z(\omega)\geq 0$ for every $\omega \in \mathbb{R}$. [guided] The goal is to prove nonnegativity of $f_Z(\omega)$ without citing a separate scalar spectral theorem. The covariance function is positive semidefinite because it comes from variances. To expose this positivity at the frequency $\omega$, we test the process against the finite exponential vector. For $n \in \mathbb{N}$ and $\omega \in \mathbb{R}$, define \begin{align*} S_n(\omega):\Omega\to\mathbb{C}, \qquad S_n(\omega) := \sum_{t=1}^{n}(Z_t-\mu_Z)e^{-it\omega}. \end{align*} This random variable is square-integrable because it is a finite linear combination of square-integrable random variables. Hence its second moment is nonnegative: \begin{align*} \mathbb{E}[|S_n(\omega)|^2]\geq 0. \end{align*} Expanding the product gives \begin{align*} \mathbb{E}[|S_n(\omega)|^2] &= \mathbb{E}\left[ \left(\sum_{t=1}^{n}(Z_t-\mu_Z)e^{-it\omega}\right) \overline{\left(\sum_{s=1}^{n}(Z_s-\mu_Z)e^{-is\omega}\right)} \right] \\ &= \sum_{s=1}^{n}\sum_{t=1}^{n} \mathbb{E}\big[(Z_t-\mu_Z)\overline{(Z_s-\mu_Z)}\big]e^{-it\omega}e^{is\omega}. \end{align*} Weak stationarity identifies the covariance term as $\gamma_Z(t-s)$, so \begin{align*} \mathbb{E}[|S_n(\omega)|^2] &= \sum_{s=1}^{n}\sum_{t=1}^{n} \gamma_Z(t-s)e^{-i(t-s)\omega}. \end{align*} Now group the summands by the lag $h=t-s$. For a fixed $h$ with $|h|<n$, there are exactly $n-|h|$ pairs $(s,t)\in\{1,\dots,n\}^2$ such that $t-s=h$. Therefore \begin{align*} \mathbb{E}[|S_n(\omega)|^2] &= \sum_{|h|<n}(n-|h|)\gamma_Z(h)e^{-ih\omega}. \end{align*} After division by $2\pi n$, we obtain \begin{align*} 0 \leq \frac{1}{2\pi n}\mathbb{E}[|S_n(\omega)|^2] = \frac{1}{2\pi}\sum_{|h|<n}\left(1-\frac{|h|}{n}\right)\gamma_Z(h)e^{-ih\omega}. \end{align*} The right-hand side is a Fejer-type average of the Fourier series defining $f_Z$. The absolute summability assumption is exactly what allows us to pass to the limit term by term. Indeed, for each fixed $h \in \mathbb{Z}$, \begin{align*} \left(1-\frac{|h|}{n}\right)\mathbb{1}_{\{|h|<n\}}\gamma_Z(h)e^{-ih\omega} \to \gamma_Z(h)e^{-ih\omega}, \end{align*} and the pointwise absolute value is bounded by $|\gamma_Z(h)|$. Since $\sum_{h\in\mathbb{Z}}|\gamma_Z(h)|<\infty$, dominated convergence with respect to counting measure on $\mathbb{Z}$ yields \begin{align*} \lim_{n\to\infty} \frac{1}{2\pi}\sum_{|h|<n}\left(1-\frac{|h|}{n}\right)\gamma_Z(h)e^{-ih\omega} = \frac{1}{2\pi}\sum_{h \in \mathbb{Z}}\gamma_Z(h)e^{-ih\omega} = f_Z(\omega). \end{align*} Since every approximating quantity is nonnegative, its limit is nonnegative. Thus $f_Z(\omega)\geq 0$. [/guided] [/step] [step:Conclude nonnegative definiteness and derive the cross-spectral inequality] Fix $\omega \in \mathbb{R}$. For every $(a,b)\in\mathbb{C}^2$, the previous two steps give \begin{align*} \begin{pmatrix} a & b \end{pmatrix} F(\omega) \begin{pmatrix} \overline{a}\\ \overline{b} \end{pmatrix} = f_Z(\omega) \geq 0. \end{align*} Together with Hermitian symmetry, this proves that $F(\omega)$ is Hermitian nonnegative definite. It remains to extract the scalar inequality. Since $F(\omega)$ is Hermitian nonnegative definite, its diagonal entries are nonnegative: \begin{align*} f_X(\omega)\geq 0, \qquad f_Y(\omega)\geq 0. \end{align*} If $f_X(\omega)=0$, then nonnegative definiteness applied to $(1,\lambda)\in\mathbb{C}^2$ gives \begin{align*} 0 \leq \lambda f_{YX}(\omega)+\overline{\lambda}f_{XY}(\omega)+|\lambda|^2f_Y(\omega) \qquad \text{for every } \lambda \in \mathbb{C}. \end{align*} Choosing $\lambda=-r f_{XY}(\omega)$ with $r>0$ gives \begin{align*} 0 \leq -2r|f_{XY}(\omega)|^2+r^2|f_{XY}(\omega)|^2f_Y(\omega) \qquad \text{for every } r>0. \end{align*} Letting $r\downarrow 0$ forces $f_{XY}(\omega)=0$, and hence \begin{align*} |f_{XY}(\omega)|^2=0=f_X(\omega)f_Y(\omega). \end{align*} If $f_X(\omega)>0$, use nonnegative definiteness on the vector $\left(-f_{YX}(\omega)/f_X(\omega),1\right)$. This gives \begin{align*} 0 &\leq \begin{pmatrix} -\frac{f_{YX}(\omega)}{f_X(\omega)} & 1 \end{pmatrix} F(\omega) \begin{pmatrix} -\frac{f_{XY}(\omega)}{f_X(\omega)}\\ 1 \end{pmatrix} \\ &= \frac{|f_{XY}(\omega)|^2}{f_X(\omega)} -\frac{|f_{XY}(\omega)|^2}{f_X(\omega)} -\frac{|f_{XY}(\omega)|^2}{f_X(\omega)} +f_Y(\omega) \\ &= f_Y(\omega)-\frac{|f_{XY}(\omega)|^2}{f_X(\omega)}. \end{align*} Multiplying by the positive number $f_X(\omega)$ and rearranging yields \begin{align*} |f_{XY}(\omega)|^2\leq f_X(\omega)f_Y(\omega). \end{align*} Therefore the desired inequality holds for every $\omega \in \mathbb{R}$. [/step]