Matrix valued positive definite kernels related to the generalized Aitken's integral for Gaussians

We introduce a method to construct general multivariate positive definite kernels on a nonempty set XX that employs a prescribed bounded completely monotone function and special multivariate functions on XX. The method is consistent with a generalized version of Aitken's integral formula for Gaussians. In the case in which XX is a cartesian product, the method produces nonseparable positive definite kernels that may be useful in multivariate interpolation. In addition, it can be interpreted as an abstract multivariate version of the well-established Gneiting's model for constructing space-time covariances commonly highly cited in the literature. Many parametric models discussed in statistics can be interpreted as particular cases of the method.