how to do this integral: $$ \int_{0}^{\infty} \int_{0}^{\infty} x y \phi(x, y) dx dy$$
where $ \phi(x,y)$ is a general pdf of bivariate normal distribution, that is:
$$\phi(x,y) = \frac1{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}}exp\left({{-\frac1{2(1-\rho^2)}} \left[\frac{(x-\mu_x)^2}{\sigma_x^2}+\frac{(y-\mu_y)^2}{\sigma_y^2}-\frac{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x\sigma_y}\right]}\right)$$
$\rho$ is the correlation between x and y.
The answer is a complicated function of the parameters $(\mu_x,\sigma_x,\mu_y,\sigma_y,\rho)$. To see why there is no reason to expect a nice formula, consider the much simpler problem of computing $$ I=\int_0^{+\infty}x\phi(x)\mathrm dx, $$ where $\phi$ denotes the PDF of a normal random variable with mean $\mu$ and variance $\sigma^2$. Then, using the fact that $X=\mu+\sigma U$ where $U$ is standard normal, one sees that $$ I=E[\mu+\sigma U;\mu+\sigma U\gt0]=\mu P[U\gt-\mu/\sigma]+\sigma E[U;U\gt-\mu/\sigma], $$ thus $$ I=\mu\Phi\left(\frac{\mu}\sigma\right)+\sigma\varphi\left(\frac{\mu}\sigma\right), $$ where $\varphi$ and $\Phi$ are the PDF and CDF of a standard normal random variable, that is, $$ \varphi(x)=\frac1{\sqrt{2\pi}}\mathrm e^{-x^2/2},\qquad\Phi(u)=\int_{-\infty}^u\varphi(x)\mathrm dx=P[U\leqslant u]. $$ This yields the $\rho=0$ case of the two-dimensional problem and suggests that the formula for $\rho\ne0$ is not easily tractable.