#Run the spss SkewKurtosis helper function one time, and load it into memory
#Code underneath the helper function shows how to obtain skewness/kurtosis + 
#standard error and z-statistic for one variable


spssSkewKurtosis=function(x) {
  w=length(x)
  m1=mean(x)
  m2=sum((x-m1)^2)
  m3=sum((x-m1)^3)
  m4=sum((x-m1)^4)
  s1=sd(x)
  skew=w*m3/(w-1)/(w-2)/s1^3
  sdskew=sqrt( 6*w*(w-1) / ((w-2)*(w+1)*(w+3)) )
  kurtosis=(w*(w+1)*m4 - 3*m2^2*(w-1)) / ((w-1)*(w-2)*(w-3)*s1^4)
  sdkurtosis=sqrt( 4*(w^2-1) * sdskew^2 / ((w-3)*(w+5)) )
  mat=matrix(c(skew,kurtosis, sdskew,sdkurtosis), 2,
             dimnames=list(c("skew","kurtosis"), c("estimate","se")))
  return(mat)
}

skewkurt<-as.data.frame(spssSkewKurtosis(Wk04$Income))

skewkurt$z.ratio<-skewkurt$estimate / skewkurt$se

skewkurt

##           estimate        se   z.ratio
## skew      4.575384 0.1199527  38.14323
## kurtosis 25.195441 0.2393378 105.27147