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?

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
}

