  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 the 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 puts forward an algorithm based on the
co-integration theory to establish ensemble pulsar timescale. It is found that
this algorithm can effectively suppress some noise sources if a co-integration
relationship between different pulsar data exist. Different from the classical
weighted average algorithm, the co-integration method provides the chances of
the pulsar with significant red noises to attend the establishment of ensemble
pulsar timescale. Based on the data from the North American Nanohertz
Observatory for Gravitational Waves, we found that the co-integration algorithm
can successfully reduce several timing noises and improve the long-term
stability of the ensemble pulsar timescale.
