This paper examines the problem of estimating linear time-invariant state-space system models. In particular, it addresses the parametrization and numerical robustness concerns that arise in the multivariable case. These difficulties are well recognised in the literature, resulting (for example) in extensive study of subspace-based techniques, as well as recent interest in 'data driven' local co-ordinate approaches to gradient search solutions. The paper here proposes a different strategy that employs the expectation-maximisation (EM) technique. The consequence is an algorithm that is iterative, with associated likelihood values that are locally convergent to stationary points of the (Gaussian) likelihood function. Furthermore, theoretical and empirical evidence presented here establishes additional attractive properties such as numerical robustness, avoidance of difficult parametrization choices, the ability to naturally and easily estimate non-zero initial conditions, and moderate computational cost. Moreover, since the methods here are maximum-likelihood based, they have associated known and asymptotically optimal statistical properties. (c) 2005 Elsevier Ltd. All rights reserved.