DATA SQUARE;                                                                   
INPUT y x1 x2;                                    
CARDS;                                                                          
457.7 17 27
568.2 22 27
452.8 16 24
783.7 34 23
629.4 30 21
216.9 13 17
1014.2 31 32
471.6 25 19
266.4 32 11
243.3 15 11
278.0 18 15
354.6 12 27
534.2 32 17
374.2 26 16
441.1 18 21
762.7 31 22
702.8 23 33
582.6 30 19
590.4 34 19
131.6 16 13
415.0 33 11
606.3 28 21
628.7 20 31
130.3 14 11
359.6 34 12
438.2 14 29
417.9 27 15
332.7 21 17
626.6 29 33
575.2 16 34
641.9 27 24
448.5 26 15
536.9 21 24
742.2 30 25
539.6 19 29
516.8 19 25
791.2 27 29
516.9 15 31
493.8 30 18
366.3 32 12
634.6 24 23
194.4 13 19
734.5 28 26
432.9 26 16
362.7 15 26
231.4 11 16
442.4 27 17
538.5 22 28
365.6 29 10
315.9 17 19
;
run; 
proc iml; use square; read all  into x;  
                                        
print x;         
 
n=nrow(x);kcol=ncol(x);                                                         
y=x[,1];z=j(n,1,1);x[,1]=z;m=n;                                                                                                                 
ADJ=I(m)-X*(INV(X`*X))*X`;                                                      
E=ADJ*Y;                                                                        
D=J(M,M,0);                                                                     
adj2=i(m)-(z*z`)/m;do kk=1 to kcol;                                          
ss=adj2*x[,kk];sd=ss`*ss/m;sd1=(sd)##(.5);x[,kk]=x[,kk]/sd1;                  
end;                                                            
DO K=1 TO M; DO L=1 TO M;do kk=1 to kcol;                                       
d[k,l]=d[k,l]+(X[K,kk]-X[L,kk])**2;end;                                         
IF( K=L | D[k,l]=0 ) THEN  D[K,L]=0;ELSE D[K,L]=1/D[K,L];                                    
END;END;                                                                        
B=ADJ*D*ADJ;                                                                    
NUM=E`*D*E;                                                                     
DEN=E`*E;                                                                       
DEN=DEN/(M-kcol);                                                              
STAT=NUM/DEN;                                                                  
MEAN=TRACE(B);                                                                  
VAR=2*((M-3)*TRACE(B**2)-(TRACE(B)**2) )/(M-1);                                
STAT=(STAT-MEAN)/SQRT(VAR);                                                    
PRINT STAT;