[proofplan]
We first define the infinite linear process as an $L^2$ limit of finite partial sums and use absolute summability of the coefficients to prove convergence. The mean is obtained by passing the zero mean of finite partial sums to the $L^2$ limit. The covariance formula is first computed for finite coefficient sets, where the white-noise covariance eliminates all non-matching innovation terms. Finally, continuity of the covariance pairing in $L^2$ permits passage to the infinite sums, and the resulting expression depends only on the lag.
[/proofplan]
[step:Construct the linear process as an $L^2$ limit]
For a finite set $F \subset \mathbb{Z}$ and $t \in \mathbb{Z}$, define the finite partial sum
\begin{align*}
S_{F,t}: \Omega &\to \mathbb{R} \\
\omega &\mapsto \sum_{j \in F} \psi_j Z_{t-j}(\omega).
\end{align*}
Since each $Z_{t-j}$ belongs to $L^2(\Omega,\mathcal{F},\mathbb{P})$, the finite sum $S_{F,t}$ belongs to $L^2(\Omega,\mathcal{F},\mathbb{P})$.
Let $F,G \subset \mathbb{Z}$ be finite. Using the triangle inequality in $L^2(\Omega,\mathcal{F},\mathbb{P})$ and the identity $\|Z_k\|_{L^2} = \sigma$ for every $k \in \mathbb{Z}$, we obtain
\begin{align*}
\|S_{F,t} - S_{G,t}\|_{L^2}
&= \left\|\sum_{j \in F \triangle G} \varepsilon_j \psi_j Z_{t-j}\right\|_{L^2} \\
&\leq \sum_{j \in F \triangle G} |\psi_j|\,\|Z_{t-j}\|_{L^2} \\
&= \sigma \sum_{j \in F \triangle G} |\psi_j|,
\end{align*}
where $\varepsilon_j \in \{-1,1\}$ records whether $j$ belongs to $F \setminus G$ or to $G \setminus F$. Since $\sum_{j \in \mathbb{Z}} |\psi_j| < \infty$, the net $(S_{F,t})_{F \subset \mathbb{Z},\, F \text{ finite}}$ is Cauchy in $L^2(\Omega,\mathcal{F},\mathbb{P})$. Because $L^2(\Omega,\mathcal{F},\mathbb{P})$ is complete, there exists $X_t \in L^2(\Omega,\mathcal{F},\mathbb{P})$ such that
\begin{align*}
S_{F,t} \to X_t \quad \text{in } L^2(\Omega,\mathcal{F},\mathbb{P})
\end{align*}
as $F$ increases through the finite subsets of $\mathbb{Z}$. This defines
\begin{align*}
X_t = \sum_{j \in \mathbb{Z}} \psi_j Z_{t-j}
\end{align*}
as an $L^2$ limit.
[guided]
The infinite sum defining $X_t$ is not a pointwise sum at the start of the proof; it must first be constructed in a normed space. For each finite set $F \subset \mathbb{Z}$ and each time $t \in \mathbb{Z}$, define
\begin{align*}
S_{F,t}: \Omega &\to \mathbb{R} \\
\omega &\mapsto \sum_{j \in F} \psi_j Z_{t-j}(\omega).
\end{align*}
This is a finite linear combination of $L^2$ random variables, hence $S_{F,t} \in L^2(\Omega,\mathcal{F},\mathbb{P})$.
To show that these finite sums converge as $F$ exhausts $\mathbb{Z}$, take two finite sets $F,G \subset \mathbb{Z}$. The difference only contains indices in the symmetric difference $F \triangle G$, so for suitable signs $\varepsilon_j \in \{-1,1\}$,
\begin{align*}
S_{F,t} - S_{G,t} = \sum_{j \in F \triangle G} \varepsilon_j \psi_j Z_{t-j}.
\end{align*}
The triangle inequality in $L^2(\Omega,\mathcal{F},\mathbb{P})$ gives
\begin{align*}
\|S_{F,t} - S_{G,t}\|_{L^2}
&\leq \sum_{j \in F \triangle G} |\psi_j|\,\|Z_{t-j}\|_{L^2}.
\end{align*}
Since $\mathbb{E}[Z_{t-j}] = 0$ and $\operatorname{Var}(Z_{t-j}) = \sigma^2$, we have
\begin{align*}
\|Z_{t-j}\|_{L^2}^2
= \mathbb{E}[Z_{t-j}^2]
= \operatorname{Var}(Z_{t-j})
= \sigma^2.
\end{align*}
Therefore
\begin{align*}
\|S_{F,t} - S_{G,t}\|_{L^2}
\leq \sigma \sum_{j \in F \triangle G} |\psi_j|.
\end{align*}
The absolute summability of $(\psi_j)_{j \in \mathbb{Z}}$ implies that the tails on the right-hand side tend to $0$ as $F$ and $G$ exhaust $\mathbb{Z}$. Thus $(S_{F,t})$ is Cauchy in $L^2(\Omega,\mathcal{F},\mathbb{P})$. Completeness of $L^2(\Omega,\mathcal{F},\mathbb{P})$ gives a limit $X_t \in L^2(\Omega,\mathcal{F},\mathbb{P})$, and this limit is precisely the meaning of
\begin{align*}
X_t = \sum_{j \in \mathbb{Z}} \psi_j Z_{t-j}.
\end{align*}
[/guided]
[/step]
[step:Pass the zero mean from finite sums to the limit]
For every finite $F \subset \mathbb{Z}$ and every $t \in \mathbb{Z}$,
\begin{align*}
\mathbb{E}[S_{F,t}]
= \sum_{j \in F} \psi_j \mathbb{E}[Z_{t-j}]
= 0.
\end{align*}
Because $S_{F,t} \to X_t$ in $L^2(\Omega,\mathcal{F},\mathbb{P})$, we also have
\begin{align*}
|\mathbb{E}[X_t] - \mathbb{E}[S_{F,t}]|
= |\mathbb{E}[X_t - S_{F,t}]|
\leq \|X_t - S_{F,t}\|_{L^2} \to 0.
\end{align*}
Hence $\mathbb{E}[X_t] = 0$ for every $t \in \mathbb{Z}$.
[/step]
[step:Compute the covariance for finite coefficient sets]
Let $F,G \subset \mathbb{Z}$ be finite sets, and define
\begin{align*}
S_{F,t+h} &= \sum_{i \in F} \psi_i Z_{t+h-i}, \\
S_{G,t} &= \sum_{j \in G} \psi_j Z_{t-j}.
\end{align*}
Since both finite sums have mean $0$, their covariance is their second moment:
\begin{align*}
\operatorname{Cov}(S_{F,t+h},S_{G,t})
&= \mathbb{E}[S_{F,t+h}S_{G,t}] \\
&= \sum_{i \in F}\sum_{j \in G} \psi_i\psi_j\,\mathbb{E}[Z_{t+h-i}Z_{t-j}].
\end{align*}
For the white-noise process, $\mathbb{E}[Z_aZ_b] = \operatorname{Cov}(Z_a,Z_b)$ because $\mathbb{E}[Z_a] = \mathbb{E}[Z_b] = 0$. Therefore
\begin{align*}
\mathbb{E}[Z_{t+h-i}Z_{t-j}]
=
\begin{cases}
\sigma^2, & t+h-i = t-j,\\
0, & t+h-i \neq t-j.
\end{cases}
\end{align*}
The equality $t+h-i = t-j$ is equivalent to $i = j+h$. Hence
\begin{align*}
\operatorname{Cov}(S_{F,t+h},S_{G,t})
= \sigma^2 \sum_{\substack{j \in G \\ j+h \in F}} \psi_{j+h}\psi_j.
\end{align*}
[guided]
The finite computation is where the convolution structure appears. We take two finite coefficient sets $F,G \subset \mathbb{Z}$ and write
\begin{align*}
S_{F,t+h} &= \sum_{i \in F} \psi_i Z_{t+h-i}, \\
S_{G,t} &= \sum_{j \in G} \psi_j Z_{t-j}.
\end{align*}
Both sums have mean $0$, because each innovation has mean $0$. Therefore the covariance equals the expectation of the product:
\begin{align*}
\operatorname{Cov}(S_{F,t+h},S_{G,t})
= \mathbb{E}[S_{F,t+h}S_{G,t}].
\end{align*}
Since both sums are finite, we may expand the product term by term:
\begin{align*}
\mathbb{E}[S_{F,t+h}S_{G,t}]
= \sum_{i \in F}\sum_{j \in G} \psi_i\psi_j\,\mathbb{E}[Z_{t+h-i}Z_{t-j}].
\end{align*}
Now use the white-noise covariance condition. Because all innovations have mean $0$,
\begin{align*}
\mathbb{E}[Z_aZ_b] = \operatorname{Cov}(Z_a,Z_b).
\end{align*}
Thus
\begin{align*}
\mathbb{E}[Z_{t+h-i}Z_{t-j}]
=
\begin{cases}
\sigma^2, & t+h-i = t-j,\\
0, & t+h-i \neq t-j.
\end{cases}
\end{align*}
The only surviving terms are the matching innovation terms. Solving the equality
\begin{align*}
t+h-i = t-j
\end{align*}
gives $i = j+h$. Therefore every surviving term has coefficient $\psi_{j+h}\psi_j$, and the finite covariance is
\begin{align*}
\operatorname{Cov}(S_{F,t+h},S_{G,t})
= \sigma^2 \sum_{\substack{j \in G \\ j+h \in F}} \psi_{j+h}\psi_j.
\end{align*}
This formula is already independent of $t$; the time variable cancels in the matching condition.
[/guided]
[/step]
[step:Pass the finite covariance formula to the infinite process]
Fix $h \in \mathbb{Z}$ and $t \in \mathbb{Z}$. Let $(F_n)_{n \in \mathbb{N}}$ be an increasing sequence of finite subsets of $\mathbb{Z}$ with $\bigcup_{n=1}^{\infty} F_n = \mathbb{Z}$. Then
\begin{align*}
S_{F_n,t+h} \to X_{t+h}
\quad \text{and} \quad
S_{F_n,t} \to X_t
\end{align*}
in $L^2(\Omega,\mathcal{F},\mathbb{P})$.
For square-integrable random variables $Y,Y_n,W,W_n$ with $Y_n \to Y$ and $W_n \to W$ in $L^2$, the covariance pairing is continuous:
\begin{align*}
|\mathbb{E}[Y_nW_n] - \mathbb{E}[YW]|
&\leq |\mathbb{E}[(Y_n-Y)W_n]| + |\mathbb{E}[Y(W_n-W)]| \\
&\leq \|Y_n-Y\|_{L^2}\|W_n\|_{L^2} + \|Y\|_{L^2}\|W_n-W\|_{L^2}.
\end{align*}
The sequence $(W_n)$ is bounded in $L^2$ because it converges in $L^2$, so the right-hand side tends to $0$. Applying this with $Y_n=S_{F_n,t+h}$, $Y=X_{t+h}$, $W_n=S_{F_n,t}$, and $W=X_t$, and using the already proved zero means, gives
\begin{align*}
\operatorname{Cov}(X_{t+h},X_t)
= \lim_{n \to \infty} \operatorname{Cov}(S_{F_n,t+h},S_{F_n,t}).
\end{align*}
By the finite covariance formula,
\begin{align*}
\operatorname{Cov}(S_{F_n,t+h},S_{F_n,t})
= \sigma^2 \sum_{\substack{j \in F_n \\ j+h \in F_n}} \psi_{j+h}\psi_j.
\end{align*}
The series $\sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j$ is absolutely summable, since $(\psi_j)_{j \in \mathbb{Z}} \in \ell^1(\mathbb{Z})$ is bounded and
\begin{align*}
\sum_{j \in \mathbb{Z}} |\psi_{j+h}\psi_j|
\leq \left(\sup_{k \in \mathbb{Z}} |\psi_k|\right)\sum_{j \in \mathbb{Z}} |\psi_j|
< \infty.
\end{align*}
Therefore
\begin{align*}
\lim_{n \to \infty}
\sum_{\substack{j \in F_n \\ j+h \in F_n}} \psi_{j+h}\psi_j
=
\sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j.
\end{align*}
Consequently,
\begin{align*}
\operatorname{Cov}(X_{t+h},X_t)
= \sigma^2 \sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j.
\end{align*}
[guided]
We now replace the finite coefficient sets by all of $\mathbb{Z}$. Fix a lag $h \in \mathbb{Z}$ and a time $t \in \mathbb{Z}$. Choose an increasing sequence of finite sets $(F_n)_{n \in \mathbb{N}}$ such that $\bigcup_{n=1}^{\infty} F_n = \mathbb{Z}$. From the construction of $X_t$,
\begin{align*}
S_{F_n,t+h} \to X_{t+h}
\quad \text{and} \quad
S_{F_n,t} \to X_t
\end{align*}
in $L^2(\Omega,\mathcal{F},\mathbb{P})$.
The operation of taking covariance is continuous under $L^2$ convergence. To see this directly, let $Y_n \to Y$ and $W_n \to W$ in $L^2(\Omega,\mathcal{F},\mathbb{P})$. Then
\begin{align*}
|\mathbb{E}[Y_nW_n] - \mathbb{E}[YW]|
&\leq |\mathbb{E}[(Y_n-Y)W_n]| + |\mathbb{E}[Y(W_n-W)]| \\
&\leq \|Y_n-Y\|_{L^2}\|W_n\|_{L^2} + \|Y\|_{L^2}\|W_n-W\|_{L^2}.
\end{align*}
The first factor $\|Y_n-Y\|_{L^2}$ tends to $0$, the second sequence $\|W_n\|_{L^2}$ is bounded because $W_n$ converges in $L^2$, and the final factor $\|W_n-W\|_{L^2}$ tends to $0$. Hence $\mathbb{E}[Y_nW_n] \to \mathbb{E}[YW]$.
Applying this to $Y_n=S_{F_n,t+h}$, $Y=X_{t+h}$, $W_n=S_{F_n,t}$, and $W=X_t$, and using that all these variables have mean tending to $0$, gives
\begin{align*}
\operatorname{Cov}(X_{t+h},X_t)
= \lim_{n \to \infty} \operatorname{Cov}(S_{F_n,t+h},S_{F_n,t}).
\end{align*}
For each $n$, the finite computation gives
\begin{align*}
\operatorname{Cov}(S_{F_n,t+h},S_{F_n,t})
= \sigma^2 \sum_{\substack{j \in F_n \\ j+h \in F_n}} \psi_{j+h}\psi_j.
\end{align*}
It remains only to justify that these finite sums converge to the full series. Since $\sum_{j \in \mathbb{Z}}|\psi_j| < \infty$, the sequence $(\psi_j)_{j \in \mathbb{Z}}$ is bounded; define
\begin{align*}
M := \sup_{k \in \mathbb{Z}} |\psi_k| < \infty.
\end{align*}
Then
\begin{align*}
\sum_{j \in \mathbb{Z}} |\psi_{j+h}\psi_j|
\leq M \sum_{j \in \mathbb{Z}} |\psi_j|
< \infty.
\end{align*}
Thus the covariance series is absolutely summable, and exhausting $\mathbb{Z}$ by finite sets gives
\begin{align*}
\lim_{n \to \infty}
\sum_{\substack{j \in F_n \\ j+h \in F_n}} \psi_{j+h}\psi_j
=
\sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j.
\end{align*}
Combining the previous limits yields
\begin{align*}
\operatorname{Cov}(X_{t+h},X_t)
= \sigma^2 \sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j.
\end{align*}
[/guided]
[/step]
[step:Conclude weak stationarity]
We have shown that $\mathbb{E}[X_t] = 0$ for every $t \in \mathbb{Z}$. For every $h,t \in \mathbb{Z}$,
\begin{align*}
\operatorname{Cov}(X_{t+h},X_t)
= \sigma^2 \sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j,
\end{align*}
and the right-hand side depends only on $h$. Therefore $(X_t)_{t \in \mathbb{Z}}$ has constant mean and autocovariance depending only on the lag. Hence $(X_t)_{t \in \mathbb{Z}}$ is weakly stationary, with autocovariance function
\begin{align*}
\gamma_X(h) = \sigma^2 \sum_{j \in \mathbb{Z}} \psi_{j+h}\psi_j.
\end{align*}
[/step]
