A general computational protocol for accurate predictions of nonlinear optical (NLO) properties of solvated molecules based on the rigorous local field (RLF) approach taking all relevant effects into account is presented. para-Nitroaniline (pNA) was taken as a model NLO system dissolved in cyclohexane, tetrahydrofuran, and 1,4-dioxane. Molecular dynamics (MD) simulations employing either non-polarizable or polarizable force fields were used to generate representative sets of structures of the solutions. The static NLO properties of the solute were calculated at the MP2/aug-cc-pVDZ level of theory with the multiplicative scaling method used to account for the frequency dispersion of the properties. Focusing on the electric field-induced second harmonic generation (EFISH) and hyper-Rayleigh scattering (HRS), a good agreement between calculated results and experimental measurements was achieved with a polarizable force field. While the solvent effects on the vibrational contributions to the static molecular properties are significant, they remain small for both EFISH and HRS. Our results show that the proposed approach yields reliable predictions of dynamic NLO properties of solvated chromophores, which paves a route to further applications of the RLF approach to study a wide range of NLO phenomena in heterogeneous environments.