  We present a new Bayesian hierarchical model (BHM) named Steve for performing
type Ia supernova (SNIa) cosmology fits. This advances previous works by
including an improved treatment of Malmquist bias, accounting for additional
sources of systematic uncertainty, and increasing numerical efficiency. Given
light curve fit parameters, redshifts, and host-galaxy masses, we fit Steve
simultaneously for parameters describing cosmology, SNIa populations, and
systematic uncertainties. Selection effects are characterised using Monte-Carlo
simulations. We demonstrate its implementation by fitting realisations of SNIa
datasets where the SNIa model closely follows that used in Steve. Next, we
validate on more realistic SNANA simulations of SNIa samples from the Dark
Energy Survey and low-redshift surveys. These simulated datasets contain more
than $60\,000$ SNeIa, which we use to evaluate biases in the recovery of
cosmological parameters, specifically the equation-of-state of dark energy,
$w$. This is the most rigorous test of a BHM method applied to SNIa cosmology
fitting, and reveals small $w$-biases that depend on the simulated SNIa
properties, in particular the intrinsic SNIa scatter model. This $w$-bias is
less than $0.03$ on average, less than half the statistical uncertainty on
$w$.These simulation test results are a concern for BHM cosmology fitting
applications on large upcoming surveys, and therefore future development will
focus on minimising the sensitivity of Steve to the SNIa intrinsic scatter
model.
