***------------------------------------------------------------***;
**input Ryan and Joiner (1994) data;
data ryan;
infile 'd:\kurt\ryan\ryan.dat';
input v1-v3;
run;
***--------------------------------------------------------------------***;
**Example for Mardia's kurtosis;
** based on marginal variables when missing variables;
***--------------------------------------------------------------------***;
options ls=75 nodate nonumber;
options nodate; title;

proc iml;
reset noname;
***--------------------------------------------------------------------***;
**** The following subroutine is to perform vech(.) function;
***--------------------------------------------------------------------***;
start vech(A,Va);
l=0;
p=nrow(A);
pstar=p*(p+1)/2;
Va=j(pstar,1,0);
do i=1 to p by 1;
        do j=i to p by 1;
                l=l+1;
                Va(|l|)=A(|j,i|);
        end;
end;
finish;
***--------------------------------------------------------------------***;
**Subroutines for impute missing values;
***--------------------------------------------------------------------***;

** E algorithm for unstructured (\mu,sigma) starts here;
start emmus1(xbar0,sbig0,xmis, sumx,sumxx);
*impute 3rd variable, p=3;
n=nrow(xmis); p=ncol(xmis); mu=xbar0; sig=sbig0;
sumx=j(p,1,0); sumxx=j(p,p,0);

do i= 1 to n;
        xi=xmis(|i,|)`;
        xio=xi(|1:2|);
        muo=mu(|1:2|);
        mum=mu(|3|);
        sigoo=sig(|1:2,1:2|);
        sigom=sig(|1:2,3|);
        sigmo=sigom`;
        sigmm=sig(|3,3|);
        sigooin=inv(sigoo);
        delti=sigmm-sigmo*sigooin*sigom;
        delt=j(p,p,0);
        delt(|3,3|)=delti;
        xm1=mum+sigmo*sigooin*(xio-muo);
        xi(|3|)=xm1;
        xxi=xi*xi`+delt;
        sumx=sumx+xi;
        sumxx=sumxx+xxi;
end;
finish;
***--------------------------------------------------------------------***;
** Computing Mardia's multivariate kurtosis;
***--------------------------------------------------------------------***;

start Mdkurt(x, mu1, sig1, mkurt);
nsamp=nrow(x);
n=nsamp;
p=ncol(x);

mkurt=0.0;
mean=mu1;
sigin=inv(sig1);

do i=1 to nsamp;
        xi=x(|i,|)`;
        xi0=xi-mean; *centralized the ith case;
        di2=xi0`*sigin*xi0; **Mahalonobis-distance;
        mkurt=mkurt+di2*di2;
end;

finish;


***--------------------------------------------------------------------***;
** creat missing data pattern, estimate mu & sigma, and kurtosis;
***--------------------------------------------------------------------***;
start misptn(x, mkurt0);
ep=.000000001;
n=nrow(x);
p=ncol(x);

***------------------------------------------------------------***;
nm=9;
nc=n-nm;
xc=x(|1:nc,|);
xm1=x(|nc+1:n,|);
***------------------------------------------------------------***;

*print "complete case=" nc;
pmi=2; **# of observed marginals;

eap2=nc*p*(p+2)+nm*pmi*(pmi+2);
varap2=8*eap2;

sumxc=j(p,1,0);
sumxxc=j(p,p,0);
do i=1 to nc;
        xi=xc(|i,|)`;
        sumxc=sumxc+xi;
        sumxxc=sumxxc+xi*xi`;
end;
xbar0=j(p,1,0);
sbig0=i(p);

** using EM to get better estimates (xbar1,sbig1);
label1:
run emmus1(xbar0,sbig0,xm1,sumxm,sumxxm);
*print "3rd variable is missing";
xbar1=(sumxc+sumxm)/n;
sbig1=(sumxxc+sumxxm-n*xbar1*xbar1`)/n;
dt=(ssq(xbar1-xbar0)+ssq(sbig1-sbig0));
**/(ssq(xbar1)+ssq(sbig1));
if dt>ep then do;
        xbar0=xbar1;
        sbig0=sbig1;
        goto label1;
end;

print "MLE of mean=" xbar1;
print "MLE of cov=" sbig1;

run Mdkurt(xc, xbar1, sbig1, mkurtc);

xm1=xm1(|,1:2|);
xbarm=xbar1(|1:2|); sbigm=sbig1(|1:2,1:2|);
run Mdkurt(xm1, xbarm, sbigm, mkurtm);

mkurt=mkurtc+mkurtm;
mkurt0=(mkurt-eap2)/sqrt(varap2);

finish;


***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;
print "--------------------------------------------------------------------";
print "Ryan and Joiner (1994) data";
print "---------------------------------------------------------";
use ryan;
read all var _num_ into xmatr;
x=xmatr;
n=nrow(x);
p=ncol(x);
nc=19;

print "Sample size n=" n;

**creat missing data;
run misptn(x, mkurt0);

print "mkurt=" mkurt0;

finish;
run main;