For a system strongly coupled to a heat bath, the quantum coherence of the system and the heat bath plays an important role in the system dynamics. This is particularly true in the case of non-Markovian noise. We rigorously investigate the influence of system-bath coherence by deriving the reduced hierarchal equations of motion (HEOM), not only in real time, but also in imaginary time, which represents an inverse temperature. It is shown that the HEOM in real time obtained when we include the system-bath coherence of the initial thermal equilibrium state possess the same form as those obtained from a factorized initial state. We find that the difference in behavior of systems treated in these two manners results from the difference in initial conditions of the HEOM elements, which are defined in path integral form. We also derive HEOM along the imaginary time path to obtain the thermal equilibrium state of a system strongly coupled to a non-Markovian bath. Then, we show that the steady state hierarchy elements calculated from the real-time HEOM can be expressed in terms of the hierarchy elements calculated from the imaginary-time HEOM. Moreover, we find that the imaginary-time HEOM allow us to evaluate a number of thermodynamic variables, including the free energy, entropy, internal energy, heat capacity, and susceptibility. The expectation values of the system energy and system-bath interaction energy in the thermal equilibrium state are also evaluated. (C) 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.