A traditional stochastic hierarchy equations of motion method is extended into the correlated real-time and imaginary-time propagations, in this paper, for its applications in calculating the equilibrium correlation functions. The central idea is based on a combined employment of stochastic unravelling and hierarchical techniques for the temperature-dependent and temperature-free parts of the influence functional, respectively, in the path integral formalism of the open quantum systems coupled to a harmonic bath. The feasibility and validity of the proposed method are justified in the emission spectra of homodimer compared to those obtained through the deterministic hierarchy equations of motion. Besides, it is interesting to find that the complex noises generated from a small portion of real-time and imaginary-time cross terms can be safely dropped to produce the stable and accurate position and flux correlation functions in a broad parameter regime.