How to compute this moment of a bivariate normal distribution? Consider that $X \sim N(\mu_X, \sigma_X^2)$ and $Y \sim N(\mu_Y, \sigma^2_Y)$ and $\text{Cov}[X, Y] = \sigma_{XY}$ where $\sigma_{XY} \ne 0$.
How can an expression for $\text{E}[Y(X-\mu_X)^p]$ in terms of $p$, $\mu_X$, $\mu_Y$, $\sigma_X$, $\sigma_Y$ and $\sigma_{XY}$ be found?
So far:

*

*I have tried to use a moment-generating function but $e^{tY(X-\mu_X)}$ does not separate the $X$ and $Y$ nicely.

*I did some numerical tests to see if it has an obvious relationship to the $p$-th central moment of $X$ with no luck.

 A: The analysis at https://stats.stackexchange.com/a/71303/919 suggests simplifying the question by expressing $X$ and $Y$ in terms of two independent standard Normal variables.
An easy way to do this begins by standardizing $X$ to express it as
$$X=\mu_X + \sigma_X Z$$
for a standard Normal variable $Z.$ We know $Y$ can be expressed in terms of $Z$ and an independent standard Normal variable $W$ as
$$Y = \mu_Y + \alpha Z + \beta W$$
Because $Z$ and $W$ are independent, the covariance of $(X,Y)$ depends only on their $Z$ coefficients, telling us that
$$\rho\sigma_X\sigma_Y = \operatorname{Cov}(X,Y)  = \operatorname{Cov}(\mu_X + \sigma_X Z, \mu_Y + \alpha Z + \beta W)= \sigma_X\alpha.$$
The solution is $\alpha = \rho\sigma_Y.$ That's all we need to know.
To answer the question, apply the basic properties of expectation to compute
$$\begin{aligned}
E[Y(X-\mu_X)^p] &= E[(\mu_Y + \sigma_Y\rho Z + \sigma_Y\sqrt{1-\rho^2}W)(\sigma_X Z)^p]\\
&=\mu_Y\sigma_X^p E[Z^p] + \rho\sigma_Y\sigma_X^p E[Z^{p+1}] + (\text{constants})E[WZ^p]\\
& = \mu_Y \sigma_X^p E[Z^p] + \rho\sigma_Y\sigma_X^p E[Z^{p+1}].
\end{aligned}\tag{*}$$
As found at https://stats.stackexchange.com/a/176814/919, when $p=2k$ is a positive even integer,
$$E[Z^p] = E[Z^{2k}] = \frac{(2k)!}{k! 2^k}$$
and $E[Z^{p+1}] = 0.$  Otherwise, when $p=2k-1$ is a positive odd integer,
$$E[Z^{p+1}] = E[Z^{2k}] = \frac{(2k)!}{k! 2^k}$$
and $E[Z^p]=0.$
Plug these into $(*)$ for the answer.
If you would like to perform numerical checks, here is an R implementation.
#
# Normal moment of non-negative integral order `p`.
#
norm.moment <- function(p) {
  ifelse(p %% 2 == 0, 
         ifelse(p == 0, 1, exp(lfactorial(p) - lfactorial(p/2) - (p/2) * log(2))), 0)
}
#
# The formula for the expectation in the question.
#
f <- function(p, mu, Sigma) {
  rho <- Sigma[1,2] / sqrt(prod(diag(Sigma)))
  Sigma[1,1] ^ (p/2) * (mu[2] * norm.moment(p) + rho * sqrt(Sigma[2,2]) * norm.moment(p+1))
}

You might, for instance, compare these values to a Monte-Carlo estimate:
library(MASS)                      # Exports mvrnom()
n <- 1e6                           # Number of (X,Y) to simulate
mu <- c(-2, 3)                     # Mean vector
Sigma <- matrix(c(4, 2, 2, 5), 2)  # Covariance matrix (symmetric!)
X <- mvrnorm(n, mu, Sigma)         # Simulated data

p <- 4
(mean(X[,2] * (X[,1] - mean(X[,1]))^p)) # Monte-Carlo estimate
f(p, mu, Sigma)                         # What the formula says

The Monte-Carlo result in this case will be close to $144,$ the value given by the formula.
