This study introduces a novel theoretical framework for analyzing heteroscedastic Gaussian processes (HGPs) that identify unknown systems in a data-driven manner. Although HGPs effectively address the heteroscedasticity of noise in complex training datasets, calculating the exact posterior distributions of the HGPs is challenging, as these distributions are no longer multivariate normal. This study derives the exact means, variances, and cumulative distributions of the posterior distributions. Furthermore, the derived theoretical findings are applied to a chance-constrained tracking controller. After an HGP identifies an unknown disturbance in a plant system, the controller can handle chance constraints regarding the system despite the presence of the disturbance.