System identification uses measurements of a dynamic system’s input and output to reconstruct a mathematical model for that system. These can be mechanical, electrical, physiological, among others. Since most of the systems around us exhibit some form of nonlinear behavior, nonlinear system identification techniques are the tools that will help us gain a better understanding of our surroundings and potentially let us improve their performance. One model that is often used to represent nonlinear systems is the polynomial NARX model, an equation error model where the output is a polynomial function of the past inputs and outputs. That said, a major disadvantage with the polynomial NARX model is that the number of parameters increases rapidly with increasing polynomial order. Furthermore, the polynomial NARX model is a black-box model, and is therefore difficult to interpret. This paper discusses a decoupling algorithm for the polynomial NARX model that substitutes the multivariate polynomial with a transformation matrix followed by a bank of univariate polynomials. This decreases the number of model parameters significantly and also imposes structure on the black-box NARX model. Since a non-convex optimization is required for this identification technique, initialization is an important factor to consider. In this paper the decoupling algorithm is developed in conjunction with several different initialization techniques. The resulting algorithms are applied to two nonlinear benchmark problems: measurement data from the Silver-Box and simulation data from the Bouc-Wen friction model, and the performance is evaluated for different validation signals in both simulation and prediction.