Correct modeling of light-emitting diode (LED) nonlinearity is vital for evaluating the performance of advanced modulation formats and nonlinear distortion mitigation methods in visible light communication systems. However, the models that can be found in the literature are either too complicated (e.g., the rate equation model), incorrect (the static polynomial model), or too simplified (the Wiener and Hammerstein models). In this paper, we demonstrate that rigorous frequency-domain factorization of the Volterra kernel up to the 2nd-order derived from the LED rate equations leads to an evenly simple but highly accurate model of the LED, which can be represented by a concatenation of 3 blocks: a filter, squaring, and a second filter. The model is parametrized by the 3 dB bandwidth of the LED, the cutoff frequency of the second filter, the optical conversion efficiency, and the nonlinearity factor. We also propose a two-point single-tone method for finding the mentioned parameters. We provide experimental evidence for the correctness of this block structure in LEDs of five different types and also disprove the previously mentioned simplified models. Finally, we validate the model for transmission of single carrier (PAM) and multicarrier (OFDM) signals by comparing it with the rate equation and the basic linear models. A mean squared error reduction of up to 5 dB compared with the linear model is experimentally observed.