  The Galactic Center Excess (GCE) $\gamma$-ray emission detected with the
Large Area Telescope onboard the {\it Fermi Gamma-ray Space Telescope} has been
considered as a possible sign for dark matter (DM) annihilation, but other
possibilities such as the millisecond pulsar (MSP) origin have also been
suggested. As a spectral fitting method, constructed based on properties of
$\gamma$-ray MSPs, has been developed, we apply this method to the study of the
GCE emission for the purpose of probing the MSP origin for the GCE. A number of
$\sim$1660 MSPs can provide a fit to the spectrum of the GCE emission upto
$\sim$10\,GeV, but the higher energy part of the spectrum requires additional
emission components. We further carry out a stacking analysis of 30--500\,GeV
data for relatively nearby $\gamma$-ray MSPs, and the resulting flux upper
limits are still lower than those of the GCE emission. We consider the single
DM annihilation channel $\tau^{+}\tau^{-}$ or channel $b\bar{b}$,or the
combination of the two for comparison, and find they generally can provide
better fits than MSPs. Combination of MSPs plus a DM channel are also tested,
and MSPs plus the DM channel $b\bar{b}$ can always provide better fits.
Comparing this combination case to the pure DM channel $b\bar{b}$, the MSP
contribution is found to be marginally needed.
