Why would initializing an ARMA process in a deterministic way, prevent it from being stationary? Let's assume that we have process $X_t$ that follows $\phi(B)X_t=\theta(B)Z_t$, where $Z_t$ is a white noise process, and $\deg(\phi)=p$, $\deg(\theta)=q$.
Assuming that $\phi$ has no unit roots, we can state that $X_t=\frac{\theta(B)}{\phi(B)}Z_t$ is the unique stationary solution to the above equation. 
However, if we don't restrict ourselves to stationarity, then we can find infinite solutions, depending on how we initialize the first $p$ values of $X_t$, in a deterministic way. This is precisely what I don't get, in an intuitive way ...
 A: There is an answer to a similar question here discussing the relationship between stationarity and the recursive equation for a model.  Some additional remarks are pertinent to your case.
Firstly, it is important to understand that the recursive equation defining the ARMA($p,q$) process is not sufficient to fully specify the process, even with specification of the distribution of the underlying noise series $\boldsymbol{Z} = \{ Z_t| t \in \mathbb{Z} \}$.  The recursive equations lock in the auto-correlation of the process, and specification of the distribution of the underlying noise series $\boldsymbol{Z}$ is then sufficient to give an asymptotic marginal distribution for the observable values.  However, the full joint distribution for the observable series $\boldsymbol{X} = \{ X_t| t \in \mathbb{Z} \}$ also depends on the distribution of the "starting values".  In some cases this is specifically defined, and in cases where it is not, it is usual to take this to mean that you are using the unique stationary distribution.
ARMA($p,q$) model with white noise: You have specified an underlying white noise series so you have distribution $Z_t \sim \text{IID N}(0, \sigma^2)$ for the underlying series.  Assuming that $\max| \phi_i | <1$ you have asymptotic stationarity with limiting distribution:
$$X_\infty \sim \text{N} \Bigg( 0, \sigma^2 \sum_{i=0}^\infty \psi_i^2 \Bigg) \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \psi(B) \equiv \frac{\phi(B)}{\theta(B)}.$$
To get an ARMA process with strict stationarity, you would set $X_1, \cdots , X_p$ to have a joint normal distribution with zero mean, the above variance, and auto-correlation that is consistent with your specified ARMA model.  This would give you strict stationarity, in the sense that the joint distribution of observable values does not change when shifted in time.
If you were to use a different joint distribution for your "starting values", this would lead to a non-stationary model.  It would still be asymptotically stationary with the above limiting distribution, but the condition of strict stationarity would not hold.
Intuitive explanation: The recursive formula for an ARMA model specifies the relationship that each observable value has with previous values.  There is an infinite class of possible models that are consistent with this recursive specification (each represented by a full joint distribution for which that recursive equation holds).  For a model with auto-regressive coefficients inside the unit circle, there is a unique strictly stationary model in this class, but there are also an infinite number of other models that are non-stationary, but still have an asymptotic limiting distribution.  The basic thing to remember for the intuitive explanation is that the ARMA equation, which is merely a recursive equation does not fully specify the model.
