A unified approach is presented in determining autoregressive moving average (ARMA) algorithms for simulating realizations of multivariate random processes with a specified (target) spectral matrix. The ARMA algorithms are derived by relying on a prior autoregressive (AR) approximation of the target matrix. Several AR to ARMA procedures are formulated by minimizing a frequency domain error. Equations which can lead to a convenient computation of the ARMA matrix coefficients for a particular problem are given. Finally, the features of the various procedures are critically assessed.