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.