New formulas for moments 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 to a non-negative integer power. Some of the 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. Furthermore, by setting the function inside expectation equal to one, we easily re-derive Isserlis theorem and its generalizations, regarding higher-order moments of a Gaussian random vector. We provide two proofs of the formula, the first being a rigorous proof via mathematical induction. The second proof 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.