For two centered (zero expectation) random variables $X$ and $Z$ I am interested in the covariance of the product $XZ$ and either $X$ or $Z$.
$$Cov(X,XZ) = E( X(XZ - E(XZ))) = E(X^2Z)$$
I think the last quantity can be non-zero. However, for any choice of $X$ and $Z$ I have tried, in simulations, I find the quantity to be approx. zero, see code. Perhaps I used examples for which the expectation above happens to be zero or the expectation simplifies further and turns out to be zero. Where do I go wrong?
#With correlated normal variables
library(MASS)
cors = -.9
n = 10000
m = c(0,0)
cova = matrix(c(1,cors,cors,1),ncol=2)
X = mvrnorm(n,m,cova)
cor( cbind(X,X[,1]*X[,2]) )
mean(X[,1]^2*X[,2])
# With independent F-distributed variables
a = rf(n,df1=34,df2=506)
a = a - mean(a)
b = rf(n,df1=34,df2=506)
b = b - mean(b)
cor( cbind(a, b, a*b ))
mean(a^2*b)