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, due to the limited performance of the code phase discriminator, multipath interference (MI), and noise, the code phase estimate can be much different from that of the true first arrival path. In this paper, we derive the statistical ED distribution and power delay spectrum (PDS) of GNSS multipath based on exponential scatterer distribution model (ESDM), and, in parallel, we investigate the ED of the first arrival path, MI and noise distributions, 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 P.681-7 ITU-R channel model. This paper introduces one of the first theoretical analysis and modeling of the GNSS CPEE distributions which can provide insights into the CPEE in multipath environments and are essential to develop algorithms against multipath.