We present a method to produce mock galaxy catalogues with efficient perturbation theory schemes, which match the number density, power spectra and bispectra in real and in redshift space from N-body simulations. The essential contribution of this work is the way in which we constrain the bias parameters of the PATCHY-code. In addition to aiming at reproducing the two-point statistics, we seek the set of bias parameters, which constrain the univariate halo probability distribution function (PDF) encoding higher order correlation functions. We demonstrate that halo catalogues based on the same underlying dark matter field with a fix halo number density, and accurately matching the power spectrum (within 2 per cent) can lead to very different bispectra depending on the adopted halo bias model. A model ignoring the shape of the halo PDF can lead to deviations up to factors of 2. The catalogues obtained additionally constraining the shape of the halo PDF can significantly lower the discrepancy in the three-point statistics, yielding closely unbiased bispectra both in real and in redshift space; which are in general compatible with those corresponding to an N-body simulation within 10 per cent (deviating at most up to 20 per cent). Our calculations show that the constant linear bias of ˜2 for luminous red galaxy (LRG) like galaxies found in the power spectrum, mainly comes from sampling haloes in high-density peaks, choosing a high-density threshold rather than from a factor multiplying the dark matter density field. Our method contributes towards an efficient modelling of the halo/galaxy distribution required to estimate uncertainties in the clustering measurements from galaxy redshift surveys. We have also demonstrated that it represents a powerful tool to test various bias models.