Why is ergodicity not a requirement for ARIMA models besides stationarity? I frequently read that ARIMA models must be fitted on stationary data. But stationarity does not ensure ergodicity, which I understand is necessary to deduce population parameters from a single time series sample. Why is ergodicity not a requirement for ARIMA modeling? Do we just assume it?
Also, is there an example of a ergodic, but non-stationary process? Can you forecast these types of series?
 A: 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.
A: Ergodicity and mean-ergodicity are not the same properties.
Ergodicity is a much stronger property than mean-ergodicity (mean-ergodicity just means an $L^2$-LLN holds). There are easy examples of ARMA processes which are not ergodic.
What was shown by previous answer is that an ARMA process is mean-ergodic. (This is simply because that $l^1$, the space of absolutely summable sequences, is closed under convolution, and this makes the autocovariances also $l^1$, which implies mean-ergodicity.)

Why is ergodicity not a requirement for ARIMA modeling?

There is no reason for it to be. These notions have different historical origins. Ergodicity was first introduced in statistical mechanics, and intended to capture the phenomenon that "time average equals ensemble average". On the other hand, ARIMA models were introduced by Box and Jenkins for time series modeling.
You can already see from the definitions that they occur in different settings. Ergodicity is a property defined for strictly stationary processes, whereas ARMA processes are considered under covariance-stationarity.
From a time series perspective, first, the strict stationarity under which ergodicity is considered is too stringent an assumption to impose on general data. Second, the weak LLN that holds for many covariance-stationary processes (e.g. under $l^1$-condition for the autocovariances) is empirically just as good as the strong ergodic LLN.
For a good while, these two literatures developed separately and did not talk to each other. Later there were attempts to link the two notions by characterizing when ARMA processes satisfy strong-mixing type of conditions, which is a strengthening of ergodicity for more general processes (by, e.g. Kolmogorov and co-authors). But the connection is still incomplete.

...is there an example of a ergodic, but non-stationary process?

As stated above, ergodic processes are by definition strictly stationary.
