Equivalent energy function method is practical and efficient in power system stochastic production simulation. However, it may result in low accuracy when applied to calculate the loss of load probability (LOLP) index. Using the step-shaping feature of the load duration curve for an actual power system, the abscissas of inflexion are rounded to integers by introducing reference unit. Furthermore, a sparse recursive convolution method is proposed for power systems stochastic production simulation. By describing the load duration curve with variable-width rectangles, it makes full use of the sparse feature of the inflection point distribution in the integer space. The computational time required by the proposed method is insensitive to the values of reference unit. It is more applicable for a large scale power system that requires a higher accuracy of LOLP index. The accuracy and efficiency of the proposed method are verified by simulations on the IEEE-RTS79 system and its modified and expanded systems. This work is supported by National Natural Science Foundation of China (No. 51377143).