proc iml;
start main;

y={5,
   20,
   6,
   10};

  X={1,
     1,
     1,
     1};

A={1  0  0  0 .5 .5  0  0,
   0  1  0  0 .5 .5  0  0,
   0  0  1  0  0  0 .5 .5,
   0  0  0  1  0  0 .5 .5,
  .5 .5  0  0  1 .5  0  0,
  .5 .5  0  0 .5  1  0  0,
   0  0 .5 .5  0  0  1 .5,
   0  0 .5 .5  0  0 .5  1};

AINV=INV(A);


Z1={ 0  0  0  0  1  0  0  0,
	 0  0  0  0  0  1  0  0,
     0  0  0  0  0  0  1  0,
     0  0  0  0  0  0  0  1};

Z2={  0  0  0  0  0  0  1  0,
      0  0  0  0  0  0  0  1,
	  0  0  0  0  1  0  0  0,
      0  0  0  0  0  1  0  0};

P={50 -5,
   -5 10};

K=inv(P)#50;

LHS=((X`*X)||(X`*Z1)||(X`*Z2))
//((Z1`*X)||(Z1`*Z1+AINV#K[1,1])||(Z1`*Z2+AINV#K[1,2]))
//((Z2`*X)||(Z2`*Z1+AINV#K[2,1])||(Z2`*Z2+AINV#K[2,2]));

RHS=(X`*Y)//(Z1`*Y)//(Z2`*Y);
C=INV(LHS);
BU=C*RHS;
RMSE=(Y`*Y-BU`*RHS)#(1/3);
print BU RMSE;
finish main;
run;
quit;