A least-squares method procedure for synthesizing the discrete time series that is characteristic of the uniform samples of the response of linear structural systems to stationary random excitation is described. The structural system is assumed to be an n-degree-of-freedom system that is representable by a set of ordinary differential equations excited by a vector white noise force. It is known that the discrete time series of uniformly spaced samples of a scalar white noise excited stationary linear differential equation can be represented as an autoregressive-moving average (AR-MA) time series and that the parameters of the AR-MA model can be computed from the covariance function of the differential equation model. The contributions of this paper are (i) the result that a scalar input scalar output AR-MA model duplicates the scalar output covariance function of a regularly sampled linear structural system with a multivariate white noise input, (ii) a computationally efficient method for computing the covariance function of a randomly excited structural system, and (iii) a demonstration of the theory and the numerical details of a two-stage least-squares procedure for the computation of the AR-MA parameters from the output covariance functions data.