In a modern power system, there is often large difference in the decay speeds of transients. This could lead to numerical problems such as heavy simulation burden and singularity when the traditional methods are used to estimate the stability region of such a dynamic system with saturation nonlinearities. To overcome these problems, a reduced-order method, based on the singular perturbation theory, is suggested to estimate the stability region of a singular system with saturation nonlinearities. In the reduced-order method, a low-order linear dynamic system with saturation nonlinearities is constructed to estimate the stability region of the primary high-order system so that the singularity is eliminated and the estimation process is simplified. In addition, the analytical foundation of the reduction method is proven and the method is validated using a test power system with 3 buses and 5 machines.