# Covariance of polynomials of random normal variables

$\newcommand{\Cov}{\operatorname{Cov}}$If $X$ and $Y$ are random variables with a bivariate normal distribution and:

• $X\sim\mathcal{N}(\mu_X,\sigma_X^2)$
• $Y\sim\mathcal{N}(\mu_Y,\sigma_Y^2)$
• $\Cov(X,Y)\neq0$

May I compute $\Cov(X^m,Y^n)$ for arbitrary positive integers $m$ and $n$?

• I don't think you can compute this object for arbitrary m and n. You would need to know more information about the random variables X, Y in the form of their joint distribution. Even if you assume that the covariance is well known, you still can't find the higher moments of the distribution. Aug 7, 2018 at 22:32
• If I assume instead that $(X,Y)$ have a bivariate normal distribution whose variance-covariance matrix is known, does the question become answerable? Aug 7, 2018 at 22:52
• This question fails to state that $X,Y$ are JOINTLY normal. If that is assumed, then, since their covariance is $0,$ they are independent. Aug 8, 2018 at 0:32
• It is sloppy notation to write $X\sim N(\mu_x, \sigma^2_x)$ instead of $X\sim N(\mu_X, \sigma^2_X).$ And in many contexts when working with only slightly more involved problems of this kind, this sort of confusion can paralyze you. Aug 8, 2018 at 0:37
• In your comment you say $X,Y$ have a bivariate normal distribution, but in your question you say only that each one separately is normally distributed and their covariance is $0.$ I can show you examples of a distribution of a pair $(X,Y)$ in which each is normally distributed and their covariance is $0$ and they are NOT JOINTLY normally distributed. But if we assume bivariate (thus joint) normality, then the answer is easy. See my answer below. Aug 8, 2018 at 0:44

$$\newcommand{\Cov}{\operatorname{Cov}}$$Given the new background that the OP provided to his question (most general normal bivariate distribution for variables (X_1, X_2) with non-zero covariance), it is possible to find a general formula for $$\Cov(X_1^m, X_2^n)$$ as follows:

Consider the generating function $$\phi(t_1,t_2)$$ for the bivariate distribution given here, in equation (57) (a clean derivation is given for it so I won't repeat it). Then:

$$E(X_1^m X_2^n)=(-i)^{m+n}(\frac{\partial}{\partial t_1})^m(\frac{\partial}{\partial t_2})^n\phi(t_1,t_2)\Big|_{(t_1, t_2)=(0,0)} \tag 1$$

After some tedious and careful algebra we wish to write the generating function $$\phi(t_1,t_2)=e^{it_1\mu_1+it_2\mu_2}e^{-\frac{1}{2}(\sigma_1^2t_1^2+2\rho\sigma_1t_1\sigma_2t_2+\sigma_2^2t_2^2)}$$ in the form:

$$\phi(t_1,t_2)=e^{-a}e^{-b(t_2+c)^2}e^{-d(t_1+gt_2+h)^2}\tag2$$

which is possible for the values:

$$a=\frac{\mu_1^2}{2\sigma_1^2}+\frac{(\mu_2-\mu_1\rho \frac{\sigma_2}{\sigma_1})^2}{2\sigma^2_2(1-\rho^2)}, \hspace{0.2cm}b=\frac{1}{2}(1-\rho^2)\sigma_{2}^2, \hspace{0.2cm}c=-i\frac{\mu_2-\mu_1\rho \frac{\sigma_2}{\sigma_1}}{2\sigma^2_2(1-\rho^2)}, \hspace{0.2cm}d=\frac{\sigma_1^2}{2},\hspace{0.2cm} g=\frac{\rho\sigma_2}{\sigma_1},\hspace{0.2cm}h=-i\frac{\mu_1}{\sigma_1^2}$$

We change variables $$y_2=\sqrt{b}(t_2+c), y_1=\sqrt{d}(t_1+gt_2+h)$$ and we find:

\begin{align}\frac{\partial}{\partial t_1}&=\sqrt{d}\frac{\partial}{\partial y_1}\\ \frac{\partial}{\partial t_2}&=g\sqrt{d}\frac{\partial}{\partial y_1}+\sqrt{b}\frac{\partial}{\partial y_2}\end{align}

Substitute into (1) for the result:

$$E(X_1^mX_2^n)=(-i)^{m+n}(\sqrt{d}\frac{\partial}{\partial y_1})^m(g\sqrt{d}\frac{\partial}{\partial y_1}+\sqrt{b}\frac{\partial}{\partial y_2})^n e^{-a}e^{-y_1^2}e^{-y^2_2}\Big|_{(y_1, y_2)=(c\sqrt{b},h\sqrt{d})}$$

Expanding the parentheses and using the Rodrigues formula for Hermite polynomials we get:

$$E(X_1^mX_2^n)=e^{-a-bc^2-dh^2}i^{m+n}g^n d^{\frac{m+n}{2}}\sum_{k=0}^n\frac{n!}{k!(n-k)!}\Big(\frac{\sqrt{b}}{g\sqrt{d}}\Big)^{n-k}H_{n-k}(h\sqrt{d})H_{m+k}(c\sqrt{b})$$

The expectation values $$E(X_1^m)$$ and $$E(X_2^n)$$ can be calculated by setting $$n=0$$, $$m=0$$ in the general formula respectively. Also, the imaginary units in the argument of the Hermite polynomials are precisely cancelled out by the imaginary prefactor, so that the final result is real.

EDIT: According to the analysis by @machfour, there is indeed an error in the expression of $$c$$, which has been corrected for. Due to this, the formula simplifies to the final expression:

$$\small{E(X_1^mX_2^n)=\Big(\frac{i}{\sqrt{2}}\Big)^{m+n}(\rho\sigma_2)^n \sigma_1^m\sum_{k=0}^n {n\choose k} \Big(\sqrt{\frac{1}{\rho^2}-1}\Big)^{n-k}H_{n-k}\Big(\frac{-i\mu_1}{\sigma_1\sqrt{2}}\Big)H_{m+k}\Big(-i\frac{\mu_2\sigma_1-\mu_1\rho}{\sqrt{2}\sigma_1\sigma_2\sqrt{1-\rho^2}}\Big)}$$

• This makes the problem far more complicated than it is, unless you don't know that JOINT normality plus uncorrelatedness entails independence. Aug 9, 2018 at 2:08
• The OP has stated that the two variables are correlated since their covariance is non-zero. The OP also updated their request with a joint probability that reflects this. This problem is valid and I answered with a calculation that addresses it. The formula is simple and a finite polynomial in all variables except for $\rho$ so I dont see how this problem is complicated. Aug 9, 2018 at 3:43
• I don't think the OP gave a good description of what the problem was in it's initial wording, but then they changed their mind. It seems to me that to them quadratic correlations are more important than solving the totally trivial independent variable problem. Aug 9, 2018 at 3:48

I don't have enough reputation to comment, but I think there is a typo/error in DinosaurEgg's answer. I repeated the derivation of (2) in that answer and obtained a slightly different value for $$c$$ which was off by a factor of 2: $$\begin{equation*} c = -i \frac{\left(\mu_2 - \mu_1\rho \frac{\sigma_2}{\sigma_1}\right)}{\sigma_2^2(1 - \rho^2)} \end{equation*}$$ while the values I obtained for $$a$$, $$b$$, $$d$$, $$g$$, and $$h$$ were the same.

Regardless of the correct value, notice also that since there is no constant term in the exponent of the original characteristic function $$\phi(t_1, t_2) = \exp\left(i t_1 \mu_1 + i t_2 \mu_2 - \frac{1}{2}(t_1^2 \sigma_1^2 + 2t_1 t_2 \sigma_1 \sigma_2 \rho + t_2^2 \sigma_2^2) \right),$$ examining (2) suggests that we should have $$- a - bc^2 - dh^2 = 0$$ as an identity. Thus the final formula simplifies to \begin{align*} E(X_1^mX_2^n) &=i^{m+n}g^n d^{\frac{m+n}{2}}\sum_{k=0}^n\frac{n!}{k!(n-k)!}\Big(\frac{\sqrt{b}}{g\sqrt{d}}\Big)^{n-k}H_{n-k}(h\sqrt{d})H_{m+k}(c\sqrt{b}) \\ &=\left(\frac{i}{\sqrt{2}}\right)^{m+n} \sigma_1^m (\rho \sigma_2)^n \sum_{k=0}^n\frac{n!}{k!(n-k)!}\Big(\sqrt{1/\rho^2 - 1}\Big)^{n-k}H_{n-k}(-\frac{i \mu_1}{\sqrt{2} \sigma_1})H_{m+k}(c\sqrt{b}) \\ \end{align*}

You have said$$\newcommand{\Cov}{\operatorname{Cov}}$$ \begin{align} & X\sim N(\mu_X,\sigma^2_X), \\[4pt] & Y \sim N(\mu_Y,\sigma^2_Y), \\[4pt] & \Cov(X,Y) = 0. \end{align} That falls short of specifying the joint distribution of $$(X,Y).$$ If it were further specified that $$X,Y$$ are jointly normally distributed, then the covariance can be $$0$$ only if $$X,Y$$ are independent. If $$X,Y$$ are independent, then so are $$X^m,Y^n,$$ so their covariance is also $$0.$$

Here is a simple example: Suppose $$X\sim N(\mu_X, \sigma_X^2),$$ let $$Z = \text{the “z-score''} = (X-\mu_X)/\sigma_X,$$ and independently of $$X$$ you toss a coin. Then let $$Y = \mu_Y \pm \sigma_Y Z,$$ where the choice between $$\text{“}\pm\text{''}$$ is determined by the coin toss. Then $$X,Y$$ have covariance $$0$$ and have just the distributions you specified in the question, but they are NOT JOINTLY normally distributed and not independent.

But if you assume joint normality, which in the case of two random variables means bivariate normality, then the answer is just as in the first paragraph above.

• Thanks, but my original question states that Cov(X,Y) is not equal to zero. I edited it to specify a joint normal distribution. Aug 8, 2018 at 6:05
• @Mathieu : In that case, my first paragraph above answers your question, and the answer by E.Malkin misses the point at best. There is no need to know the variance of $X^m$ or even the expected value of $X^m$ in order to answer your question. Aug 8, 2018 at 12:41
• Sorry but I must be missing something: my question states that Cov(X,Y) is not zero and your first paragraph above states the opposite. Not sure what to make of that. Aug 8, 2018 at 13:14
• @MichaelHardy I also don't see how your first paragraph answers the question and could you please explain how I've missed the point ("at best")? Knowing $E(X^m)$ means that finding $\Cov (X^m,Y^n)$ is equivalent to finding $E(X^mY^n)$, which surely isn't invalid? Aug 8, 2018 at 22:04
• @Malkin : If $X,Y$ are not just separately normal but jointly normal, then having zero covariance entails that they are independent. If $X,Y$ are independent, then $X^m,Y^n$ are independent. If $X^m, Y^n$ are independent, then $\operatorname E(X^m Y^n) = \operatorname E(X^m) \operatorname E(Y^n).$ Therefore $$\operatorname{cov}(X^m, Y^n) = \operatorname E(X^m Y^n) - \operatorname E(X^m) \operatorname E(Y^n) = \operatorname E(X^m) \operatorname E(Y^n) - \operatorname E(X^m) \operatorname E(Y^m) = 0.$$ Aug 9, 2018 at 2:04