In this paper, a new methodology for identifying multiple inputs multiple outputs Hammerstein systems is presented. The proposed method aims at incorporating the impulse response of the system into a least-squares support vector machine (LS-SVM) formulation and therefore the regularisation capabilities of LS-SVM are applied to the system as a whole. One of the main advantages of this method comes from the fact that it is flexible concerning the class of problems it can model and that no previous knowledge about the underlying non-linearities is required except for very mild assumptions. Also, it naturally adapts to handle different numbers of inputs/outputs and performs well in the presence of white Gaussian noise. Finally, the method incorporates information about the structure of the system but still the solution of the model follows from a linear system of equations. The performance of the proposed methodology is shown through three simulation examples and compared with other methods in the literature.