An alternative approach to constructing linear spectra in the condensed phase is the cumulant method, where the linear response function is directly constructed from the fluctuations of the excitation energy sampled along ground-state molecular dynamics (MD) simulations in thermal equilibrium. Direct solute-solvent interactions, anharmonic effects, collective chromophore-environment motion and environmental polarization effects 10 are all accounted for in the MD sampling, making it a highly promising, albeit computationally expensive, approach in biological systems and pigment-protein complexes. Additionally, the computational cost of sampling energy gap fluctuations can be significantly reduced using machine learning (ML) techniques. However, in practical calculations, the cumulant expansion is generally truncated at second order, corresponding to mapping energy gap fluctuations onto a bath of linearly coupled harmonic oscillators (Brownian Oscillator Model or BOM). This truncation is only exact for systems where energy-gap fluctuations follow Gaussian statistics, and already leads to errors in the harmonic GBOM Hamiltonian including changes in PES curvature and Duschinsky mode mixing effects, that are captured exactly by the FC approach. The inclusion of higher order cumulants is vital in computing the correct lineshapes many systems, and we have recently demonstrated that a third order correction term can be constructed directly from MD, yielding improvements of spectral shapes in systems with small non-Gaussian contributions to the energy gap fluctuations. However, the approach has the significant shortcoming of being numerically unstable, yielding unphysical spectra with negative absorbances and divergent lineshapes for moderate to strong non-Gaussian fluctuations. These divergences can be traced to the discarding of higher order terms in the cumulant expansion and have been observed previously in 1D anharmonic model potentials. Since it is in general difficult to determine a priori whether the inclusion of a third order cumulant correction will improve the spectral lineshape or cause unphysical spectral contributions, the benefit of extending the cumulant method beyond second order in realistic condensed phase systems is questionable.