The vibrational contributions to the average polarizability (ᾱ), to the second harmonic scattering (SHS) first hyperpolarizability (β SHS), and depolarization ratio (DR SHS), as well as to the third harmonic scattering (THS) second hyperpolarizability (γ THS) and depolarization ratio (DR THS), have been evaluated for the water molecule using the Bishop and Kirtman perturbative theory approach, in combination with finite differentiation techniques to evaluate the higher-order derivatives. From a hierarchy of coupled cluster techniques and extended atomic basis sets, the CCSD/d-aug-cc-pVTZ level has been selected to assess the importance of the zero-point vibrational average (ZPVA) contributions and of the pure vibrational contributions with respect to their electronic counterparts. This is the first investigation demonstrating electronic and vibrational SHS, and THS responses can be computed for small molecules, with the perspective of performing comparisons with recent experimental data [Van Steerteghem et al., Anal. Chem. 89, 2964 (2017) and V. Rodriguez, J. Phys. Chem. C 121, 8510 (2017)]. Numerical results on the water molecule highlight that (i) the vibrational contributions to the dynamic ᾱ, β SHS, and γ THS are small but non negligible; (ii) they amount to 3%, 10%, and 4% at the typical 1064 nm wavelength, respectively; (iii) the mechanical anharmonicity term dominates the ZPVA contribution; (iv) the double harmonic terms dominate the pure vibrational contributions; (v) the stretching vibrations provide the largest contributions to the dynamic (hyper)polarizabilities; and (vi) these conclusions are strongly impacted in the static limit where the vibrational contributions are much larger, in particular the double harmonic pure vibrational terms, and even more in the case of the first hyperpolarizability.