R help - Dec 2009 - How to control the skewness of a heteroscedastic variable? - A Correction

When going through my earlier post I find a mistake in the example 
that I provided. The correct version is provided below. I also start 
to suspect that my problem is that although the cumulant of a sum of 
independent variable is the sum of the cumulants, the moments of a 
sum is not the sum of the moments. But that might not be the only 
flaw in my application.

Regards,
Karl-Oskar 

#An example:

library(moments)
set.seed(1234)

#create two uncorrelated gamma variates
z1<-rgamma(100000,5,scale=sqrt(1/5))
z1<-z1-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 gamma 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)

#The moments seem to differ for the two z's
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)