It is shown how response functions can be derived as derivatives of a time-averaged quasienergy for molecules surrounded by a polarizable dielectric continuum. The vacuum quasienergy is modified with additional solvent terms which depend on the chosen solvent configuration. We consider equilibrium and nonequilibrium solvent models. Coupled cluster response functions are derived by introducing a corresponding coupled cluster quasienergy variational Lagrangian. Computational tractable expressions are derived for linear, quadratic, and cubic response functions. Sample calculations are reported for solvation effects on excitation energies, oscillator strengths, dipole moments, and frequency-dependent polarizabilities of H2S and furan.