The orientation time correlation functions (OTCFs) of molecular liquids are related to many experimental observations but subtle in theoretical calculations. Usually, the OTCFs at different timescales are described by distinct approaches. The OTCFs at short times were expanded into a series in an order of molecular angular velocity components, and the second cumulant (SC) of the series gave a single-frequency integral over one-mode propagation function weighted by a rotational density of states (DOS) of the liquid, whereas the OTCFs at long times displayed an exponential decay according to the Debye theory. In this paper, we derived the series of the OTCFs up to the fourth cumulant (FC), which gave a sum of twofold frequency integrals over two-mode propagation functions weighted by a product of two rotational DOSs associated with one or two molecular axes. With no need of any fitting parameter, the OTCFs calculated up to the SC and up to the FC were compared with simulation results for liquid water. With the power spectra of angular velocity autocorrelation function, which have a finite value at zero frequency, as the inputs, the OTCFs in both SC and FC schemes resulted in an exponential decay at long times, while the decay times were smaller than the simulation values. The positive results were the calculations with stable rotational instantaneous-normal-mode spectra, which have a linear behavior extrapolated to zero at low frequencies. As the stable rotational INM spectra were used as the inputs, the OTCFs in the SC scheme produced a slow power-law decay, but the OTCFs in the FC scheme captured the simulation behaviors from short to long timescale, including an exponential decay, with the decay times of the 2nd-rank OTCFs close to the simulation values.