[proofplan] We test the autocovariance kernel against an arbitrary finite vector of real coefficients. The key construction is the centered linear combination $Y=\sum_{i=1}^n a_i(X_i-\mu)$, whose square has non-negative expectation. Expanding $\mathbb{E}[Y^2]$ and using weak stationarity identifies that expectation with the required quadratic form. [/proofplan]

[step:Form the centered linear combination associated to the coefficient vector] Fix $n\in\mathbb{N}$ and fix real numbers $a_1,\dots,a_n\in\mathbb{R}$. Let $(\Omega,\mathcal{F},\mathbb{P})$ denote the probability space on which the process $(X_t)_{t\in\mathbb{Z}}$ is defined, let $\mathcal{B}(\mathbb{R})$ denote the Borel $\sigma$-algebra on $\mathbb{R}$, and let $L^2(\Omega,\mathcal{F},\mathbb{P})$ denote the space of real-valued square-integrable random variables on this probability space, modulo equality $\mathbb{P}$-almost surely. Let $\mu\in\mathbb{R}$ denote the common mean $\mu=\mathbb{E}[X_t]$, which is independent of $t$ by weak stationarity. For each $i\in\{1,\dots,n\}$, define the centered random variable \begin{align*} Z_i:(\Omega,\mathcal{F})&\to(\mathbb{R},\mathcal{B}(\mathbb{R}))\\ \omega&\mapsto X_i(\omega)-\mu. \end{align*} Since $X_i$ has finite second moment, $Z_i\in L^2(\Omega,\mathcal{F},\mathbb{P})$. Define the finite linear combination \begin{align*} Y:(\Omega,\mathcal{F})&\to(\mathbb{R},\mathcal{B}(\mathbb{R}))\\ \omega&\mapsto \sum_{i=1}^{n}a_i Z_i(\omega). \end{align*} Because $L^2(\Omega,\mathcal{F},\mathbb{P})$ is closed under finite linear combinations, $Y\in L^2(\Omega,\mathcal{F},\mathbb{P})$. [/step]

[step:Use the non-negativity of the second moment] Since $Y^2\ge 0$ pointwise on $\Omega$ and $Y\in L^2(\Omega,\mathcal{F},\mathbb{P})$, the expectation $\mathbb{E}[Y^2]$ is finite and satisfies \begin{align*} \mathbb{E}[Y^2]\ge 0. \end{align*} [/step]

[step:Expand the square into the covariance quadratic form]Using finite distributivity and linearity of expectation, we compute \begin{align*} \mathbb{E}[Y^2] &=\mathbb{E}\left[\left(\sum_{i=1}^{n}a_i Z_i\right)\left(\sum_{j=1}^{n}a_j Z_j\right)\right]\\ &=\mathbb{E}\left[\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j Z_i Z_j\right]\\ &=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\,\mathbb{E}[Z_i Z_j]. \end{align*} For each pair $i,j\in\{1,\dots,n\}$, the definition of $Z_i$ and $Z_j$ gives \begin{align*} \mathbb{E}[Z_i Z_j] =\mathbb{E}[(X_i-\mu)(X_j-\mu)]. \end{align*}[/step]

[guided]The point of introducing $Y$ is that its square automatically produces all pairwise covariance terms. We first expand the square as a finite algebraic identity: \begin{align*} Y^2 &=\left(\sum_{i=1}^{n}a_i Z_i\right)\left(\sum_{j=1}^{n}a_j Z_j\right)\\ &=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j Z_i Z_j. \end{align*} There is no convergence issue here because both sums are finite. Since each $Z_i\in L^2(\Omega,\mathcal{F},\mathbb{P})$, the products $Z_iZ_j$ are integrable by the Cauchy-Schwarz inequality for random variables: \begin{align*} \mathbb{E}[|Z_iZ_j|]\le \mathbb{E}[Z_i^2]^{1/2}\mathbb{E}[Z_j^2]^{1/2}<\infty. \end{align*} Therefore linearity of expectation applies to the finite sum, giving \begin{align*} \mathbb{E}[Y^2] &=\mathbb{E}\left[\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j Z_i Z_j\right]\\ &=\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\,\mathbb{E}[Z_i Z_j]. \end{align*} Finally, because $Z_i=X_i-\mu$ and $Z_j=X_j-\mu$, each expectation in the last expression is exactly \begin{align*} \mathbb{E}[Z_i Z_j] =\mathbb{E}[(X_i-\mu)(X_j-\mu)]. \end{align*}[/guided]

[step:Identify each covariance term with the autocovariance function] By weak stationarity, the covariance between $X_i$ and $X_j$ depends only on the lag $i-j$. Since $\gamma(h)=\mathbb{E}[(X_{t+h}-\mu)(X_t-\mu)]$ for every $h,t\in\mathbb{Z}$, taking $h=i-j$ and $t=j$ gives \begin{align*} \mathbb{E}[(X_i-\mu)(X_j-\mu)] =\mathbb{E}[(X_{j+(i-j)}-\mu)(X_j-\mu)] =\gamma(i-j). \end{align*} Substituting this identity into the expansion of $\mathbb{E}[Y^2]$ yields \begin{align*} \mathbb{E}[Y^2] =\sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\gamma(i-j). \end{align*} Combining this equality with $\mathbb{E}[Y^2]\ge0$ proves \begin{align*} \sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j\gamma(i-j)\ge0. \end{align*} This is the desired positive definiteness of the autocovariance function. [/step]