We propose a novel nonparametric approach to ARMAX identification with missing data relying upon recent work on predictor estimation via Gaussian regression. The Bayesian setup allows one to compute explicitly an input-output marginal density where the model dependence has been integrated out. This turns out to be a key step in facilitating the imputation of missing variables. Thus, this approach has the advantage that no classical ??model selection?? (or model order estimation) has to be performed. Model ??complexity?? is described by means of hyperparameters which are estimated as part of the identification procedure. The new approach is shown to perform better than standard prediction error methods (PEM), also when the full data set is made available to the latter, in terms of both predictive capability on new data and accuracy in predictor coefficients reconstruction.