A lower bound error for free-run simulation of the polynomial NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous input) is introduced. The ultimate goal of the polynomial NARMAX is to predict an arbitrary number of steps ahead. Free-run simulation is also used to validate the model. Although free-run simulation of the polynomial NARMAX is essential, little attention has been given to the error propagation to round off in digital computers. Our procedure is based on the comparison of two pseudo-orbits produced from two mathematical equivalent models, but different from the point of view of floating point representation. We apply successfully our technique for three identified models of the systems: sine map, Chua’s circuit and Duffing–Ueda oscillator. This technique may be used to reject a simulation, if a required precision is greater than the lower bound error, increasing the numerical reliability in free-run simulation of the polynomial NARMAX.