This paper focuses on time-domain identification issues of multi-input multi-output (MIMO) fractional order Hammerstein systems which are the extension of traditional Hammerstein type models by allowing linear part to be fractional order systems. The principal component analysis (PCA) method in subspace family is extended to identify coefficient matrixes of fractional order systems. Singular value decomposition (SVD) is utilized to estimate the unknown parameters of nonlinear part of system directly. A proper instrumental variable is chosen to eliminate the bias of identification results. Numerical simulation validates the proposed method.