Recent Hubble Space Telescope observations have unveiled the white dwarf cooling sequence of the Galactic bulge. Although the degenerate sequence can be well fitted employing the most up-to-date theoretical cooling sequences, observations show a systematic excess of red objects that cannot be explained by the theoretical models of single carbon-oxygen white dwarfs of the appropriate masses. Here we present a population synthesis study of the white dwarf cooling sequence of the Galactic bulge that takes into account the populations of both single white dwarfs and binary systems containing at least one white dwarf. These calculations incorporate state-of-the-art cooling sequences for white dwarfs with hydrogen-rich and hydrogen-deficient atmospheres, for both white dwarfs with carbon-oxygen and helium cores, and also take into account detailed prescriptions of the evolutionary history of binary systems. Our Monte Carlo simulator also incorporates all the known observational biases. This allows us to model with a high degree of realism the white dwarf population of the Galactic bulge. We find that the observed excess of red stars can be partially attributed to white dwarf plus main sequence binaries, and to cataclysmic variables or dwarf novae. Our best fit is obtained with a higher binary fraction and an initial mass function slope steeper than standard values, as well as with the inclusion of differential reddening and blending. Our results also show that the possible contribution of double degenerate systems or young and thick-disk bulge stars is negligible.