A fixed set of output cumulants of order greater than two guarantees unique identification of known-order causal ARMA models, which are driven by unobservable non-Gaussian i.i.d. noise. The models are allowed to be nonminimum phase, and their outputs may be corrupted by additive colored Gaussian noise of unknown covariance. The ARMA parameters can be estimated either via linear equations and closed-form expressions, or by minimizing quadratic cumulant matching criteria. The latter approach requires computation of cumulants in terms of the ARMA parameters which is carried out in the state space using Kronecker products.