A new method is presented for the numerical identification of a pure random stationary process, a realization of which is given in the form of a time series. The sample covariance function (in the form of serial correlation coefficients) serves as input to an estimation programme which fits it to the autocovariance function of a mixed autoregressive, moving-average (AKMA) model. The fitting is done by maximum likelihood (ML) estimation of the model parameters. An asymptotic theory has been developed, valid for a somewhat more general model. Solution of the ML equation is accomplished by a multivariable Newton—Raphson procedure preceded by a strategic computer search. The estimates, which exist with a probability approaching unity, have the properties of uniqueness, consistency, efficiency, and a joint normal multivariate distribution. The ML method is shown to be asymptotically equivalent to a weighted least squares procedure, accomplished by the minimization of a quadratic form having chi-square distribution, a property which is used for hypothesis testing. In the case of rejection, ARMA models of different order must be tried until a hypothesis is accepted. The computer algorithms have been tried on simulated time series of ARMA models having up to seven parameters. The numerical consequences of various aspects of the asymptotic theory have been investigated and found to corroborate the theory.