Si la hipótesis crucial (que falta en tu post) que $(X,Y)$ es conjuntamente normal se mantiene, un enfoque simple es realizar $X$ y $Y$ a través de variables aleatorias normales i.i.d. A saber, $$ X=\mu_X+\sigma_X\xi,\qquad Y=\mu_Y+\sigma_Y\cos(t)\xi+\sigma_Y\sin(t)\eta, $$ donde $(\xi,\eta)$ son variables aleatorias normales i.i.d. y el ángulo $t$ es tal que $\sigma_X\sigma_Y\cos(t)=C_{X,Y}$ donde $C_{X,Y}=\mathrm{Cov}(X,Y)$ .
( Prueba: Definir las variables aleatorias $\xi$ y $\eta$ a través de las identidades anteriores y comprueba que $\mathrm E(\xi^2)=\mathrm E(\eta^2)=1$ y $\mathrm E(\xi)=\mathrm E(\eta)=\mathrm E(\xi\eta)=0$ . Fin de la prueba).
Entonces, se puede expresar todo en términos de momentos de $\xi$ y $\eta$ sólo. Por ejemplo, $$ X^2Y=(\mu_X^2+2\mu_X\sigma_X\xi+\sigma_X^2\xi^2)\cdot(\mu_Y+\sigma_Y\cos(t)\xi+\sigma_Y\sin(t)\eta), $$ y $\mathrm E(\xi)=\mathrm E(\eta)=\mathrm E(\xi\eta)=\mathrm E(\xi^3)=0$ mientras que $\mathrm E(\xi^2)=\mathrm E(\eta^2)=1$ . El desarrollo del producto y la identificación de la expectativa de cada término produce $$ \mathrm E(X^2Y)=\mu_X^2\mu_Y+2\mu_X\sigma_X\sigma_Y\cos(t)+\sigma_X^2\mu_Y, $$ es decir, $$ \color{red}{\mathrm E(X^2Y)=\mu_X^2\mu_Y+2\mu_XC_{X,Y}+\sigma_X^2\mu_Y}. $$ Comprobación de cordura: Si $\mathrm{Cov}(X,Y)=0$ el resultado es $\mathrm E(X^2)\mathrm E(Y)$ como debe ser.
Es posible que quieras calcular $\mathrm E(X^3Y)$ utilizando la misma técnica, un ingrediente no utilizado anteriormente es el valor numérico $\mathrm E(\xi^4)=3$ .
Editar: Un enfoque similar, para calcular cada $\mathrm E(X^nY^k)$ al mismo tiempo, es considerar la función $\varphi$ definido por $$ \varphi(u,v)=\mathrm E(\mathrm e^{uX+vY})=\sum\limits_{n\geqslant0}\sum\limits_{k\geqslant0}\frac{u^n}{n!}\frac{v^k}{k!}\mathrm E(X^nY^k). $$ Dado que la variable aleatoria $uX+vY$ es normal con media y varianza $$ m(u,v)=u\mu_X+v\mu_Y,\qquad s(u,v)^2=u^2\sigma_X^2+2uvC_{X,Y}+v^2\sigma_Y^2, $$ uno sabe que $\varphi(u,v)=\exp(m(u,v)+\frac12s(u,v)^2)$ . La tarea ahora es desarrollar esto como una serie de potencias en $(u,v)$ y para identificar el coeficiente de $u^nv^k$ . Esto nos lleva a la fórmula $$ \mathrm E(X^nY^k)=n!\,k!\,\sum_\ast\frac{\mu_X^i}{i!}\frac{\mu_Y^j}{j!}\frac{\sigma_X^{2r}}{2^rr!}\frac{\sigma_X^{2s}}{2^ss!}\frac{C_{X,Y}^t}{t!}, $$ donde la suma $\sum\limits_\ast$ es sobre cada no negativo $(i,j,r,s,t)$ tal que $$ i+2r+t=n,\qquad j+2s+t=k. $$ Para $(n,k)=(2,1)$ los posibles valores de $(i,j,r,s,t)$ son $$ (2,1,0,0,0),\quad (0,1,1,0,0),\quad (1,0,0,0,1), $$ que corresponden respectivamente a los monomios $$\mu_X^2\mu_Y,\quad \mu_Y\sigma_X^2,\quad \mu_XC_{X,Y}. $$