[step:Define the causal process and verify weak stationarity]
For each $t\in\mathbb{Z}$, define the [random variable](/page/Random%20Variable)
\begin{align*}
X_t=\sum_{j=0}^{\infty}\psi_j Z_{t-j}.
\end{align*}
Let $L^2(\Omega,\mathcal F,\mathbb P)$ denote the [Hilbert space](/page/Hilbert%20Space) of real-valued square-integrable random variables, modulo equality almost surely, with norm $\|W\|_{L^2}=(\mathbb{E}[|W|^2])^{1/2}$. Let $\ell^1(\mathbb{N}_0)$ denote the space of absolutely summable real sequences indexed by $\mathbb{N}_0=\{0,1,2,\dots\}$, and let $\ell^2(\mathbb{N}_0)$ denote the space of square-summable real sequences indexed by $\mathbb{N}_0$. Since $\sum_{j=0}^{\infty}|\psi_j|<\infty$, the series converges in $L^2(\Omega,\mathcal F,\mathbb P)$. Indeed, for $m>n$,
\begin{align*}
\mathbb{E}\left[\left|\sum_{j=n+1}^{m}\psi_j Z_{t-j}\right|^2\right]=\sigma^2\sum_{j=n+1}^{m}|\psi_j|^2,
\end{align*}
using the orthogonality of distinct white-noise innovations. Since $\ell^1(\mathbb{N}_0)\subset \ell^2(\mathbb{N}_0)$, the right-hand side tends to $0$ as $m,n\to\infty$.
The mean is
\begin{align*}
\mathbb{E}[X_t]=\sum_{j=0}^{\infty}\psi_j\mathbb{E}[Z_{t-j}]=0.
\end{align*}
For $h\in\mathbb{Z}$, define the covariance function $\gamma:\mathbb{Z}\to\mathbb{R}$ by
\begin{align*}
\gamma(h)=\operatorname{Cov}(X_{t+h},X_t).
\end{align*}
The $L^2$ convergence of the defining series and orthogonality of white noise give
\begin{align*}
\gamma(h)=\sigma^2\sum_{\substack{j, k\ge 0, t+h-j=t-k}}\psi_j\psi_k.
\end{align*}
Equivalently,
\begin{align*}
\gamma(h)=\sigma^2\sum_{j=\max\{0,h\}}^{\infty}\psi_j\psi_{j-h}
\end{align*}
for $h\ge 0$, with the analogous expression for $h<0$. This depends only on $h$, not on $t$. Hence $(X_t)_{t\in\mathbb{Z}}$ is weakly stationary and causal.
[/step]