We apply random matrix theory to derive spectral density of large sample covariance matrices generated by multivariate VMA(q), VAR(q) and VARMA(q1,q2) processes. In particular, we consider a limit where the number of random variables N and the number of consecutive time measurements T are large but the ratio N/T is fixed. In this regime the underlying random matrices are asymptotically equivalent to Free Random Variables (FRV). We apply the FRV calculus to calculate the eigenvalue density of the sample covariance for several VARMA-type processes. We explicitly solve the VARMA(1,1) case and demonstrate a perfect agreement between the analytical result and the spectra obtained by Monte Carlo simulations. The proposed method is purely algebraic and can be easily generalized to q1>1 and q2>1.