How do I restrict coefficient values to ensure stationarity and invertibility of ARMA(p,q)? I am trying to simulate an ARMA series but I am concerned about it being conformed to ensure invertibility of $\text{ARMA}(p,q)$ in addition to stationarity. I know of the following conditions for MA series:

*

*For an $\text{MA}(1)$ model:  $−1 < \theta_1 < 1$.

*For an $\text{MA}(2)$ model:  $−1 < \theta_2 < 1$, $\theta_2 + \theta_1 > −1$, $ \theta_1 − \theta_2 < 1$.

Moving average models in Forecasting Principles and Practice by Rob J Hyndman and George Athanasopoulos
I also know of the following conditions for AR series:


*For an $\text{AR}(1)$ model:  $−1 < \phi_1 < 1$.

*For an $\text{AR}(2)$ model:  $−1 < \phi_2 < 1$, $\phi_2 + \phi_1 > −1$, $ \phi_1 − \phi_2 < 1$.

Autoregressive models in Forecasting Principles and Practice by Rob J Hyndman and George Athanasopoulos
What conditions holds for $\text{ARMA}(1,1)$ and more?
 A: This is worked out in the book 'Time Series Analysis and Its Applications With R Examples, Third Edition' by Robert Shumway and David Stoffer.
Following their notation, let an $ARMA(p, q)$ model be expressed as $\phi(B)x_t = \theta(B)w_t$ (see equation (3.21) in page 93 of the book), where $\phi(z) = 1 - \phi_1 z − \cdots − \phi_p z^p$ with $\phi_p \neq 0$, $\theta(z) = 1 + \theta_1 z + \cdots + \theta_q z^q$ with $\theta_q \neq 0$, $z$ is a complex number and $B$ is the backshift operator (see Definition 3.6, page 94). Then, an $ARMA(p, q)$ model is invertible if and only if $\theta(z) \neq 0$ for $|z| \le 1$ (see Property 3.2 in page 95). The way of proving this statement is also discussed there.
Also check out these lectures notes. The stationarity of $ARMA(p, q)$ is ensured when all the (real and complex) roots $z_j$ of the polynomial $\phi(z)$ satisfy $|z_j| > 1$ (see page 11 of the lecture notes), i.e., $\phi(z) \neq 0$ for $|z| \le 1$.
In the particular case of $ARMA(1,1)$, $\phi(z) = 1 - \phi_1 z$ and $\theta(z) = 1 + \theta_1 z$. So, it is invertible and stationary if $0 < |\phi_1|, |\theta_1| < 1$.
A: I found a very useful material that answers this question the way I expected. Lecture 2: ARMA(p,q) models(part 3)
On page 6 precisely, it goes as follows:

*

*The (stability) stationarity condition is the one of an AR(1) process (or
ARMA(1,0) process) :
$|\phi| < 1$.

*The invertibility condition is the one of a MA(1) process (or
ARMA(0,1) process) :
$|\theta| < 1$.

*The representation of an ARMA(1,1) process is fundamental or causal
if :
$|\phi| < 1$ and $|\theta| < 1$.

*The representation of an ARMA(1,1) process is said to be minimal and
causal if :
$|\phi| < 1$, $|\theta| < 1$ and $\phi \neq \theta$.

A: In general you can reparameterize an ARMA model in terms of the partial autocorrelations associated with the AR part of the model and an analogous reparameterisation of the MA part, see Barndorff-Nielsen & Schou (1973).  These new parameters only need to satisfy simple box constraints and there is a one-to-one correspondance between the two parameterizations.
Thus, to check that the AR coefficients give you a stationary model, all you need to do is to transform the coeffients to the corresponding partial autocorrelations and check if these are smaller than 1 in absolute values.  Recursive formulas for this transformation as well as its inverse are concisely given in Monahan 1984.
The MA coefficients can be checked for invertibility in the same way.
Checking for stationarity and invertibility via this transformation is in practice considerably more efficient than computing the roots of the AR and MA polynomials numerically as suggested elsewhere.  It applies to ARMA models of any order.  This reparameterization also has other practical applications, see for example Sørby and Rue (2017).
A simple R implementation of this transformation and its inverse following the notation used by Monahan is given below.  As noted by Yves, these functions are also available in the FitAR R-package as functions ARToPacf and PacfToAR.
# Reparameterisation from AR or MA coeffients to "partial autocorrelations"
y_to_r <- function(y) {
  r <- numeric(length(y))
  if (length(y)==1) return(y[1])
  for (k in length(y):1) {
    r[k] <- y[k]
    y <- (y[1:(k-1)] + y[k]*y[(k-1):1])/(1 - y[k]^2)
  }
  r
}
# The inverse transformation
r_to_y <- function(r) {
  y <- r[1]
  if (length(r)==1) return(y)
  for (k in 2:length(r)) {
    y <- c(y - r[k]*rev(y), r[k])
  }
  y
}

