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.
