[R] How to control the skewness of a heteroscedastic variable?

[Karl-Oskar Lindgren]

Dear listusers,

I don't know whether my problem is statistical or computational, but 
I hope I could recieve some help in either case. 

I'm currently working on a MC-simulation in which I would like to 
control the skewness of a heteroscedastic dependent variable defined 
as:

y=d*z+sqrt(.5+.5*x^2)*e 					(eq.1)

where d is a parameter and, z, x, and e are gamma r.vs. The variables 
x (the one creating the heteroscedasticity) and z are assumed to be 
positively correlated. 

I thought that since the two terms on the rhs of eq.1 are 
uncorrelated the 3rd central moment of y should equal the sum of the 
3rd central moments of the two terms on the rhs. This seems to be 
correct as long as x and z in eq1 are uncorrelated. But if I make x 
and z correlated the 3rd moment of y exceeds the sum of the 3rd 
moments of the terms on the rhs. 

My problem is that I cannot understand why this is the case (there 
seems to be no linear correlation between z and the multiplicative 
error term). Is it my statistical understanding or my computational 
set-up that is flawed? Basically what I want to do is to control the 
skew of y in eq.1 in my simulations by varying the skew of e. Is that 
possible to do, and if so, how would that best be implemented in R?

The code below provides an illustration of my problem in case my 
verbal explanation was difficult to follow.

Regards,
Karl-Oskar Lindgren
Researcher 
Department of Government
Uppsala University
 
#An example:

library(moments)
set.seed(1234)

#create two uncorrelated gamma variates
z1<-rgamma(100000,5,scale=sqrt(1/5))
z1<-z-5*sqrt(1/5)
x1<-rgamma(100000,10,scale=sqrt(1/10))
x1<-x1-10*sqrt(1/10)

#create two correlated gamma variates
R<-matrix(c(1,.5,.5,1),2,2)
Y<-cbind(x1,z1)%*%chol(R)
x2<-Y[,1]
z2<-Y[,2]

#create a gamma dist. error term
e<-rgamma(100000,2,scale=sqrt(1/2))
e<-e-2*sqrt(1/2)

#create the heteroscedasticity functions
h1<-sqrt(.5+.5*x1^2)
h2<-sqrt(.5+.5*x2^2)

#create the heteroscedastic dependent variables
y1<-.5*z1+h1*e
y2<-.5*z2+h2*e

#The 3rd moments of y1 and y2 differ
moment(y1,3,central=T)
moment(y2,3,central=T)

#But the moments of terms on rhs seem to be the same
moment(.5*z1,3,central=T)
moment(.5*z2,3,central=T)
moment(h1*e,3,central=T)
moment(h2*e,3,central=T)

#The corr bw z and the het. error terms
#seems to be the same in the two cases
var(h1*e,.5*z1)
var(h2*e,.5*z2)