[guided]We first build the candidate process as a limit of finite moving averages, because finite moving averages are unambiguous random variables. For $n\ge0$ and $t\in\mathbb{Z}$, set
\begin{align*}
X_{t,n}=\sum_{j=0}^{n}\psi_j Z_{t-j}.
\end{align*}
To prove that these partial sums converge in $L^2$, take $m>n$. The difference is
\begin{align*}
X_{t,m}-X_{t,n}
=
\sum_{j=n+1}^{m}\psi_j Z_{t-j}.
\end{align*}
When we square and take expectation, all cross terms vanish because distinct white-noise variables are orthogonal:
\begin{align*}
\mathbb{E}\left[\left|X_{t,m}-X_{t,n}\right|^2\right]
=
\sigma^2\sum_{j=n+1}^{m}|\psi_j|^2.
\end{align*}
The coefficient sequence is absolutely summable, hence square summable, so the right-hand side tends to $0$ as $m,n\to\infty$. Thus the partial sums are Cauchy in $L^2$. By [Completeness of $L^2$](/page/L2%20Space), this Cauchy sequence has a limit in $L^2$, and we denote it by
\begin{align*}
X_t=\sum_{j=0}^{\infty}\psi_j Z_{t-j}.
\end{align*}
This representation is causal because $X_t$ is built only from present and past innovations $Z_t,Z_{t-1},Z_{t-2},\dots$. To check stationarity, compute the covariance through finite approximations and then pass to the limit. For each $n$, orthogonality gives
\begin{align*}
\mathbb{E}[X_{t+h,n}\overline{X_{t,n}}]
=
\sigma^2\sum_{\substack{0\le i,j\le n\\ t+h-i=t-j}}\psi_i\overline{\psi_j}.
\end{align*}
Because $X_{t+h,n}\to X_{t+h}$ and $X_{t,n}\to X_t$ in $L^2$, the Cauchy-Schwarz inequality justifies passing the covariance to the limit:
\begin{align*}
\mathbb{E}[X_{t+h,n}\overline{X_{t,n}}]
\to
\mathbb{E}[X_{t+h}\overline{X_t}].
\end{align*}
The finite coefficient sums converge absolutely since $(\psi_j)$ is absolutely summable. Therefore
\begin{align*}
\mathbb{E}[X_{t+h}\overline{X_t}]
=
\sigma^2\sum_{\substack{i,j\ge0\\ t+h-i=t-j}}\psi_i\overline{\psi_j}.
\end{align*}
The condition $t+h-i=t-j$ is equivalent to $j=i-h$, so the resulting sum depends only on $h$. Therefore $(X_t)$ is weakly stationary.[/guided]