We present a formal analysis of nonlinear response functions in terms of correlation functions in real-and imaginary-time domains. In particular, we show that causal nonlinear response functions, expressed in terms of nested commutators in real time, can be obtained from the analytic continuation of time ordered response functions, which are more easily amenable to diagrammatic calculation. This generalizes the well-known result of linear response theory. We then use gauge invariance arguments to derive exact relations between second-order response functions in density and current channels. These identities, which are non-perturbative in the strength of inter-particle interactions, allow us to establish exact connections between nonlinear optics calculations done in different electromagnetic gauges.