/* Computes the eigenvalues and eigenvectors of S and R */
/* In effect, it performs a principal components */
  proc iml;
use _last_;
read all into y;
   n=nrow(y);
   p=ncol(y);
   m=y[ :,];
   d=y-j(n,1)*m;
   s=d`*d/(n-1);
   di=vecdiag(s); di=1/Sqrt(di);dhinv=Diag(di);
   print s;
r=dhinv*s*dhinv; print r;
val=eigval(s); print val; vec= eigvec(s); print vec;
valr=eigval(r); print valr; vecr=eigvec(r); print vecr;
quit; run;

/* Multiple Regression Example  */
data example; input x1 x2 x3 y1 y2 x4;
x1=x1/100; x2=x2/100; x3=x3/100;x4=x4/100; y1=y1/100;y2=y2/100; cards;
  3387  2200  1181  14861  236  6187
  3109  875  3532  11367  310  7336
  2670  957  2502  13329  1182  6988
  3125  1758  4510  12328  1208  6964
  3469  868  3032  12847  1385  8425
  3120  398  2130  13979  1053  6778
  3671  1603  1982  13528  1046  5922
  4531  523  4675  12699  1100  7307
  3678  2034  2354  13534  1349  7679
  3238  1136  4606  11609  1150  8259
  3135  5326  3044  14189  1216  10954
  5271  1658  3340  15052  660  9353
  3728  1945  2111  12236  299  6291
  3506  344  1291  15482  206  4969
  3824  807  1365  14900  239  4825
  3516  1223  1175  15078  161  6019
proc glm; model y1= x1 x2 x3 x4;
proc glm; model y1= x1 x3;
proc reg; model y1= x1 x2 x3 x4; test x2, x4;
/* Multivariate Regression */
prog reg; model y1 y2 = x1 x2 x3 x4;
/* Overall test, the B matrix =0 */
mtest x1, x2, x3, x4;
/* Subtest */
mtest x2, x4;
run;