For the chaotic time series single-step and multi-step prediction, Gaussian processes (GPs) method based on composite covariance function is proposed. GP priors over functions are determined mainly by covariance functions, and through learning from training data sets, all hyperparameters that define the covariance function and mean function can be estimated by using matrix operations and optimal algorithms within evidence maximum bayesian framework. As a probabilistic kernel machine, the number of tunable parameters for a GP model is greatly reduced compared with those for neural networks and fuzzy models. GP models with different composite covariance functions are applied to chaotic time series single-step and multi-step ahead prediction and compared with other models such as standard GP model with single covariance function, standard support vector machines, least square support vector machine, radial basis functional (RBF) neural networks, etc. Simulation results reveal that GP method with using different composite covariance functions can be used to accurately predict the chaotic time series and show stable performance with robustness. Hence,it provides an effective approach to studying the properties of complex nonlinear system modeling and control.