A method to estimate the (positive) largest Lyapunov exponent (LLE) from data using interval extensions is proposed. The method differs from the ones available in the literature in its simplicity since it is only based on three rather simple steps. Firstly, a polynomial NARMAX is used to identify a model from the data under investigation. Secondly, interval extensions, which can be easily extracted from the identified model, are used to calculate the lower bound error. Finally, a simple linear fit to the logarithm of lower bound error is obtained and then the LLE is retrieved from it as the third step. To illustrate the proposed method, the LLE is calculated for the following well-known benchmarks: sine map, Rössler Equations, and Mackey-Glass Equations from identified models given in the literature and also from two identified NARMAX models: a chaotic jerk circuit and the tent map. In the latter, a Gaussian noise has been added to show the robustness of the proposed method.