/*Example 3.6.5 ---Skewness and kurtosis tests*/
/* p=4;n=28; alpha=0.05*/
proc iml;
north={ 72 60 56 41 32 30 39 42 37 33 32 63 54 47
        91 56 79 81 78 46 39 32 60 35 39 50 43 48};
east={  66 53 57 29 32 35 39 43 40 29 30 45 46 51
        79 68 65 80 55 38 35 30 50 37 36 34 37 54};
south={ 76 66 64 36 35 34 31 31 31 27 34 74 60 52
       100 47 70 68 67 37 34 30 67 48 39 37 39 57};
west={  77 63 58 38 36 26 27 25 25 36 28 63 52 45
        75 50 61 58 60 38 37 32 54 39 31 40 50 43};
x=north`||east`||south`||west`;
n=28; p=4; alpha=0.05;
xb=J(1,n,1/n)*x;
xxb=J(n,1,1)*J(1,n,1/n)*x;
S=(x-xxb)`*(x-xxb)/(n-1);
call eigen(Du,H,S);
yb=H`*X`*J(n,1,1/n);
b1p=0; b2p=0;
do i=1 to p;
y1=(H[,i]`*X`-yb[i,]);
y3=y1#y1#y1*J(n,1,1/n);
y4=y1#y1#y1#y1*J(n,1,1/n)/(p*Du[i,]**2);
b1p=b1p+(y3/sqrt(Du[i,])**3)**2/p;
b2p=b2p+y4;
end;
bsk=n*p*b1p/6;
bku=abs(sqrt(n*p/24)*(b2p-3));
chi=cinv(1-0.05, p);
zahalf=probit(1-alpha/2);
pv1=1-probchi(bsk,p);
pv2=1-probnorm(bku);
print bsk chi pv1;
if bsk>chi then print "bskew > chi:   reject";
         else print "bskew < chi:    accept";
print bku zahalf pv2;
If bku>zahalf then print "bkurt > zahalf:    reject";
         else print "bkurt < zahalf:    accept";
quit;
                 BSKEW           CHI            PV1
             7.3105334      9.487729      0.1203601
              bskew < chi:    accept

                 BKURT        ZAHALF            PV2
             0.6925403      1.959964       0.244299
            bkurt < zahalf:    accept