Can we simplify $Cov(XY,X)$? Can we express $Cov(XY,X)$ in terms of moments of $X$ and $Y$ (instead of joint moments)?
Is there an alternative simplification than $$Cov(XY,X)=E[X^2Y]-E[XY]E[X]$$
If not, is there perhaps an approximation for $Cov(XY,X)$?
 A: If you calculate $E(Y|X)=g(X)$ (Since E(Y|X) is a function of $X$)
so
\begin{align}
Cov(XY,X)&=E(X^2Y)-E(XY)E(X)
\\ &= E\color{blue}(E(X^2Y\mid X)\color{blue})-E\color{blue}(E(XY|X)
\color{blue})E(X)
\\ &=E\color{blue}(X^2E(Y\mid X)\color{blue})-E\color{blue}(XE(Y|X)\color{blue})E(X)
\\ &= E(X^2g(X))-E(Xg(X))E(X)
\end{align}
For example consider $$(X,Y)\sim Normal(\mu_x,\mu_y,\sigma_x,\sigma_y,\rho)$$
$E(Y|X)=\mu_y+\rho \frac{\sigma_y}{\sigma_x}(X-\mu_x)$
so
\begin{align}
Cov(XY,X)&= E(X^2g(X))-E(Xg(X))E(X)
\\ &= E(X^2\color{red}(\color{blue}{\mu_y+\rho \frac{\sigma_y}{\sigma_x}(X-\mu_x)}\color{red}))
-E(X\color{red}(\color{blue}{\mu_y+\rho \frac{\sigma_y}{\sigma_x}(X-\mu_x)}\color{red}))E(X)
\end{align}
that can be calculated.
Consider situation that $(X,Y)$ has a complicated distribution but you can calculate $E(Y|X)=g(X)$. if you simulate $x_1,\cdots ,x_N $ from $X$ distribution, by Monte Carlo methods
you can simply approximate
$$Cov(XY,X)=E(X^2g(X))-E(Xg(X))E(X)$$ by
$$\frac{1}{N} \sum_{i=1}^{N} x_i^2 g(x_i)-\color{blue}( \frac{1}{N} \sum_{i=1}^{N} x_i g(x_i) \color{blue})
\color{green}( \frac{1}{N} \sum_{i=1}^{N} x_i \color{green})$$
