  We compute the $S$-factor of the proton-proton ($pp$) fusion reaction using
chiral effective field theory ($\chi$EFT) up to next-to-next-to-leading order
(NNLO) and perform a rigorous uncertainty analysis of the results. We quantify
the uncertainties due to (i) the computational method used to compute the $pp$
cross section in momentum space, (ii) the statistical uncertainties in the
low-energy coupling constants of $\chi$EFT, (iii) the systematic uncertainty
due to the $\chi$EFT cutoff, and (iv) systematic variations in the database
used to calibrate the nucleon-nucleon interaction. We also examine the
robustness of the polynomial extrapolation procedure, which is commonly used to
extract the threshold $S$-factor and its energy-derivatives. By performing a
statistical analysis of the polynomial fit of the energy-dependent $S$-factor
at several different energy intervals, we eliminate a systematic uncertainty
that can arise from the choice of the fit interval in our calculations. In
addition, we explore the statistical correlations between the $S$-factor and
few-nucleon observables such as the binding energies and point-proton radii of
$^{2,3}$H and $^3$He as well as the $D$-state probability and quadrupole moment
of $^2$H, and the $\beta$-decay of $^{3}$H. We find that, with the
state-of-the-art optimization of the nuclear Hamiltonian, the statistical
uncertainty in the threshold $S$-factor cannot be reduced beyond 0.7%.
