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.
  • 1
    $\begingroup$ Start by simplifying the problem. You may assume both means are zero because $$E[Y(X-\mu_X)^p]=E[(Y-\mu_Y)(X-\mu_X)^p]+\mu_YE[(X-\mu_X)^p].$$ Then you can change the units of measure of $X$ and $Y$ (leading to a simple, predictable change in this expectation) so that only $\rho$ is left to figure in the calculation. Finally, $Y$ can be expressed as a linear combination of $X$ and a variable uncorrelated with it. With these simplifications, the answer can be found as you propose, as shown at stats.stackexchange.com/questions/315983. $\endgroup$
    – whuber
    Jul 25, 2022 at 15:48
  • $\begingroup$ Thank you! I would upvote if I could. $\endgroup$
    – SeanBrooks
    Jul 25, 2022 at 15:59
  • 1
    $\begingroup$ Your question title assumes a bivariate Normal distribution, but the content of your question does not. $\endgroup$
    – wolfies
    Jul 25, 2022 at 22:54

1 Answer 1


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.

  • 1
    $\begingroup$ I appreciate this answer, after moving forward with my work, I am seeing that expressing normally distributed random variables in terms of independent standard normal distributions is very useful in general. $\endgroup$
    – SeanBrooks
    Jul 28, 2022 at 20:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.