Supermassive black hole (BH) binaries would comprise the strongest sources of gravitational waves (GW) once they reach <<1 pc separations, for both pulsar timing arrays (PTAs) and space based (SB) detectors. While BH binaries coalescences constitute a natural outcome of the cosmological standard model and galaxy mergers, their dynamical evolution is still poorly understood and therefore their abundances at different stages. We use a dynamical model for the decay of BH binaries coupled with a cosmological simulation and semi-empirical approaches to the occupation of haloes by galaxies and BHs, in order to follow the evolution of the properties distribution of galaxies hosting BH binaries candidates to decay due to GWs emission. Our models allow us to relax simplifying hypothesis about the binaries occupation in galaxies and their mass, as well as redshift evolution. Following previously proposed electromagnetic (EM) signatures of binaries in the subpc regime, that include spectral features and variability, we model possible distributions of such signatures and also set upper limits to their lifespan. We found a bimodal distribution of hosts properties, corresponding to BH binaries suitable to be detected by PTA and the ones detectable only from space missions, as eLISA. Although it has been discussed that the peak of eLISA sources may happen at high z, we show that there must be a population of such sources in the nearby Universe that might show detectable EM signatures, representing an important laboratory for multimessenger astrophysics. We found a weak dependence of galaxy host properties on the binaries occupation, that can be traced back to the BH origin. The combination of the host correlations reported here with the expected EM signal, may be helpful to verify the presence of nearby GW candidates, and to distinguish them from 'regular' intrinsic AGN variability.