Methods based on path integral molecular dynamics (PIMD) are a family of chemical dynamics techniques which are a common way of treating nuclear quantum effects in condensed phase. One of the approximations shared by all of these is the neglect of quantum coherence. This is usually a reasonable assumption, because the large number of interactions present in a solid or a liquid typically leads to quick decoherence. To better understand the role that this assumption may have, we decided to explore a method for exact quantum dynamics in a dissipative environment. Hierarchical equations of motion (HEOM) are a computational method that describes time evolution of a system quantum density matrix coupled to a bath of harmonic oscillators, which in the appropriate limits yields the exact results. In this work, quantum HEOM were implemented both in their high-temperature form, omitting Matsubara frequency terms present in the bath time autocorrelation function, as well as the full form, which results in a many-dimensional set of coupled differential equations. The behaviour of HEOM was studied at various limits, including the low-temperature breakdown, and the importance of the Matsubara terms has been demonstrated. Carrying out a Wigner transform of HEOM, one can obtain a position-momentum phase space analogue, which at the limit of high temperature gives completely classical HEOM. These were implemented and compared to the quantum version. In addition, the computational program was tested by reproducing the results of Sakurai and Tanimura on the Morse oscillator and preliminary work has been done on the OH potential of a hydrated oxonium ion. Drawing from this work, we propose future studies of the blue shift, which currently plagues PIMD-based dynamics methods in strongly anharmonic potentials, and of the role of an explicit bath and quantum coherence in condensed phase dynamics.