  This paper develops recurrence relations for integrals that relate the
density of multivariate extended skew-normal (ESN) distribution, including the
well-known skew-normal (SN) distribution introduced by Azzalini and Dalla-Valle
(1996) and the popular multivariate normal distribution. These recursions offer
a fast computation of arbitrary order product moments of the multivariate
truncated extended skew-normal and multivariate folded extended skew-normal
distributions with the product moments as a byproduct. In addition to the
recurrence approach, we realized that any arbitrary moment of the truncated
multivariate extended skew-normal distribution can be computed using a
corresponding moment of a truncated multivariate normal distribution, pointing
the way to a faster algorithm since a less number of integrals is required for
its computation which result much simpler to evaluate. Since there are several
methods available to calculate the first two moments of a multivariate
truncated normal distribution, we propose an optimized method that offers a
better performance in terms of time and accuracy, in addition to consider
extreme cases in which other methods fail. The R MomTrunc package provides
these new efficient methods for practitioners.
