In multipath environments, a global navigation satellite system (GNSS) receiver can obtain the most correct code phase estimate from the resolved first arrival path, which is expected to have the smallest excess delay (ED). However, because of the limited performance of the code phase discriminator, multipath interference (MI), and noise, the code phase estimate can be different from that of the true first arrival path. In this paper, we derive the statistical ED distribution and power delay spectrum of GNSS multipath components based on exponential scatterer distribution model (ESDM). In parallel, we investigate the ED distributions of the first arrival path, MI, and noise to develop mathematical expressions for the code phase estimation error (CPEE) distribution for wide, narrow, and strobe correlators in various multipath channels. The mathematical models of CPEE distributions have good match with the ESDM-based CPEE distributions and the CPEE distributions obtained from Monte Carlo simulations using the International Telecommunications Union Recommendations Section recommendation P.681-7 channel model. This paper introduces one of the first theoretical analyses and models of the GNSS CPEE distributions, which can provide insights into the CPEE in multipath environments and are essential to develop algorithms against multipath distortion.