System Identification (SI) is a control engineering discipline concerned with the discovery of mathematical models based on dynamic measurements collected from the system. It is an important discipline in the construction and design of controllers, as SI can be used for understanding the properties of the system as well as to forecast its behavior under certain past inputs and/or outputs. The NARMAX model and its derivatives (Nonlinear Auto-Regressive with Exogenous Inputs (NARX) and Nonlinear Auto-Regressive Moving Average (NARMA)) are powerful, efficient and unified representations of a variety of nonlinear models. The identification process of NARX/NARMA/NARMAX involves structure selection and parameter estimation, which can be simultaneously performed using the widely accepted Orthogonal Least Squares (OLS) algorithm. Several criticisms have been directed towards OLS for its tendency to select excessive or sub-optimal terms. The suboptimal selection of regressor terms leads to models that are non-parsimonious in nature. This thesis proposes the application of a stochastic optimization algorithm called Binary Particle Swarm Optimization algorithm for structure selection of polynomial NARX/NARMA/NARMAX models. The algorithm searches the solution space by selecting various model structures and evaluating its fitness. A MySQL database was created to analyze the optimization results and speed up computations of the optimization algorithm. The proposed optimization algorithm was tested on several benchmark datasets, namely the Direct Current Motor (DCM), Mackey-Glass Differential Equation (MG) and Flexible Robot Arm (FRA). The DCM motor was the least complication dataset, followed by the FRA (medium complexity) and MG (most complexity). The results suggest that the proposed method can reduce the number of correlation violations down to between 28.57% and 69.23% at the expense of increased model size (requirement of additional regressor terms to explain the behavior of the system)