This study presents a multiple-input multiple-output (MIMO) approach for multi-step-ahead time series prediction with a Gaussian process regression (GPR) model. We assess the forecasting performance of the GPR model with respect to several neural network architectures. The MIMO setting allows modelling the cross-correlations between all regions simultaneously. We find that the radial basis function (RBF) network outperforms the GPR model, especially for long-term forecast horizons. As the memory of the models increases, the forecasting performance of the GPR improves, suggesting the convenience of designing a model selection criteria in order to estimate the optimal number of lags used for concatenation