Observations of exoplanets and protoplanetary disks show that binary stellar systems can host planets in stable orbits. Given the high binary fraction among stars, the contribution of binary systems to Galactic habitability should be quantified. Therefore, we have designed a suite of Monte Carlo experiments aimed at generating large (up to $10^6$) samples of binary systems. For each system randomly extracted we calculate the intersection between the radiative habitable zones and the regions of dynamical stability using published empirical formulations that account for the dynamical and radiative parameters of both stars of the system. We also consider constraints on planetary formation in binary systems. We find that the habitability properties of circumstellar and circumbinary regions are quite different and complementary with respect to the binary system parameters. Circumbinary HZs are, generally, rare ($simeq 4%$) in the global population of binary systems, even if they are common for stellar separations $lesssim 0.2$ AU. Conversely, circumstellar HZs are frequent ($ge 80%$) in the global population, but are rare for stellar separations $lesssim 1$ AU. These results are robust against variations of poorly constrained binary systems parameters. We derive ranges of stellar separations and stellar masses for which HZs in binary systems can be wider than the HZs around single stars; the widening can be particularly strong (up to one order of magnitude) for circumstellar regions around M-type secondary stars. The comparison of our statistical predictions with observational surveys shows the impact of selection effects on the habitability properties of detected exoplanets in binary systems.