There is an important reason why the ARIMA might be preferred when the series are stationary. And this reason is the Wold's decomposition theorem - any covariance stationary process has a linear representation: a linear deterministic component ($V_t$) and a linear indeterministic components ($\varepsilon_t$)
Suppose that ${X_t}$ is a covariance stationary process with $\mathbb{E}[X_t] = 0$ and
covariance function, $\gamma(j) = \mathbb{E}[X_t X_{t−j}]$ , $ \forall j$. Then
$$X_t = \sum_{j=0}^{\infty} \psi_j \varepsilon_{t−j} + V_t$$
where
- $\psi_0=1$, $\sum_{j=0}^{\infty} \psi_j^2<\infty$
- $\varepsilon_{t−j} \sim WN(0, \sigma_{\varepsilon}^2)$
- $\mathbb{E}[\varepsilon_t V_s] = 0, \forall s,t>0$
- $\varepsilon_t = X_t - \mathbb{E}[X_t|X_{t-1},X_{t-2},...]$
As you may see, the first part of the representation looks like an $MA(\infty)$ process with square summable moving average terms. The second part is the deterministic part of $X_t$ because $V_t$ is perfectly predictable based on
past observations on $X_t$. And we know that models of $MA(\infty)$ representations are in their most general form $ARMA(p,q)$ representations: as long as the roots of the autoregressive part of an ARMA process are less than unity in absolute value, the process has a $MA(\infty)$ representation.
However, note, while an ARMA process generates an $MA(\infty)$ with square summable weights, it is not the only form that does this. A process that is square summable is not necessarily absolutely summable. $ARMA(p,q)$ models have
‘short memory’ relative to the entire class representations envisioned by the Wold representation. But Wold representation - despite covering more general cases- provides us with a strong argument of why modelling with ARMA is justifiable on stationary, short memory series.
Note also, another example of a stationary process is the periodic processes.If $Z_1,Z_2$ independent $N(0,\sigma^2)$ and $\omega$ constant then the process
$$X_t = Z_1 \cos (t \omega) + Z_2 \sin (t \omega)$$
is second order stationary with mean zero and autocovariance $cov(X_t,X_{t+h}) = \sigma^2 cos (\omega h) $.