We suggest a new method for building data-driven dynamical models from observed multidimensional time series. The method is based on a recurrent neural network with specific structure, which allows for the joint reconstruction of both a low-dimensional embedding for dynamical components in the data and an evolution operator. The key link of the method is a Bayesian optimization of both model structure and the hypothesis about the data generating law, which is needed for constructing the cost function for model learning. First, the performance of the method is successfully tested in the situation when a signal from a low-dimensional dynamical system is hidden in noisy multidimensional observations. Second, the method is used for building the data-driven model of the low frequency variability (LFV) in the quasigeostrophic model of the Earth's midlatitude atmosphere-a high-dimensional chaotic system. It is demonstrated that the key regimes of the atmospheric LFV are reproduced correctly in data simulations by means of the obtained model.