# 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.
• 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.
– whuber
Jul 25, 2022 at 15:48
• Thank you! I would upvote if I could. Jul 25, 2022 at 15:59
• Your question title assumes a bivariate Normal distribution, but the content of your question does not. Jul 25, 2022 at 22:54

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 * 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.

• 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. Jul 28, 2022 at 20:39