[proofplan]
The root condition places every zero of $\phi$ strictly outside the closed unit disc, so the rational function $\theta/\phi$ is analytic on a disc of radius larger than $1$. Its Taylor coefficients therefore decay geometrically, which gives absolute summability. The corresponding infinite moving-average series is then well-defined in $L^2$ and satisfies the ARMA equation by matching coefficients in the identity $\phi(z)\psi(z)=\theta(z)$. Finally, applying the absolutely summable inverse autoregressive filter to the difference of two causal weakly stationary solutions proves uniqueness of the causal representation.
[/proofplan]
[step:Construct the analytic quotient and its Taylor coefficients]
Let $R_\phi \in (1,\infty]$ denote the distance from $0$ to the nearest zero of $\phi$, with $R_\phi=\infty$ if $\phi$ has no zero. The hypothesis $\phi(z)\neq 0$ for $|z|\leq 1$ gives $R_\phi>1$. Choose $\rho$ with $1<\rho<R_\phi$. Then $\phi(z)\neq 0$ for $|z|\leq \rho$, so the function
\begin{align*}
\psi: \{z\in\mathbb{C}: |z|<\rho\} &\to \mathbb{C}\\
z &\mapsto \frac{\theta(z)}{\phi(z)}
\end{align*}
is analytic.
Let $(\psi_j)_{j\geq 0}$ be the Taylor coefficients of $\psi$ at $0$, so that
\begin{align*}
\psi(z)=\sum_{j=0}^{\infty}\psi_j z^j
\end{align*}
for $|z|<\rho$. Define
\begin{align*}
M_\rho:=\sup_{|z|=\rho}|\psi(z)|.
\end{align*}
By the [Cauchy integral formula](/page/Cauchy%20Integral%20Formula), applied to the analytic map $\psi$ on the circle $\{z\in\mathbb{C}: |z|=\rho\}$, the Taylor coefficient satisfies
\begin{align*}
\psi_j
=
\frac{1}{2\pi i}\oint_{|z|=\rho}\frac{\psi(z)}{z^{j+1}}\,dz.
\end{align*}
Estimating this contour integral by the length $2\pi\rho$ of the circle and the bound $|\psi(z)|\leq M_\rho$ on the contour gives, for every $j\geq 0$,
\begin{align*}
|\psi_j|\leq M_\rho \rho^{-j}.
\end{align*}
Therefore
\begin{align*}
\sum_{j=0}^{\infty}|\psi_j|
\leq
M_\rho \sum_{j=0}^{\infty}\rho^{-j}
=
M_\rho\frac{\rho}{\rho-1}
<\infty.
\end{align*}
Thus $(\psi_j)_{j\geq 0}$ is absolutely summable.
[guided]
The root condition is used first to create analytic room around the unit disc. Let $R_\phi$ be the modulus of the closest zero of $\phi$ to the origin, and set $R_\phi=\infty$ if $\phi$ has no zero. Since no zero lies in the closed unit disc, we have $R_\phi>1$. Pick a number $\rho$ satisfying $1<\rho<R_\phi$. On the closed disc $|z|\leq\rho$, the denominator $\phi(z)$ never vanishes, so the quotient
\begin{align*}
\psi: \{z\in\mathbb{C}: |z|<\rho\} &\to \mathbb{C}\\
z &\mapsto \frac{\theta(z)}{\phi(z)}
\end{align*}
is analytic.
Because $\psi$ is analytic at $0$, it has a Taylor expansion
\begin{align*}
\psi(z)=\sum_{j=0}^{\infty}\psi_j z^j
\end{align*}
for $|z|<\rho$. The key point is that the radius $\rho$ is larger than $1$, so the Taylor coefficients decay at least geometrically. Define the finite boundary bound
\begin{align*}
M_\rho:=\sup_{|z|=\rho}|\psi(z)|.
\end{align*}
The [Cauchy integral formula](/page/Cauchy%20Integral%20Formula) on the circle $|z|=\rho$ gives
\begin{align*}
\psi_j
=
\frac{1}{2\pi i}\oint_{|z|=\rho}\frac{\psi(z)}{z^{j+1}}\,dz.
\end{align*}
The contour has length $2\pi\rho$, and on the contour we have $|\psi(z)|\leq M_\rho$ and $|z|^{-(j+1)}=\rho^{-(j+1)}$. Hence
\begin{align*}
|\psi_j|
\leq
\frac{1}{2\pi}(2\pi\rho)M_\rho\rho^{-(j+1)}
=
M_\rho\rho^{-j}
\end{align*}
for every $j\geq 0$. Summing this geometric bound yields
\begin{align*}
\sum_{j=0}^{\infty}|\psi_j|
\leq
M_\rho \sum_{j=0}^{\infty}\rho^{-j}
=
M_\rho\frac{\rho}{\rho-1}
<\infty.
\end{align*}
This proves absolute summability of the impulse response coefficients.
[/guided]
[/step]
[step:Verify that the moving-average series is well-defined in $L^2$]
For each $t\in\mathbb{Z}$, define the partial sums
\begin{align*}
Y_{t,N}:=\sum_{j=0}^{N}\psi_j Z_{t-j},\qquad N\geq 0.
\end{align*}
Since $(Z_t)$ is white noise, the random variables $Z_{t-j}$ and $Z_{t-k}$ are uncorrelated whenever $j\neq k$. Hence, for $M>N$,
\begin{align*}
\mathbb{E}\left[\left|Y_{t,M}-Y_{t,N}\right|^2\right]
&=
\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*}
Because $\sum_{j=0}^{\infty}|\psi_j|<\infty$, also $\sum_{j=0}^{\infty}|\psi_j|^2<\infty$. Therefore $(Y_{t,N})_{N\geq 0}$ is Cauchy in $L^2$, and the completeness of $L^2(\Omega,\mathcal{F},\mathbb{P})$ gives a random variable $Y_t\in L^2$ such that
\begin{align*}
Y_t=\sum_{j=0}^{\infty}\psi_j Z_{t-j}
\end{align*}
with convergence in $L^2$.
The process $(Y_t)_{t\in\mathbb{Z}}$ is weakly stationary. Indeed, $\mathbb{E}[Y_t]=0$ because each $Z_{t-j}$ has mean $0$ and the partial sums converge in $L^2$, hence in $L^1$. For $h\in\mathbb{Z}$, let $Y_{t,N}=\sum_{j=0}^{N}\psi_jZ_{t-j}$ and $Y_{t+h,N}=\sum_{j=0}^{N}\psi_jZ_{t+h-j}$. Since $Y_{t,N}\to Y_t$ and $Y_{t+h,N}\to Y_{t+h}$ in $L^2$, continuity of the $L^2$ inner product gives
\begin{align*}
\mathbb{E}[Y_{t+h}\overline{Y_t}]
=
\lim_{N\to\infty}\mathbb{E}[Y_{t+h,N}\overline{Y_{t,N}}].
\end{align*}
For each $N$, expanding the finite sums and using the white-noise covariance relation $\mathbb{E}[Z_r\overline{Z_s}]=\sigma^2\mathbb{1}_{\{r=s\}}$ gives
\begin{align*}
\mathbb{E}[Y_{t+h,N}\overline{Y_{t,N}}]
&=
\sigma^2\sum_{j=0}^{N}\sum_{k=0}^{N}\psi_j\overline{\psi_k}\,\mathbb{1}_{\{t+h-j=t-k\}}.
\end{align*}
The finite sums converge as $N\to\infty$ to the absolutely convergent series
\begin{align*}
\mathbb{E}[Y_{t+h}\overline{Y_t}]
&=
\sigma^2\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}\psi_j\overline{\psi_k}\,\mathbb{1}_{\{t+h-j=t-k\}}\\
&=
\sigma^2\sum_{j=\max(0,h)}^{\infty}\psi_j\overline{\psi_{j-h}},
\end{align*}
where terms with negative indices are omitted. Absolute convergence follows from
\begin{align*}
\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}|\psi_j|\,|\psi_k|\,\mathbb{1}_{\{t+h-j=t-k\}}
\leq
\left(\sum_{j=0}^{\infty}|\psi_j|\right)^2
<\infty.
\end{align*}
The final expression depends only on $h$, not on $t$. Thus $(Y_t)$ is weakly stationary.
[/step]
[step:Show that the constructed process satisfies the ARMA equation]
Write
\begin{align*}
\varphi_0:=1,\qquad \varphi_k:=-\phi_k \text{ for } 1\leq k\leq p,
\end{align*}
so that
\begin{align*}
\phi(z)=\sum_{k=0}^{p}\varphi_k z^k.
\end{align*}
Also define $\vartheta_0:=1$, $\vartheta_k:=\theta_k$ for $1\leq k\leq q$, and $\vartheta_m:=0$ for $m>q$.
The identity $\phi(z)\psi(z)=\theta(z)$ implies equality of Taylor coefficients:
\begin{align*}
\sum_{k=0}^{\min(p,m)}\varphi_k\psi_{m-k}=\vartheta_m
\end{align*}
for every $m\geq 0$. Since $(\psi_j)$ is absolutely summable, for each fixed $k\in\{0,\dots,p\}$ the shifted series $\sum_{j=0}^{\infty}\psi_jZ_{t-k-j}$ converges in $L^2$. The outer sum over $k$ has only $p+1$ terms, so linearity and continuity of addition in $L^2$ permit interchange of this finite outer sum with the $L^2$ limits of the shifted series. Thus
\begin{align*}
\phi(B)Y_t
&=
\sum_{k=0}^{p}\varphi_k Y_{t-k}\\
&=
\sum_{k=0}^{p}\varphi_k\sum_{j=0}^{\infty}\psi_j Z_{t-k-j}\\
&=
\sum_{m=0}^{\infty}\left(\sum_{k=0}^{\min(p,m)}\varphi_k\psi_{m-k}\right)Z_{t-m}\\
&=
\sum_{m=0}^{q}\vartheta_m Z_{t-m}\\
&=
\theta(B)Z_t.
\end{align*}
Therefore the process $(Y_t)_{t\in\mathbb{Z}}$ is an $\operatorname{ARMA}(p,q)$ solution with the asserted infinite moving-average representation.
[guided]
We now check that the formal filter calculation is valid as an identity in $L^2$. Introduce the coefficients
\begin{align*}
\varphi_0:=1,\qquad \varphi_k:=-\phi_k \text{ for } 1\leq k\leq p,
\end{align*}
so that
\begin{align*}
\phi(z)=\sum_{k=0}^{p}\varphi_k z^k.
\end{align*}
Similarly define the moving-average coefficients by $\vartheta_0:=1$, $\vartheta_k:=\theta_k$ for $1\leq k\leq q$, and $\vartheta_m:=0$ for $m>q$.
The analytic identity
\begin{align*}
\phi(z)\psi(z)=\theta(z)
\end{align*}
holds for $|z|<\rho$. Equating Taylor coefficients gives, for every $m\geq 0$,
\begin{align*}
\sum_{k=0}^{\min(p,m)}\varphi_k\psi_{m-k}=\vartheta_m.
\end{align*}
This coefficient identity is exactly the filter identity needed to verify the ARMA equation.
Because $\phi$ has only finitely many coefficients and $\sum_j|\psi_j|<\infty$, each shifted moving-average series converges in $L^2$. The sum over $k\in\{0,\dots,p\}$ is finite, so we are only using linearity and continuity of finite addition in $L^2$, not an unconditional rearrangement of a double infinite series. Hence
\begin{align*}
\phi(B)Y_t
&=
\sum_{k=0}^{p}\varphi_k Y_{t-k}\\
&=
\sum_{k=0}^{p}\varphi_k\sum_{j=0}^{\infty}\psi_j Z_{t-k-j}.
\end{align*}
Now collect the coefficient of the same shock $Z_{t-m}$. The index relation is $m=k+j$, so $j=m-k$, and the admissible indices are $0\leq k\leq \min(p,m)$. Therefore
\begin{align*}
\phi(B)Y_t
&=
\sum_{m=0}^{\infty}\left(\sum_{k=0}^{\min(p,m)}\varphi_k\psi_{m-k}\right)Z_{t-m}\\
&=
\sum_{m=0}^{\infty}\vartheta_m Z_{t-m}\\
&=
\sum_{m=0}^{q}\vartheta_m Z_{t-m}\\
&=
\theta(B)Z_t.
\end{align*}
Thus the process built from the impulse response coefficients satisfies the original $\operatorname{ARMA}(p,q)$ equation.
[/guided]
[/step]
[step:Prove that any causal weakly stationary solution equals the constructed one]
Let $(X_t)_{t\in\mathbb{Z}}$ be any causal weakly stationary solution of
\begin{align*}
\phi(B)X_t=\theta(B)Z_t.
\end{align*}
The constructed process $(Y_t)_{t\in\mathbb{Z}}$ is causal because each $Y_t$ is the $L^2$ limit of finite linear combinations of $Z_t,Z_{t-1},\dots$. Define
\begin{align*}
W_t:=X_t-Y_t,\qquad t\in\mathbb{Z}.
\end{align*}
Since both $(X_t)$ and $(Y_t)$ solve the same ARMA equation, subtracting the two equations gives the homogeneous equation
\begin{align*}
\phi(B)W_t=0.
\end{align*}
Define
\begin{align*}
\pi: \{z\in\mathbb{C}: |z|<\rho\} &\to \mathbb{C}\\
z &\mapsto \frac{1}{\phi(z)}.
\end{align*}
Let $(\pi_j)_{j\geq 0}$ be its Taylor coefficients:
\begin{align*}
\pi(z)=\sum_{j=0}^{\infty}\pi_j z^j.
\end{align*}
Repeating the Cauchy integral formula estimate from the first step, with $\psi$ replaced by $\pi$, gives
\begin{align*}
\sum_{j=0}^{\infty}|\pi_j|<\infty.
\end{align*}
For $N\geq 0$, define the polynomial map
\begin{align*}
P_N: \mathbb{C} &\to \mathbb{C}\\
z &\mapsto \sum_{j=0}^{N}\pi_j z^j.
\end{align*}
Since $\pi(z)\phi(z)=1$, the coefficients of $P_N(z)\phi(z)$ equal $1$ at degree $0$, equal $0$ at degrees $1,\dots,N$, and may be nonzero only in degrees $N+1,\dots,N+p$. Thus
\begin{align*}
P_N(B)\phi(B)W_t
=
W_t+\sum_{m=N+1}^{N+p}a_{m,N}W_{t-m},
\end{align*}
where
\begin{align*}
a_{m,N}:=\sum_{k=0}^{p}\varphi_k\pi_{m-k},
\end{align*}
with the convention $\pi_\ell=0$ for $\ell<0$.
Since $\phi(B)W_t=0$, the left-hand side is $0$. Hence
\begin{align*}
W_t=-\sum_{m=N+1}^{N+p}a_{m,N}W_{t-m}.
\end{align*}
Define the finite uniform bound
\begin{align*}
C_W:=\sup_{s\in\mathbb{Z}}\left\|W_s\right\|_{L^2}.
\end{align*}
This number is finite because
\begin{align*}
\left\|W_s\right\|_{L^2}
\leq
\left\|X_s\right\|_{L^2}+\left\|Y_s\right\|_{L^2},
\end{align*}
and weak stationarity of $(X_t)$ and $(Y_t)$ makes the two terms on the right independent of $s$. The triangle inequality in $L^2$ gives
\begin{align*}
\left\|W_t\right\|_{L^2}
\leq
\sum_{m=N+1}^{N+p}|a_{m,N}|\,\left\|W_{t-m}\right\|_{L^2}
\leq
C_W\sum_{m=N+1}^{N+p}|a_{m,N}|.
\end{align*}
For $m\in\{N+1,\dots,N+p\}$ and $k\in\{0,\dots,p\}$, the index $m-k$ is at least $N+1-p$. Therefore
\begin{align*}
\sum_{m=N+1}^{N+p}|a_{m,N}|
&\leq
\sum_{m=N+1}^{N+p}\sum_{k=0}^{p}|\varphi_k|\,|\pi_{m-k}|\\
&\leq
\left(\sum_{k=0}^{p}|\varphi_k|\right)
\sum_{\ell=N+1-p}^{\infty}|\pi_\ell|,
\end{align*}
where terms with $\ell<0$ are harmless for fixed small $N$ and the limit as $N\to\infty$ is governed by the displayed tail. Since $(\pi_j)$ is absolutely summable, the final tail tends to $0$ as $N\to\infty$. The left-hand side $\|W_t\|_{L^2}$ is independent of $N$, so letting $N\to\infty$ gives $\|W_t\|_{L^2}=0$. Hence $W_t=0$ in $L^2$ for every $t\in\mathbb{Z}$. Hence $X_t=Y_t$ in $L^2$ for every $t$, and
\begin{align*}
X_t=\sum_{j=0}^{\infty}\psi_j Z_{t-j}.
\end{align*}
The representation is causal because the right-hand side uses only $Z_t,Z_{t-1},\dots$, and the coefficients are uniquely determined as the Taylor coefficients of $\theta/\phi$.
[guided]
We must show that no other causal weakly stationary solution can differ from the process already constructed. Let $(X_t)_{t\in\mathbb{Z}}$ be a causal weakly stationary solution of
\begin{align*}
\phi(B)X_t=\theta(B)Z_t.
\end{align*}
The process $(Y_t)$ constructed above is also causal: for each $t$, the partial sums $Y_{t,N}=\sum_{j=0}^{N}\psi_jZ_{t-j}$ are finite linear combinations of present and past noise variables, and $Y_{t,N}\to Y_t$ in $L^2$. Thus $Y_t$ lies in the closed $L^2$ linear span of $Z_t,Z_{t-1},\dots$.
Define the difference process by
\begin{align*}
W_t:=X_t-Y_t,\qquad t\in\mathbb{Z}.
\end{align*}
Since both $X$ and $Y$ satisfy the same ARMA equation, subtracting the two equations gives the homogeneous equation
\begin{align*}
\phi(B)W_t=0.
\end{align*}
We do not need to claim that $W$ is weakly stationary. The goal is instead to prove that this homogeneous equation, together with the uniform $L^2$ bounds inherited from the weak stationarity of $X$ and $Y$, forces $W$ to vanish.
Define the inverse autoregressive transfer function
\begin{align*}
\pi: \{z\in\mathbb{C}: |z|<\rho\} &\to \mathbb{C}\\
z &\mapsto \frac{1}{\phi(z)}.
\end{align*}
This map is analytic on the same disc because $\phi$ has no zero there. Let $(\pi_j)_{j\geq 0}$ denote its Taylor coefficients:
\begin{align*}
\pi(z)=\sum_{j=0}^{\infty}\pi_j z^j.
\end{align*}
Repeating the Cauchy integral formula estimate used for $\psi$ shows that $(\pi_j)$ is absolutely summable:
\begin{align*}
\sum_{j=0}^{\infty}|\pi_j|<\infty.
\end{align*}
This absolute summability is the stability input: it lets us approximate the inverse filter by finite filters whose error is a small tail.
For $N\geq 0$, define
\begin{align*}
P_N: \mathbb{C} &\to \mathbb{C}\\
z &\mapsto \sum_{j=0}^{N}\pi_j z^j.
\end{align*}
Since $\pi(z)\phi(z)=1$, the full infinite product has coefficient $1$ at degree $0$ and coefficient $0$ at every positive degree. Truncating $\pi$ at degree $N$ preserves those coefficient identities through degree $N$, and the only possible remaining terms in $P_N(z)\phi(z)$ have degrees $N+1,\dots,N+p$, because $\phi$ has degree $p$. Therefore
\begin{align*}
P_N(B)\phi(B)W_t
=
W_t+\sum_{m=N+1}^{N+p}a_{m,N}W_{t-m},
\end{align*}
where
\begin{align*}
a_{m,N}:=\sum_{k=0}^{p}\varphi_k\pi_{m-k},
\end{align*}
and we use the convention $\pi_\ell=0$ for $\ell<0$.
Now use the homogeneous equation. Since $\phi(B)W_t=0$, applying the finite filter $P_N(B)$ gives $P_N(B)\phi(B)W_t=0$. Hence
\begin{align*}
W_t=-\sum_{m=N+1}^{N+p}a_{m,N}W_{t-m}.
\end{align*}
Define
\begin{align*}
C_W:=\sup_{s\in\mathbb{Z}}\left\|W_s\right\|_{L^2}.
\end{align*}
This supremum is finite. Indeed, the triangle inequality gives
\begin{align*}
\left\|W_s\right\|_{L^2}
\leq
\left\|X_s\right\|_{L^2}+\left\|Y_s\right\|_{L^2},
\end{align*}
and the two terms on the right are independent of $s$ because $(X_t)$ and $(Y_t)$ are weakly stationary. Applying the triangle inequality in $L^2$ gives
\begin{align*}
\left\|W_t\right\|_{L^2}
\leq
\sum_{m=N+1}^{N+p}|a_{m,N}|\,\left\|W_{t-m}\right\|_{L^2}
\leq
C_W\sum_{m=N+1}^{N+p}|a_{m,N}|.
\end{align*}
It remains to check that the coefficient sum tends to $0$. For $m\in\{N+1,\dots,N+p\}$ and $k\in\{0,\dots,p\}$, the index $m-k$ satisfies $m-k\geq N+1-p$. Thus
\begin{align*}
\sum_{m=N+1}^{N+p}|a_{m,N}|
&\leq
\sum_{m=N+1}^{N+p}\sum_{k=0}^{p}|\varphi_k|\,|\pi_{m-k}|\\
&\leq
\left(\sum_{k=0}^{p}|\varphi_k|\right)
\sum_{\ell=N+1-p}^{\infty}|\pi_\ell|.
\end{align*}
The last expression tends to $0$ as $N\to\infty$ because $(\pi_j)$ is absolutely summable and $p$ is fixed. In the previous $L^2$ bound the left-hand side $\|W_t\|_{L^2}$ is independent of $N$, so letting $N\to\infty$ yields $\|W_t\|_{L^2}=0$. Hence $W_t=0$ in $L^2$ for every $t\in\mathbb{Z}$.
Therefore $X_t=Y_t$ in $L^2$ for every $t$, and the already constructed representation gives
\begin{align*}
X_t=\sum_{j=0}^{\infty}\psi_jZ_{t-j}.
\end{align*}
The right-hand side is causal because it only involves present and past shocks. The coefficients are uniquely determined by the Taylor expansion of $\psi(z)=\theta(z)/\phi(z)$, so the impulse response function is unique.
[/guided]
[/step]