Control systems are increasingly applied in domains where an analytic description of the system dynamics does not exist or is difficult to obtain. Example applications include autonomous robots in unstructured environments, human behavior modeling for prediction and action recognition in human-machineinteraction, and chemical process industry. In many of these cases, classical system identification is challenging, because a parametric model structure is unknown. Data-driven nonparametric models such as Gaussian process statespace models (GPSSMs) offer a suitable alternative: GPSSMs are known for their data-efficiency and rely on Bayesian principles to include prior knowledge.However, properties like stability or boundedness are often known a priori, but rarely exploited during modeling. We therefore propose a novel approach for learning GPSSMs subject to stability constraints. Our approach enforces the convergence using control Lyapunov functions which are also obtained in a datadriven fashion. We analyze the resulting dynamics with respect to convergence radius and data collection. In simulation, we illustrate the precision of the identified model on a real-world dataset of goal-directed motions. 1 .