We present density-functional theory for linear and nonlinear response functions using an explicit exponential parametrization of the density operator. The response functions are derived using two alternative variation principles, namely, the Ehrenfest principle and the quasienergy principle, giving different but numerically equivalent formulas. We present, for the first time, calculations of dynamical hyperpolarizabilities for hybrid functionals including exchange-correlation functionals at the general gradient-approximation level and fractional exact Hartree-Fock exchange. Sample calculations are presented of the first hyperpolarizability of the para-nitroaniline molecule and of a porphyrin derived push-pull molecule, showing good agreement with available experimental data.