It was found in the Fermi-LAT data that there is an extended $gamma$-ray excess in the Galactic center region. The proposed sources to be responsible for the excess include the dark matter annihilation or an astrophysical alternative from a population of millisecond pulsars (MSPs). Whether or not the MSP scenario can explain the data self-consistently has very important implications for the detection of particle dark matter, which is however, subject to debate in the literature. In this work we study the MSP scenario in detail, based on the detected properties of the MSPs by Fermi-LAT. We build a model of the Milky Way MSPs which can reproduce the $gamma$-ray properties of the Fermi-LAT MSPs, and derive the intrinsic luminosity function of the MSPs. The model is then applied to a bulge population of MSPs. We find that the extended $gamma$-ray excess can be well explained by the bulge MSPs without violating the detectable flux distribution of MSPs by Fermi-LAT. The spatial distribution of the bulge MSPs as implied by the distribution of low mass X-ray binaries follows a $r^{-2.4}$ profile, which is also consistent with the $gamma$-ray excess data. We conclude that the MSP model can explain the Galactic center $gamma$-ray excess self-consistently, satisfying all the current observational constraints.