In the area of black-box identication, NARMAX models are of great interest. The main diculty faced when working with such models is the selection of the correct structure to represent the underlying system in the data. Orthogonal Least Squares (OLS) methods are widely used for this task, however, there are systems with a high degree of non-linearity and long term dependencies, which makes the use of traditional OLS methods computationally impracticable. In this sense, this paper studies the use of Multi-Gene Genetic Programing (MGGP) together with the traditional OLS method to increase the search space and turn the structure selection practicable for average performance computer. It is shown that, in real-life problem data, the algorithm can nd better models than previous works' models. The MGGP found a model for a hydraulic pumping system with a better one-step-ahead prediction error (0:058 mlc2 against 0:070 mlc2) using PEM technique and better free-run simulation error (0:997 mlc2 against 1:120 mlc2) using SEM technique. The MGGP found a model with such a degree of non-linearity and maximum input-output lags that totalizes 142505 candidate terms for traditional OLS analysis, which is impracticable for average performance computers.