In this paper, an identification method based on the recursive auxiliary variable least squares algorithm is proposed for a multi-input-multi-output Hammerstein-Wiener system with process noise. In the proposed identification method, the system is converted into the multivariate regression form under the condition that the nonlinear block in the output part is invertible. Then, the auxiliary variable is constructed, the parameters of the regression equations are identified, and the system parameter matrices can be obtained by matrix composition of the parameter product matrix. A theoretical analysis showed that the proposed method has uniform convergence when the process noise is white and has a finite variance. The effectiveness of the proposed method is validated through the experiments.