This paper presents a method for the identification of multiple-input-multiple-output (MIMO) Hammerstein systems for the goal of prediction. The method extends the numerical algorithms for subspace state space system identification (N4SID), mainly by rewriting the oblique projection in the N4SID algorithm as a set of componentwise least squares support vector machines (LS-SVMs) regression problems. The linear model and static nonlinearities follow from a low-rank approximation of a matrix obtained from this regression problem.