We present a Bayesian identification procedure for polynomial NARMAX models based on message passing on a factor graph. Our probabilistic perspective lets us treat the noise instances as latent random variables and infer posterior distributions. These have the structure of precision-weighted prediction errors, as opposed to the point prediction errors obtained by classical NARMAX estimators. We empirically show that precision-weighting improves the estimates of the coefficients and ultimately the performance of the model during 1-step ahead prediction and simulation.