We study two linear estimators for stationary invertible VARMA models in echelon form – to achieve identification (model parameter unicity) – with known Kronecker indices. Such linear estimators are much simpler to compute than Gaussian maximum-likelihood estimators often proposed for such models, which require highly nonlinear optimization. The first estimator is an improved two-step estimator which can be interpreted as a generalized-least-squares extension of the two-step least-squares estimator studied in Dufour and Jouini (2005). The setup considered is also more general and allows for the presence of drift parameters. The second estimator is a new relatively simple three-step linear estimator which is asymptotically equivalent to ML, hence asymptotically efficient, when the innovations of the process are Gaussian. The latter is based on using modified approximate residuals which better take into account the truncation error associated with the approximate long autoregression used in the first step of the method. We show that both estimators are consistent and asymptotically normal under the assumption that the innovations are a strong white noise, possibly non-Gaussian. Explicit formulae for the asymptotic covariance matrices are provided. The proposed estimators are computationally simpler than earlier “efficient” estimators, and the distributional theory we supply does not rely on a Gaussian assumption, in contrast with Gaussian maximum likelihood or the estimators considered by Hannan and Kavalieris (1984b) and Reinsel, Basu and Yap (1992). We present simulation evidence which indicates that the proposed three-step estimator typically performs better in finite samples than the alternative multi-step linear estimators suggested by Hannan and Kavalieris (1984b), Reinsel et al. (1992), and Poskitt and Salau (1995).