We present density-functional theory for time-dependent response functions up to and including cubic response. The working expressions are derived from an explicit exponential parametrization of the density operator and the Ehrenfest principle, alternatively, the quasienergy ansatz. While the theory retains the adiabatic approximation, implying that the time-dependency of the functional is obtained only implicitly-through the time dependence of the density itself rather than through the form of the exchange-correlation functionals-it generalizes previous time-dependent implementations in that arbitrary functionals can be chosen for the perturbed densities (energy derivatives or response functions). In particular, general density functionals beyond the local density approximation can be applied, such as hybrid functionals with exchange correlation at the generalized-gradient approximation level and fractional exact Hartree-Fock exchange. With our implementation the response of the density can always be obtained using the stated density functional, or optionally different functionals can be applied for the unperturbed and perturbed densities, even different functionals for different response order. As illustration we explore the use of various combinations of functionals for applications of nonlinear optical hyperpolarizabilities of a few centrosymmetric systems; molecular nitrogen, benzene, and the C(60) fullerene. Considering that vibrational, solvent, and local field factors effects are left out, we find in general that very good experimental agreement can be obtained for the second dynamic hyperpolarizability of these systems. It is shown that a treatment of the response of the density beyond the local density approximation gives a significant effect. The use of different functional combinations are motivated and discussed, and it is concluded that the choice of higher order kernels can be of similar importance as the choice of the potential which governs the Kohn-Sham orbitals.