A bit technical maybe, but stationary ARMA processes are by construction mean-ergodic (as the other answer correctly pointed out, a previous version of my answer did not spell that out clearly and wrote ergodic as mean-ergodicity is maybe the most important "flavor" of ergodicity and hence sometimes treated synonymously with erdogicity, which, as this discussion shows, it should indeed not).
First, here is a sufficient condition for mean ergodicity:
Theorem:
Let $Y_t$ be covariance stationary with $E(Y_t)=\mu$ and $Cov(Y_t,Y_{t-j})=\gamma_j$ such that $\sum_{j=0}^\infty|\gamma_j|<\infty$. Then
$$\bar{Y}_T\to_p \mu$$
Proof:
We shall actually prove that $\bar{Y}_T$ converges to $\mu$ in mean square, which implies convergence in probability. Write
\begin{eqnarray*}
E(\bar{Y}_T- \mu)^2&=&E\left[(1/T)\sum_{t=1}^T(Y_t- \mu)\right]^2\\
&=&1/T^2E[\{(Y_1- \mu)+(Y_2- \mu)+\ldots+(Y_T- \mu)\}\\
&&\quad\{(Y_1- \mu)+(Y_2- \mu)+\ldots+(Y_T- \mu)\}]\\
&=&1/T^2\{[\gamma_0+\gamma_1+\ldots+\gamma_{T-1}]+[\gamma_1+\gamma_0+\gamma_1+\ldots+\gamma_{T-2}]\\
&&\quad+\ldots+[\gamma_{T-1}+\gamma_{T-2}+\ldots+\gamma_1+\gamma_0]\}
\end{eqnarray*}
Thus,
\begin{eqnarray*}
E(\bar{Y}_T- \mu)^2&=& 1/T^2\{T\gamma_0+2(T-1)\gamma_1+2(T-2)\gamma_2+\ldots+2\gamma_{T-1}\}
\end{eqnarray*}
Put differently,
\begin{eqnarray*}
E(\bar{Y}_T- \mu)^2&=& 1/T\{\gamma_0+2(T-1)\gamma_1/T+2(T-2)\gamma_2/T+\ldots+2\gamma_{T-1}/T\}
\end{eqnarray*}
This expression tends to zero as $T\to\infty$, as $TE(\bar{Y}_T- \mu)^2$ remains bounded, because
\begin{eqnarray*}
TE(\bar{Y}_T- \mu)^2&=& |\gamma_0+2(T-1)\gamma_1/T+2(T-2)\gamma_2/T+\ldots+2\gamma_{T-1}/T|\\
&\leqslant&|\gamma_0|+2(T-1)|\gamma_1|/T+2(T-2)|\gamma_2|/T+\ldots+2|\gamma_{T-1}|/T\\
&\leqslant&|\gamma_0|+2|\gamma_1|+2|\gamma_2|+\ldots+2|\gamma_{T-1}|\\
&\to&c<\infty,
\end{eqnarray*}
using summability of the autocovariances.
That is, if the autocovariances decay sufficiently quickly, mean ergodicity follows.
We next show that any causal $ARMA(p,q)$ process is ergodic, as it has the required summable autocovariances.
Let us look at the $MA(\infty)$ representation and use the triangle inequality to bound the sufficient condition for mean ergodicity of a stationary/causal process from above.
Stationarity implies that a causal, or $MA(\infty)$ with summable coefficients, representation of the process exists.
The claim is therefore shown if we can show that summability of the $MA(\infty)$ coefficients $\sum_{j=0}^\infty|\psi_j|<\infty$ implies $\sum_{k=0}^\infty|\gamma_k|<\infty$ where $\gamma_k=\sigma^2\sum_{j=0}^{\infty}\psi_j\psi_{j+k}$ is the $k$th autocovariance of an $MA(\infty)$-process.
We write
\begin{eqnarray*}
\sum_{k=0}^\infty|\gamma_k|&=&\sum_{k=0}^\infty\left|\sigma^2\sum_{j=0}^{\infty}\psi_j\psi_{j+k}\right|\\
&=&\sigma^2\sum_{k=0}^\infty\left|\sum_{j=0}^{\infty}\psi_j\psi_{j+k}\right|\\
&\leqslant&\sigma^2\sum_{k=0}^\infty\sum_{j=0}^{\infty}\left|\psi_j\psi_{j+k}\right|\\
&=&\sigma^2\sum_{k=0}^\infty\sum_{j=0}^{\infty}\left|\psi_j\right|\left|\psi_{j+k}\right|\\
&=&\sigma^2\sum_{j=0}^{\infty}\left|\psi_j\right|\sum_{k=0}^\infty\left|\psi_{j+k}\right|\\
&\leqslant&\sigma^2\sum_{j=0}^{\infty}\left|\psi_j\right|\sum_{k=0}^\infty\left|\psi_{k}\right|\\
&<&\infty
\end{eqnarray*}
Here, the first inequality uses the triangle inequality. Summability of the coefficients permits interchanging the order of summation in fourth equality (and hence taking out $|\psi_j|$ which does not depend on $k$). The second inequality follows because the second summation additionally has the terms $\psi_0,\ldots,\psi_{j-1}$ for $j>0$. The last inequality then follows from summability of the coefficients.