We consider large scale covariance estimation using a small number of samples in applications where there is a natural ordering between the random variables. The two classical approaches to this problem rely on banded covariance and banded inverse covariance structures, corresponding to time varying moving average (MA) and autoregressive (AR) models, respectively. Motivated by this analogy to spectral estimation and the well known modeling power of autoregressive moving average (ARMA) processes, we propose a novel time varying ARMA covariance structure. Similarly to known results in the context of AR and MA, we address the completion of an ARMA covariance matrix from its main band, and its estimation based on random samples. Finally, we examine the advantages of our proposed methods using numerical experiments.