We propose a novel method of time series decomposition based on the non-negative factorization of the Hankel matrix of time series and apply this method for time series modelling and prediction. An interim (surrogate) model of time series is built from the components of the time series using random cointegration, while the best cointegration is selected using a nature-inspired optimization method (Artificial Bee Colony). For modelling of cointegrated time series we use the ARX (AutoRegressive with eXogenous inputs) model. The results of modelling using the historical data (daily highest price) of S&P 500 stocks from 2009 are presented and compared against stand-alone ARX models. The results are evaluated using a variety of metrics (RMSE, MAE, MAPE, Pearson correlation, Nash-Suttcliffe efficiency coefficient, etc.) as well as illustrated graphically using Taylor and Target diagrams. The results show a 51–98% improvement of prediction accuracy (depending upon accuracy metric used). The proposed time series modelling method can be used for variety applications (time series denoising, prediction, etc.).