[guided]The power series was chosen so that it is the reciprocal of $\psi$. Using the notation $\psi(z)=\sum_{k=0}^{p}\psi_k z^k$ fixed above, the reciprocal identity is
\begin{align*}
\psi(z)\sum_{j=0}^{\infty}\pi_j z^j = 1
\end{align*}
for $|z|<r$. Because $\psi$ is a finite polynomial and $(\pi_j)_{j=0}^{\infty}\in\ell^1$, multiplying the finite polynomial into the absolutely convergent series is legitimate. Comparing Taylor coefficients gives
\begin{align*}
\sum_{k=0}^{\min\{p,m\}}\psi_k\pi_{m-k}
=
\begin{cases}
1, & m=0,\\
0, & m\geq 1,
\end{cases}
\end{align*}
where $\pi_j:=0$ for $j<0$.
We now translate this coefficient identity back into the stochastic equation using finite truncations. For $N\in\mathbb{N}$, set
\begin{align*}
\widetilde{Y}_{t,N}:=\sum_{j=0}^{N}\pi_jZ_{t-j}.
\end{align*}
The construction above gives $\widetilde{Y}_{t,N}\to\widetilde{Y}_t$ in $L^2(\Omega,\mathcal{F},\mathbb{P})$ for every $t$. Since $\psi(B)=\sum_{k=0}^{p}\psi_kB^k$ is a finite sum of shifts, applying $\psi(B)$ preserves $L^2$ convergence, so
\begin{align*}
\psi(B)\widetilde{Y}_{t,N}\to\psi(B)\widetilde{Y}_t
\end{align*}
in $L^2(\Omega,\mathcal{F},\mathbb{P})$. For each fixed $N$, all sums are finite, and therefore
\begin{align*}
\psi(B)\widetilde{Y}_{t,N}
&= \sum_{m=0}^{N+p}\left(\sum_{k=0}^{\min\{p,m\}}\psi_k\pi_{m-k}\right)Z_{t-m},
\end{align*}
where $\pi_j:=0$ for $j<0$ and for $j>N$ in this finite truncation. The coefficient identity kills all coefficients except the coefficient of $Z_t$, but the finite truncation creates boundary terms near the cutoff. More precisely,
\begin{align*}
\psi(B)\widetilde{Y}_{t,N}=Z_t+R_{t,N},
\end{align*}
where
\begin{align*}
R_{t,N}:=\sum_{m=N+1}^{N+p}\left(\sum_{k=m-N}^{p}\psi_k\pi_{m-k}\right)Z_{t-m}.
\end{align*}
The tail $R_{t,N}$ vanishes in $L^2$. Indeed, orthogonality of the white-noise variables and the triangle inequality give
\begin{align*}
\mathbb{E}[|R_{t,N}|^2]
&=\sigma^2\sum_{m=N+1}^{N+p}\left|\sum_{k=m-N}^{p}\psi_k\pi_{m-k}\right|^2 \\
&\leq \sigma^2\left(\sum_{m=N+1}^{N+p}\sum_{k=m-N}^{p}|\psi_k|\,|\pi_{m-k}|\right)^2 \\
&\leq \sigma^2\left(\sum_{k=0}^{p}|\psi_k|\right)^2\left(\sum_{j=N+1-p}^{N}|\pi_j|\right)^2.
\end{align*}
Because $(\pi_j)_{j=0}^{\infty}\in\ell^1$, the final tail tends to $0$ as $N\to\infty$. Since $\psi(B)\widetilde{Y}_{t,N}\to\psi(B)\widetilde{Y}_t$ in $L^2(\Omega,\mathcal{F},\mathbb{P})$, we conclude
\begin{align*}
\psi(B)\widetilde{Y}_t=Z_t.
\end{align*}
This is the precise reason the infinite filter may be convolved with the finite autoregressive polynomial: the convolution is first performed at the finite level, the boundary tail is estimated, and only then the $L^2$ limit is taken.
It remains to explain uniqueness in the stated causal absolutely summable linear-process class without appealing to an external criterion. Let $(V_t)_{t\in\mathbb{Z}}$ be another such solution. Then, by definition of this solution class, there is a sequence $(a_j)_{j=0}^{\infty}\in\ell^1$ such that
\begin{align*}
V_t=\sum_{j=0}^{\infty}a_j Z_{t-j}
\end{align*}
in $L^2(\Omega,\mathcal{F},\mathbb{P})$ and $\psi(B)V_t=Z_t$. There is one degenerate case to separate. If $\sigma^2=0$, then $Z_t=0$ in $L^2(\Omega,\mathcal{F},\mathbb{P})$ for every $t$, so every causal linear process driven by $(Z_t)$ is the zero process in $L^2$. In that case $V_t=\widetilde{Y}_t=0$ for every $t$.
Assume now that $\sigma^2>0$. Define
\begin{align*}
c_m:=\sum_{k=0}^{\min\{p,m\}}\psi_k a_{m-k},
\end{align*}
with $a_j:=0$ for $j<0$. The sequence $(c_m)_{m=0}^{\infty}$ is absolutely summable because it is the convolution of the finite sequence $(\psi_0,\dots,\psi_p)$ with the $\ell^1$ sequence $(a_j)_{j=0}^{\infty}$. Hence applying $\psi(B)$ to $V$ gives
\begin{align*}
\psi(B)V_t=\sum_{m=0}^{\infty}c_mZ_{t-m}
\end{align*}
in $L^2(\Omega,\mathcal{F},\mathbb{P})$. Since $\psi(B)V_t=Z_t$, define $d_0:=c_0-1$ and $d_m:=c_m$ for $m\geq1$; then
\begin{align*}
\sum_{m=0}^{\infty}d_mZ_{t-m}=0
\end{align*}
in $L^2(\Omega,\mathcal{F},\mathbb{P})$. To identify the coefficients, fix $\ell\geq0$ and take covariance with $Z_{t-\ell}$. Continuity of covariance under $L^2$ convergence justifies passing through the infinite sum, and the white-noise covariance relation yields
\begin{align*}
0=\operatorname{Cov}\left(\sum_{m=0}^{\infty}d_m Z_{t-m},Z_{t-\ell}\right)=\sigma^2 d_\ell.
\end{align*}
Because $\sigma^2>0$, this gives $d_\ell=0$. Since $\ell$ was arbitrary, $c_0=1$ and $c_m=0$ for every $m\geq1$. Therefore
\begin{align*}
\sum_{m=0}^{\infty}\left(\sum_{k=0}^{\min\{p,m\}}\psi_k a_{m-k}\right)z^m=1
\end{align*}
for $|z|<1$. Since the root condition gives $\psi(z)\neq0$ on $|z|<1$, division by $\psi(z)$ is valid there, and hence
\begin{align*}
\sum_{m=0}^{\infty}a_m z^m=\frac{1}{\psi(z)}=\sum_{m=0}^{\infty}\pi_m z^m
\end{align*}
for $|z|<1$. Uniqueness of Taylor coefficients now forces $a_m=\pi_m$ for every $m\geq0$. Thus any causal absolutely summable linear-process solution agrees with $\widetilde{Y}$ in $L^2(\Omega,\mathcal{F},\mathbb{P})$ at every time.
We have therefore verified the reduced ARMA description: the equation $\psi(B)Y_t=Z_t$ has the unique causal absolutely summable linear-process solution
\begin{align*}
\widetilde{Y}_t = \sum_{j=0}^{\infty} \pi_j Z_{t-j}.
\end{align*}
The covariance computation above shows that this solution is weakly stationary. The first step proved separately that the algebraic differencing operation sends $\phi(B)X_t=Z_t$ to $\psi(B)\Delta X_t=Z_t$. Thus differencing cancels the simple unit-root factor $1-B$, and, in the stated causal linear-process class, the stationary reduced dynamics are governed by $\psi(B)$.[/guided]