We present a theory of the four-wave-mixing (FWM) spectroscopy for the multilevel system coupled to the bath with arbitrary time scale, in contrast to most previous work that assumes the fast bath dynamic limit (motionally narrowed limit). Exact quantum mechanical expressions for the relevant nonlinear response tensor and the nonlinear polarization vector are obtained for a model Hamiltonian. Effects of molecular rotation and the polarization directions of the input radiation fields are also fully taken into account. The results show how the cross correlation as well as the autocorrelation functions of the fluctuations in the energy-gaps or transition frequencies between eigenstates contribute to the FWM spectroscopy. The result is applied to three-pulse echo spectroscopy. We examine the results for a variety of model spectral density functions, and recover the results of the previously reported theories as limiting cases.