We present the exact joint likelihood of pseudo-$C_\ell$ power spectrum estimates measured from an arbitrary number of Gaussian cosmological fields. Our method is applicable to both spin-0 fields and spin-2 fields, including a mixture of the two, and is relevant to Cosmic Microwave Background, weak lensing and galaxy clustering analyses. We show that Gaussian cosmological fields are mixed by a mask in such a way that retains their Gaussianity, without making any assumptions about the mask geometry. We then show that each auto- or cross-pseudo-$C_\ell$ estimator can be written as a quadratic form, and apply the known joint distribution of quadratic forms to obtain the exact joint likelihood of a set of pseudo-$C_\ell$ estimates in the presence of an arbitrary mask. Considering the polarisation of the Cosmic Microwave Background as an example, we show using simulations that our likelihood recovers the full, exact multivariate distribution of $EE$, $BB$ and $EB$ pseudo-$C_\ell$ power spectra. Our method provides a route to robust cosmological constraints from future Cosmic Microwave Background and large-scale structure surveys in an era of ever-increasing statistical precision.