  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 timescale and black hole - host
galaxy relations. We construct observational priors and compute the allowed
range of the characteristic spectrum $h_c$ of the gravitational wave background
(GWB) to be $10^{-16}<h_c<10^{-15}$ at a frequency of $f=1/{\rm yr}$. 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
$h_c(f=1/{\rm yr})=10^{-17}$, our current understanding of galaxy and SMBHB
mergers is disfavoured at a $5\sigma$ level, indicating a combination of severe
binary stalling, over-estimating of the SMBH -- host galaxy relations, and
extreme dynamical properties of merging SMBHBs. Conversely, future detections
of a Square Kilometre Array (SKA)-type array 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.
