  Weak gravitational lensing, or weak lensing, is one of the most powerful
probes for dark matter and dark energy science, although it faces increasing
challenges in controlling systematic uncertainties as \edit{the statistical
errors become smaller}. The Point Spread Function (PSF) needs to be precisely
modeled to avoid systematic error on the weak lensing measurements. The weak
lensing biases induced by errors in the PSF model second moments, i.e., its
size and shape, are well-studied. However, Zhang et al. (2021) showed that
errors in the higher moments of the PSF may also be a significant source of
systematics for upcoming weak lensing surveys. Therefore, the goal of this work
is to comprehensively investigate the modeling quality of PSF moments from the
$3^{\text{rd}}$ to $6^{\text{th}}$ order, and estimate their impact on
cosmological parameter inference. We propagate the \textsc{PSFEx} higher
moments modeling error in the HSC survey dataset to the weak lensing
\edit{shear-shear correlation functions} and their cosmological analyses. We
find that the overall multiplicative shear bias associated with errors in PSF
higher moments can cause a $\sim 0.1 \sigma$ shift on the cosmological
parameters for LSST Y10. PSF higher moment errors also cause additive biases in
the weak lensing shear, which, if not accounted for in the cosmological
parameter analysis, can induce cosmological parameter biases comparable to
their $1\sigma$ uncertainties for LSST Y10. We compare the \textsc{PSFEx} model
with PSF in Full FOV (\textsc{Piff}), and find similar performance in modeling
the PSF higher moments. We conclude that PSF higher moment errors of the future
PSF models should be reduced from those in current methods to avoid a need to
explicitly model these effects in the weak lensing analysis.
