ABSTRACT We present an analytic model to describe the supermassive black hole binary (SMBHB) merger rate in the Universe with astrophysical observables: galaxy stellar mass function, pair fraction, merger time-scale, and black hole–host galaxy relations. We construct observational priors and compute the allowed range of the characteristic spectrum hc of the gravitational wave background (GWB) to be 10−16 < hc < 10−15 at a frequency of f = 1 yr−1. We exploit our parametrization to tackle the problem of astrophysical inference from pulsar timing array (PTA) observations. We simulate a series of upper limits and detections and use a nested sampling algorithm to explore the parameter space. Corroborating previous results, we find that the current PTA non-detection does not place significant constraints on any observables; however, either future upper limits or detections will significantly enhance our knowledge of the SMBHB population. If a GWB is not detected at a level of hc(f = 1 yr−1) = 10−17, our current understanding of galaxy and SMBHB mergers is disfavoured at a 5σ level, indicating a combination of severe binary stalling, overestimating of the SMBH–host galaxy relations, and extreme dynamical properties of merging SMBHBs. Conversely, future detections of a Square Kilometre Array (SKA)-type instrument will allow to constrain the normalization of the SMBHB merger rate in the Universe, the time between galaxy pairing and SMBHB merging, the normalization of the SMBH–host galaxy relations and the dynamical binary properties, including their eccentricity and density of stellar environment.