[proofplan] We prove necessity by applying non-negativity of second moments to a finite complex linear combination of centred process values. For sufficiency, positive definiteness makes every finite Toeplitz matrix $(\gamma(t_i-t_j))$ a valid covariance matrix, so we define consistent centred Gaussian laws on all finite coordinate sets. The Kolmogorov extension theorem then produces a process on $\mathbb{R}^{\mathbb{Z}}$, and its finite-dimensional Gaussian laws give the prescribed autocovariance. [/proofplan] [step:Derive positive definiteness from non-negativity of a centred square] Assume that $\gamma$ is the autocovariance function of a real-valued weakly stationary process. Thus there are a probability space $(\Omega,\mathcal{F},\mathbb{P})$, measurable maps $X_t:(\Omega,\mathcal{F}) \to (\mathbb{R},\mathcal{B}(\mathbb{R}))$ for $t \in \mathbb{Z}$, and a constant $\mu \in \mathbb{R}$ such that $\mathbb{E}[X_t]=\mu$, $\mathbb{E}[X_t^2]<\infty$, and \begin{align*} \gamma(h)=\mathbb{E}[(X_{t+h}-\mu)(X_t-\mu)] \end{align*} for all $t,h \in \mathbb{Z}$. Fix $n \in \mathbb{N}$, indices $t_1,\dots,t_n \in \mathbb{Z}$, and coefficients $c_1,\dots,c_n \in \mathbb{C}$. Define the complex-valued random variable $Y:(\Omega,\mathcal{F}) \to (\mathbb{C},\mathcal{B}(\mathbb{C}))$ by \begin{align*} Y(\omega)=\sum_{i=1}^{n} c_i\bigl(X_{t_i}(\omega)-\mu\bigr). \end{align*} Since each $X_{t_i}$ has finite second moment, $Y \in L^2(\Omega;\mathbb{C})$. Therefore $\mathbb{E}[|Y|^2]\ge 0$. Expanding the product and using linearity of expectation for finite sums gives \begin{align*} 0 &\le \mathbb{E}[|Y|^2] \\ &= \mathbb{E}\left[\left(\sum_{i=1}^{n} c_i(X_{t_i}-\mu)\right) \overline{\left(\sum_{j=1}^{n} c_j(X_{t_j}-\mu)\right)}\right] \\ &= \sum_{i=1}^{n}\sum_{j=1}^{n} c_i\overline{c_j}\, \mathbb{E}[(X_{t_i}-\mu)(X_{t_j}-\mu)] \\ &= \sum_{i=1}^{n}\sum_{j=1}^{n} c_i\overline{c_j}\,\gamma(t_i-t_j). \end{align*} Hence $\gamma$ is positive definite. Also \begin{align*} \gamma(0)=\mathbb{E}[(X_0-\mu)^2]<\infty, \end{align*} because $X_0$ has finite second moment. [guided] Assume first that $\gamma$ really comes from a weakly stationary process. We must show that every finite quadratic form built from $\gamma$ is non-negative. Fix $n \in \mathbb{N}$, times $t_1,\dots,t_n \in \mathbb{Z}$, and complex coefficients $c_1,\dots,c_n \in \mathbb{C}$. The natural object to test is the finite centred linear combination \begin{align*} Y:(\Omega,\mathcal{F}) &\to (\mathbb{C},\mathcal{B}(\mathbb{C})) \\ \omega &\mapsto \sum_{i=1}^{n} c_i\bigl(X_{t_i}(\omega)-\mu\bigr). \end{align*} This random variable is square-integrable because it is a finite linear combination of square-integrable random variables. The quantity $\mathbb{E}[|Y|^2]$ is non-negative, and its expansion is exactly the positive-definiteness expression: \begin{align*} 0 &\le \mathbb{E}[|Y|^2] \\ &= \mathbb{E}\left[\left(\sum_{i=1}^{n} c_i(X_{t_i}-\mu)\right) \overline{\left(\sum_{j=1}^{n} c_j(X_{t_j}-\mu)\right)}\right] \\ &= \sum_{i=1}^{n}\sum_{j=1}^{n} c_i\overline{c_j}\, \mathbb{E}[(X_{t_i}-\mu)(X_{t_j}-\mu)]. \end{align*} Weak stationarity says that the covariance depends only on the time difference. Since \begin{align*} \mathbb{E}[(X_{t_i}-\mu)(X_{t_j}-\mu)] = \gamma(t_i-t_j), \end{align*} we obtain \begin{align*} \sum_{i=1}^{n}\sum_{j=1}^{n} c_i\overline{c_j}\,\gamma(t_i-t_j) \ge 0. \end{align*} This is precisely positive definiteness. Finally, $\gamma(0)$ is the variance at any time, for example \begin{align*} \gamma(0)=\mathbb{E}[(X_0-\mu)^2], \end{align*} and this is finite by the finite-second-moment requirement in weak stationarity. [/guided] [/step] [step:Build finite-dimensional Gaussian laws from the positive definite matrices] Conversely, assume that $\gamma:\mathbb{Z}\to\mathbb{R}$ is positive definite and $\gamma(0)<\infty$. Let $F=\{s_1,\dots,s_m\}\subset \mathbb{Z}$ be a finite set with a fixed ordering. Define the matrix $C_F \in \mathbb{R}^{m\times m}$ by \begin{align*} (C_F)_{ij}=\gamma(s_i-s_j). \end{align*} Positive definiteness gives \begin{align*} \sum_{i=1}^{m}\sum_{j=1}^{m} c_i\overline{c_j}(C_F)_{ij}\ge 0 \end{align*} for all $c_1,\dots,c_m \in \mathbb{C}$. Since $\gamma$ is real-valued, this also implies $(C_F)_{ij}=(C_F)_{ji}$: otherwise the Hermitian condition forced by non-negativity of all complex quadratic forms would fail. Thus $C_F$ is a real symmetric positive semidefinite matrix. By the finite-dimensional Gaussian construction, there exists a centred Gaussian probability measure $\nu_F$ on $(\mathbb{R}^{F},\mathcal{B}(\mathbb{R}^{F}))$ with covariance matrix $C_F$. Equivalently, for the coordinate maps $x_s:\mathbb{R}^{F}\to\mathbb{R}$, $s\in F$, this measure satisfies \begin{align*} \int_{\mathbb{R}^{F}} x_{s_i}\,d\nu_F(x)&=0,\\ \int_{\mathbb{R}^{F}} x_{s_i}x_{s_j}\,d\nu_F(x)&=\gamma(s_i-s_j). \end{align*} [guided] Now assume that $\gamma$ is positive definite and $\gamma(0)<\infty$. The goal is to prescribe the joint distribution of finitely many random variables first. Let $F=\{s_1,\dots,s_m\}$ be a finite subset of $\mathbb{Z}$ with a chosen ordering, and define \begin{align*} C_F \in \mathbb{R}^{m\times m}, \qquad (C_F)_{ij}=\gamma(s_i-s_j). \end{align*} This is the only possible covariance matrix for the vector $(X_{s_1},\dots,X_{s_m})$ if the desired autocovariance is $\gamma$. We must verify that $C_F$ is a valid covariance matrix. Positive definiteness says that for every $c=(c_1,\dots,c_m)\in\mathbb{C}^m$, \begin{align*} \sum_{i=1}^{m}\sum_{j=1}^{m} c_i\overline{c_j}(C_F)_{ij}\ge 0. \end{align*} A matrix whose complex quadratic form is always real and non-negative is Hermitian. Since the entries of $C_F$ are real, Hermitian means symmetric, so $C_F$ is a real symmetric positive semidefinite matrix. The diagonal entries are finite because \begin{align*} (C_F)_{ii}=\gamma(0)<\infty. \end{align*} A real symmetric positive semidefinite matrix is a valid covariance matrix for a centred Gaussian vector. Therefore there is a centred Gaussian probability measure $\nu_F$ on $(\mathbb{R}^{F},\mathcal{B}(\mathbb{R}^{F}))$ such that, for the coordinate maps $x_s:\mathbb{R}^{F}\to\mathbb{R}$, \begin{align*} \int_{\mathbb{R}^{F}} x_{s_i}\,d\nu_F(x)&=0,\\ \int_{\mathbb{R}^{F}} x_{s_i}x_{s_j}\,d\nu_F(x)&=\gamma(s_i-s_j). \end{align*} Thus positive definiteness is exactly the condition that makes all the finite covariance prescriptions legitimate. [/guided] [/step] [step:Verify consistency of the finite-dimensional laws] Let $E\subset F$ be finite subsets of $\mathbb{Z}$, and let $\pi_{F,E}:\mathbb{R}^{F}\to\mathbb{R}^{E}$ be the coordinate projection. We claim that \begin{align*} (\pi_{F,E})_{\#}\nu_F=\nu_E. \end{align*} Indeed, the pushforward $(\pi_{F,E})_{\#}\nu_F$ is a centred Gaussian measure on $\mathbb{R}^{E}$. Its covariance between the coordinates indexed by $r,u\in E$ is the same covariance computed under $\nu_F$, namely $\gamma(r-u)$. This is exactly the covariance matrix defining $\nu_E$. Since centred Gaussian measures on finite-dimensional Euclidean spaces are determined by their covariance matrices, the two measures are equal. [guided] The finite-dimensional laws must agree when we forget coordinates. Let $E\subset F$ be finite subsets of $\mathbb{Z}$, and define the coordinate projection \begin{align*} \pi_{F,E}:\mathbb{R}^{F} &\to \mathbb{R}^{E} \\ (x_s)_{s\in F} &\mapsto (x_s)_{s\in E}. \end{align*} We need to prove that the marginal of $\nu_F$ on the coordinates in $E$ is precisely $\nu_E$: \begin{align*} (\pi_{F,E})_{\#}\nu_F=\nu_E. \end{align*} The pushforward of a centred Gaussian measure under a linear map is again centred Gaussian. The map $\pi_{F,E}$ is linear, so $(\pi_{F,E})_{\#}\nu_F$ is a centred Gaussian measure on $\mathbb{R}^{E}$. Its covariance between coordinate functions $x_r$ and $x_u$, for $r,u\in E$, is computed using the same coordinates inside $\mathbb{R}^{F}$: \begin{align*} \int_{\mathbb{R}^{E}} x_r x_u\,d((\pi_{F,E})_{\#}\nu_F)(x) = \int_{\mathbb{R}^{F}} x_r x_u\,d\nu_F(x) = \gamma(r-u). \end{align*} But this is exactly the covariance entry used to define $\nu_E$. Since finite-dimensional centred Gaussian measures are uniquely determined by their covariance matrices, the marginal law agrees with $\nu_E$. This proves consistency of the family $(\nu_F)$. [/guided] [/step] [step:Apply Kolmogorov extension to obtain the process] By the Kolmogorov extension theorem (citing a result not yet in the wiki: Kolmogorov extension theorem), applied to the consistent family $(\nu_F)_{F\subset\mathbb{Z},\,F\text{ finite}}$, there exists a probability measure $\mathbb{P}$ on \begin{align*} \Omega=\mathbb{R}^{\mathbb{Z}}, \qquad \mathcal{F}=\mathcal{B}(\mathbb{R})^{\otimes \mathbb{Z}}, \end{align*} such that for every finite $F\subset\mathbb{Z}$ the coordinate projection $\pi_F:\Omega\to\mathbb{R}^{F}$ has law $\nu_F$. For each $t\in\mathbb{Z}$, define the coordinate random variable \begin{align*} X_t:(\Omega,\mathcal{F}) &\to (\mathbb{R},\mathcal{B}(\mathbb{R})) \\ \omega &\mapsto \omega_t. \end{align*} Taking $F=\{t\}$ gives $\mathbb{E}[X_t]=0$ and \begin{align*} \mathbb{E}[X_t^2]=\gamma(0)<\infty. \end{align*} Taking $F=\{s,t\}$ gives \begin{align*} \mathbb{E}[X_sX_t]=\gamma(s-t) \end{align*} for all $s,t\in\mathbb{Z}$. [guided] The consistency just proved is exactly the hypothesis needed for the Kolmogorov extension theorem. Applying that theorem to the finite-dimensional laws $(\nu_F)_{F\subset\mathbb{Z},\,F\text{ finite}}$ gives a probability measure $\mathbb{P}$ on the product space \begin{align*} \Omega=\mathbb{R}^{\mathbb{Z}}, \qquad \mathcal{F}=\mathcal{B}(\mathbb{R})^{\otimes \mathbb{Z}}, \end{align*} with the property that every finite coordinate projection has the prescribed law. More precisely, if $F\subset\mathbb{Z}$ is finite and \begin{align*} \pi_F:\Omega &\to \mathbb{R}^{F} \\ \omega &\mapsto (\omega_s)_{s\in F}, \end{align*} then $(\pi_F)_{\#}\mathbb{P}=\nu_F$. Now define the process by the coordinate maps: \begin{align*} X_t:(\Omega,\mathcal{F}) &\to (\mathbb{R},\mathcal{B}(\mathbb{R})) \\ \omega &\mapsto \omega_t. \end{align*} These maps are measurable by the definition of the product $\sigma$-algebra. Because the one-dimensional marginal of $X_t$ is the centred Gaussian law with variance $\gamma(0)$, we have \begin{align*} \mathbb{E}[X_t]=0, \qquad \mathbb{E}[X_t^2]=\gamma(0)<\infty. \end{align*} Because the two-dimensional marginal of $(X_s,X_t)$ is the centred Gaussian law with covariance matrix determined by $\gamma$, we also have \begin{align*} \mathbb{E}[X_sX_t]=\gamma(s-t) \end{align*} for all $s,t\in\mathbb{Z}$. [/guided] [/step] [step:Identify the constructed process as weakly stationary with autocovariance $\gamma$] The process $(X_t)_{t\in\mathbb{Z}}$ has constant mean $\mu=0$ and finite second moments. For every $t,h\in\mathbb{Z}$, \begin{align*} \mathbb{E}[(X_{t+h}-0)(X_t-0)] = \mathbb{E}[X_{t+h}X_t] = \gamma((t+h)-t) = \gamma(h). \end{align*} Thus the covariance depends only on the lag $h$, so $(X_t)_{t\in\mathbb{Z}}$ is weakly stationary and has autocovariance function $\gamma$. This proves the converse implication and completes the proof. [/step]