Bayesian regression, a nonparametric identification technique with several appealing features, can be applied to the identification of NARX (nonlinear ARX) models. However, its computational complexity scales as $O(N^3)$ where $N$ is the data set size. In order to reduce complexity, the challenge is to obtain fixed-order parametric models capable of approximating accurately the nonparametric Bayes estimate avoiding its explicit computation. In this work we derive, optimal finite-dimensional approximations of complexity $O(N^2)$ focusing on their use in the parametric identification of NARX models.