Raman intensities in molecular spectra are usually computed within double harmonic approximation. This procedure relies on treating a vibrating molecule as a collection of harmonic oscillators and on the assumption that polarisability tensor invariants display linear variations around the molecular equilibrium geometry. This methodology, originally formulated by Placzek, constitutes the theoretical foundation for computing Raman intensities in standard quantum chemistry packages. However, the two assumptions underlying double harmonic approximation have not been sufficiently tested. In this work, we employed exact anharmonic ro-vibrational wave functions and distance-dependent polarisability invariants together with their harmonic approximants to investigate the discrepancies in Raman intensities of the fundamental transitions in 12 diatomic molecules, caused by double harmonic approximation. We found that: (i) the errors in total Raman intensities were between −8.2% and +9.5%, (ii) the largest discrepancy was observed for F2, where the polarisability invariants could not be adequately modelled by their linear approximants, and (iii) quantum chemical methods fail to predict reliable polarisability invariants at non-equilibrium molecular geometries; the associated errors in Raman intensities are huge and completely overshadow the shortcomings of double harmonic approximation. We communicate here an urgent need for developing accurate methods capable of computing reliable polarisabilities also at distorted geometries.