In this paper we develop a method for learning nonlinear system models with multiple outputs and inputs. We begin by modeling the errors of a nominal predictor of the system using a latent variable framework. Then using the maximum likelihood principle we derive a criterion for learning the model. The resulting optimization problem is tackled using a majorization-minimization approach. Finally, we develop a convex majorization technique and show that it enables a recursive identification method. The method learns parsimonious predictive models and is tested on both synthetic and real nonlinear systems. (C) 2018 Elsevier Ltd. All rights reserved.