Employing multiple pulsars and using an appropriate algorithm to establish ensemble pulsar timescale can reduce the influences of various noises on the long-term stability of pulsar timescale, compared to a single pulsar. However, due to the low timing precision and significant red noises of some pulsars, their participation in the construction of ensemble pulsar timescale is often limited. Inspired by the principle of solving non-stationary sequence modeling using co-integration theory, we put forward an algorithm based on co-integration theory to establish an ensemble pulsar timescale. It is found that this algorithm can effectively suppress some noise sources if a co-integration relationship between different pulsar data exists. Different from the classical weighted average algorithm, the co-integration method provides the chance for a pulsar with significant red noises to be included in the establishment of an ensemble pulsar timescale. Based on data from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav), we found that the co-integration algorithm can successfully reduce several timing noises and improve the long-term stability of the ensemble pulsar timescale.