A generalized Hubbard - Stratonovitch transformation relating an integral over random unitary matrices to an integral over Efetov's unitary -model manifold, is introduced. This transformation adapts the supersymmetry method to disordered and chaotic systems that are modelled not by a Hamiltonian but by their scattering matrix or time-evolution operator. In contrast to the standard method, no saddle-point approximation is made, and no massive modes have to be eliminated. This first paper on the subject applies the generalized Hubbard - Stratonovitch transformation to Dyson's circular unitary ensemble. It is shown how a supersymmetric variant of the Harish-Chandra - Itzykson - Zuber formula can be used to compute, in the large-N limit, the n-level correlation function for any n. Non-trivial applications to random network models, quantum chaotic maps, and lattice gauge theory, are expected.