* 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.
preserve.
set printback=none.
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.
release x.
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.