Although models based on the Gaussian process has been broadly used for flexible nonparametric modeling, they are not suitable for modeling abrupt changes of the smoothness of the target function and relationships with heteroscedastic errors. The heteroscedastic Gaussian process (HeGP) regression attempts to overcome these limitations by assuming that residual variances of the regression model vary over covariates. We here generalize the idea of HeGP so that it is applicable to not only regression problems but also classification and state-space models. We let the Gaussian process be coupled with a covariate-induce precision matrix process that takes a mixture form so as to model the heteroscedastic covariance function over covariates. To cope with excessive computational burdens from sampling, we resort to the variational inference for the posterior approximation in evaluating the posterior predictive model and in training via the EM algorithm with closed-form M-step updates for evaluating the heteroscedastic covariance function. Our model works consistently on the multivariate responses, even if they are of different types (either continuous or categorical). We demonstrate its advantages by both simulations and applications to real-data examples from climatology.