To investigate the role of quantum effects in vibrational spectroscopies, we have carried out numerically exact calculations of linear and nonlinear response functions for an anharmonic potential system nonlinearly coupled to a harmonic Oscillator bath. Although one cannot carry out the quantum calculations of the response functions with full molecular dynamics (MD) simulations for a realistic system which cansists of many molecules, it is possible to grasp the essence of the quantum effects on the vibrational spectra by employing a model Hamiltonian that describes an intra- or intermolecular vibrational motion in a condensed phase. The present model fully includes vibrational relaxation, while the stochastic model often used to simulate infrared spectra does not. We have employed the reduced quantum hierarchy equations of motion approach in the Wigner space representation to deal with nonperturbative, non-Markovian, and nonsecular system-path interactions. Taking the classical limit of the hierarchy equations of motion, we have obtained the classical equations of motion that describe the classical dynamics under the same physical conditions as in the quantum case. By comparing the classical Old quantum mechanically calculated linear and multidimensional spectra, we found that the profiles of spectra for a fast modulation case were similar, but different for a slow modulation case. In both the classical and quantum cases, we identified the resonant oscillation peak in the spectra, but the quantum peak shifted to the red compared with the classical one if the potential is an harmonic. The prominent quantum effect is the 1-2 transition peak, which appears only in the quantum mechanically calculated spectra as a result of anharmonicity in the potential or nonlinearity of the system bath coupling. While the contribution of the 1-2 transition is negligible in the fast modulation case, it becomes important in the slow modulation case a long as the amplitude of the frequency fluctuation is small. Thus, we observed a distinct difference between the classical and quantum mechanically calculated multidimensional spectra in the slow modulation case where spectral diffusion plays a role. This fact indicates that one may not reproduce the experimentally obtained multidimensional spectrum for high-frequency vibrational modes based on classical molecular dynamics simulations if the modulation that arises from surrounding molecules is weak and slow. A practical way to overcome the difference between the classical and quantum simulations was discussed.