New formulas for moments and functions of the multivariate normal distribution extending Stein's lemma and Isserlis theorem

We prove a formula for the evaluation of expectations containing a scalar function of a Gaussian random vector multiplied by a product of the random vector components, each one raised at a non-negative integer power. Some powers could be of zeroth-order, and, for expectations containing only one vector component to the first power, the formula reduces to Stein's lemma for the multivariate normal distribution. On the other hand, by setting the said function inside expectation equal to one, we easily derive Isserlis theorem and its generalizations, regarding higher order moments of a Gaussian random vector. We provide two proofs of the formula, with the first being a rigorous proof via mathematical induction. The second is a formal, constructive derivation based on treating the expectation not as an integral, but as the consecutive actions of pseudodifferential operators defined via the moment-generating function of the Gaussian random vector.