We present a variational Bayesian identification procedure for polynomial NARMAX models based on message passing on a factor graph. Message passing allows us to obtain full posterior distributions for regression coefficients, precision parameters and noise instances by means of local computations distributed according to the factorization of the dynamic model. The posterior distributions are important to shaping the predictive distribution for outputs, and ultimately lead to superior model performance during 1-step ahead prediction and simulation.