Calculate $\int_{\mathbb{R}^2} uvf(u,v)dudv$ For my probability course I have to work out the following integral and I want to be sure that I did it properly since my last calculus-related course was at least 2 years ago.
The integral is this one:
$$ I = \int_{\mathbb{R}^2} uvf(u,v)\mathrm{d}u\mathrm{d}v$$
$$ f(u,v) = \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}g(u,v)$$
$$ g(u,v) = \exp\left( -\frac{1}{2} \frac{1}{1-\rho^2}\left\{\left(\frac{u}{\sigma_1}\right)^2 -2\rho\left(\frac{u}{\sigma_1}\right)\left(\frac{v}{\sigma_2}\right)+\left(\frac{v}{\sigma_2}\right)^2\right\}\right) $$
First, if I define: $U = \frac{u}{\sigma_1}$ amd $V = \frac{v}{\sigma_2}$ I can rewrite the integral the following way:
$$ I = (\sigma_1\sigma_2)^2\int_{\mathbb{R}^2}UVf(U\sigma_1,V\sigma_2)\mathrm{d}U\mathrm{d}V $$

The current problem is that we can't separate the two integrals because of the term $UV$ in $g(U\sigma_1,V\sigma_2).$ So, that means that I have to find a change of variables so that I can separate $U$ and $V$.
For that purpose, I want to find a base in which $h(U,V) = U^2 - 2\rho UV + V^2$ is free from the term $UV$.
$$h(U,V) = \pmatrix{U & V} \pmatrix{1 & -\rho \\ -\rho & 1} \pmatrix{U \\ V}$$
Its eigenvalues are:
  $$ \lambda_1 = 1+\rho \text{ and } \lambda_2 = {1-\rho}$$
and its eigenvectors are:
  $$ \vec{\lambda_1} = \pmatrix{1 \\ 1} \text{ and } \vec{\lambda_2} = \pmatrix{1 \\ -1}$$
From there I understand that the variable change $X = U + V$ and $Y = U - V$ will do what I want. Indeed,
$$ h(X+Y, X-Y) = 2(X^2(1-\rho) + Y^2(1+\rho)) $$
Since $U$ and $V$ ranged over $\mathbb{R}$ it will be the same for $X+Y$ and $X - Y$.

I switch back to lowercase. The integral can now be written:
$$ J = \frac{2\pi\sqrt{1-\rho^2}}{\sigma_1\sigma_2}I = \int_{\mathbb{R}^2}(x^2 - y^2)\exp\left\{- \frac{x^2}{1+\rho} - \frac{y^2}{1-\rho} \right\}\mathrm{d}x\mathrm{d}y  $$
I can then write,
$$J = J1 - J2 $$
$$ J1 = \int_{\mathbb{R}^2} x^2\exp\left\{- \frac{x^2}{1+\rho} - \frac{y^2}{1-\rho} \right\}\mathrm{d}x\mathrm{d}y$$
$$ J2 = \int_{\mathbb{R}^2} y^2\exp\left\{- \frac{x^2}{1+\rho} - \frac{y^2}{1-\rho} \right\}\mathrm{d}x\mathrm{d}y$$
Before working out $J1$ and $J2$ I remind that:
$$ \int_{\mathbb{R}^2} \exp(-ax^2)\mathrm{d}x = \sqrt{\frac{\pi}{a}} $$
and
$$ \int_{\mathbb{R}^2} x^2\exp(-ax^2)\mathrm{d}x = \sqrt{\frac{\pi}{4a^3}}$$
First let's work out $J1$.
$$\begin{align*}
J1 &= \int_{\mathbb{R}^2} x^2\exp\left\{- \frac{x^2}{1+\rho} - \frac{y^2}{1-\rho} \right\}\mathrm{d}x\mathrm{d}y \\
&= \int_{\mathbb{R}} x^2\exp\left\{- \frac{x^2}{1+\rho}\right\}\mathrm{d}x\int_{\mathbb{R}} \exp\left\{-\frac{y^2}{1-\rho} \right\}\mathrm{d}y \\
&= \frac{\sqrt{(1+\rho)^3\pi}\sqrt{(1-\rho)\pi}}{2}
\end{align*}
$$ 
Then, I obtain:
$$ J2 = \frac{\sqrt{(1-\rho)^3\pi}\sqrt{(1+\rho)\pi}}{2} $$
Thus,
$$ J = \frac{\pi (1-\rho^2) (\sqrt{1-\rho} + \sqrt{1+\rho})}{2}$$

Finally,
  $$ I = \frac{\sigma_1\sigma_2}{4} \left( (1-\rho)\sqrt{1+\rho} + (1+\rho)\sqrt{1-\rho} \right) $$

Is this calculation right? I found that question really mean for a starter in a probability class where the lecturer do not know anything about the background of the students. Maybe there is a shorter and simplier solution but I did not saw it.
 A: My last calculus course was about 5 or 6 years ago, so I can't verify your work.
However, it appears that $I = E[UV]$, so I believe we can proceed as follows.
It appears that $\mu_u = 0$, $\mu_v = 0$, and I will call $\sigma_u = \sigma_1$, $\sigma_2 = \sigma_v$. Hence, we have
$$I = E[UV] = \text{Cov}(U,V)+E[U]E[V] = \text{Corr}(U,V)\text{SD}(U)\text{SD}(V)+0\cdot0 = \rho\sigma_1\sigma_2$$
A: Using $| \cdot|$ to denote the determinant we have
$$f(u,v) = \frac{|\Sigma|^{-1/2}}{2 \pi}\exp\left(-\frac{1}{2}\mathbb{x}'\Sigma^{-1}\mathbb{x}\right),$$
where
$$ \mathbb{x} =\pmatrix{u \\ v},$$
and
$$\Sigma = \pmatrix{\sigma_1^2 & \rho\sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma_2^2} \\  \Sigma^{-1} = \frac{1}{1 - \rho^2}\pmatrix{\sigma_1^{-2} & -\rho\sigma_1^{-1} \sigma_2^{-1} \\ -\rho \sigma_1^{-1} \sigma_2^{-1} & \sigma_2^{-2}}.$$
Since $\Sigma^{-1}$ is a positive definite matrix, there is a lower triangular, nonsingular matrix $C$ such that
$$C'\Sigma^{-1}C = I \\ |C'||\Sigma^{-1}||C| = 1 \implies |C| = |\Sigma|^{1/2} \\ I = I^{-1} = (C'\Sigma^{-1}C)^{-1} =  C^{-1} \Sigma (C')^{-1} \implies CC' = \Sigma.$$
The actual form is
$$C = \pmatrix{ \sigma_1 & 0 \\ \rho \sigma_2  & \sigma_2\sqrt{1-\rho^2}}
.$$
You are trying to find the off-diagonal component of
$$ K =\frac{|\Sigma|^{-1/2}}{2 \pi}\int \mathbb{x} \mathbb{x}'\exp\left(-\frac{1}{2}\mathbb{x}'\Sigma^{-1}\mathbb{x}\right) \, d \mathbb{x}.$$
Make the change of variables $\mathbb{x} = C\mathbb{z}.$  The integral transforms to 
$$K = \frac{|\Sigma|^{-1/2}}{2 \pi}C\int \mathbb{z}' \mathbb{z}\exp\left(-\frac{1}{2}\mathbb{z}'\mathbb{z}\right) \, |C|d \mathbb{z}C' = C\frac{1}{2 \pi}\int \mathbb{z}' \mathbb{z}\exp\left(-\frac{1}{2}\mathbb{z}'\mathbb{z}\right) \, d \mathbb{z}C'.$$
The double integral factors into single integrals, which are easily computed as
$$\frac{1}{2 \pi}\int \mathbb{z}' \mathbb{z}\exp\left(-\frac{1}{2}\mathbb{z}'\mathbb{z}\right) \, d \mathbb{z} = I.$$
Hence, 
$$K = CIC' = CC' = \Sigma.$$
The off-diagonal component is $\rho \sigma_1 \sigma_2.$ 
