<p>In this paper, we develop a methodology of Multivariate Generalized Linear Models (MGLMs) for Markov kernels with known as well as unknown link and variance functions. Within the framework of multivariate exponential families, the marginal distribution involves interesting properties of a Markov kernel or a mixture distribution like normal and Poisson models. Theoretical results are first obtained if these link and (generalized) variance functions have closed forms for some complex but practical models as the multivariate normal inverse Gaussian distribution. In order to reduce the assumptions in a full parametric MGLM, we assume both link and variance functions to be unknown but smooth. These functions are then estimated through the multivariate classical Gaussian kernel. Therefore, we obtain a three-level approach to this semiparametric model from link function, variance function, and the vector of regression coefficients in the linear predictor of the model. Additionally, we consider total deviances criteria, for large responses, covariates, and sample sizes of datasets. Finally, through simulations and applications, we corroborate that the resulting parameter estimates are satisfactory for both parametric and semiparametric modes, despite the complexity of the considered models.</p>