preserve. set printback=none. * Mardia's multivariate skew (b1p) and multivariate kurtosis (b2p) * Author: Lawrence T. DeCarlo, 11/97 * Email: decarlo@exchange.tc.columbia.edu * * Multivariate skew is provided in a separate macro because it * is more computationally intense, particularly for large datasets * * Note: increase mxloops if n>50000. define mardia (vars=!charend('/')). set mxloops=50000. matrix. get x /variables=!vars /names=varnames /missing=omit. compute n=nrow(x). compute p=ncol(x). compute xbar=csum(x)/n. compute j=make(n,1,1). compute xdev=x-j*xbar. compute dev=csum(xdev). compute devsq=(dev&*dev)/n. compute ss=csum(xdev&*xdev). compute m2=(ss-devsq)/n. compute sdev=sqrt(m2). compute m3=csum(xdev&**3)/n. compute m4=csum(xdev&**4)/n. compute sqrtb1=m3/(m2&*sdev). compute b2=m4/(m2&**2). compute g1=((sqrt(n*(n-1)))*sqrtb1)/(n-2). compute g2=(b2-((3*(n-1))/(n+1)))*((n**2-1)/((n-2)*(n-3))). compute seskew=sqrt(6*n*(n-1) / ((n-2)*(n+1)*(n+3))). compute sekurt=sqrt(4*((n&**2)-1)*(seskew&**2) / ((n-3)*(n+5))). release x. print {n}/title"Number of observations:" /format=f5. print {p}/title"Number of variables:" /format=f5. print {g1,seskew}. /title"Univariate Skewness" /clabels=!vars,"Std.Err"/format=f10.4. print {g2,sekurt}. /title"Univariate Kurtosis" /clabels=!vars, "Std.Err" /format=f10.4. compute s=sscp(xdev)/n. compute sinv=inv(s). compute gii=make(n,1,0). compute gij=make(n,1,0). compute gsum=make(n,1,0). loop i=1 to n. + compute gii(i)=xdev(i,:)*sinv*t(xdev(i,:)). + loop j=1 to n. + compute gij(j)=xdev(i,:)*sinv*t(xdev(j,:)). + end loop. + compute gsum(i)=csum(gij&**3). end loop. compute b1p=csum(gsum)/(n*n). compute chib1p=(n*b1p)/6. compute sm=((p+1)*(n+1)*(n+3))/(n*((n+1)*(p+1)-6)). compute chism=(n*b1p*sm)/6. compute df=(p*(p+1)*(p+2))/6. compute pb1p=1-chicdf(chib1p,df). compute pb1psm=1-chicdf(chism,df). print {b1p,chib1p,pb1p,chism,pb1psm} /title"Mardia's multivariate skew (small sample adjustment: Mardia 1974 Sankya)". /clabels="b1p","Chi(b1p)","p-value","adj-Chi","p-value"/format=f10.4. compute b2p=csum(gii&**2)/n. compute nb2p=(b2p-p*(p+2))/sqrt(8*p*(p+2)/n). compute pnb2p=2*(1-cdfnorm(abs(nb2p))). print {b2p,nb2p,pnb2p}. /title"Mardia's multivariate kurtosis" /clabels="b2p","N(b2p)","p-value"/format=f10.4. end matrix. !enddefine. restore.