options ps=54 ls=80 nodate nonumber;                                                                                                   
                                                                                                                                        
proc iml;                                                                                                                               
                                                                                                                                        
XY1 = {43.7 19 177,  43.7 43 279, 43.7 56 346, 54.6 13 160,  54.6 19 193,  54.6 43 280, 54.6 56 335,                                    
55.7 13 169, 55.7 26 212, 55.7 34.5 244, 55.7 43 285, 58.8 13 181, 58.8 43 298, 60.5 19 212, 60.5 43 317,                               
60.5 56 347, 61.9 13 186, 61.9 19 216, 61.9 34.5 265, 61.9 43 306, 61.9 56 348, 66.7 13 209,                                            
66.7 43 324,  66.7 56 352};                                                                                                             
                                                                                                                                        
                                                     
                                                                                                                                        
X0 = j(nrow(XY1),1,1);                                       
                                                                                                                                        
X = X0||XY1(|,1|)||XY1(|,2|);                                 
                                                                                                                                        
Y = XY1(|,3|); 

Yhat_i = j(nrow(Y),1,0);

do i=1 to nrow(Y);

if i=1 then do; 
X_i = X(|2:nrow(Y),|);   Y_i = Y(|2:nrow(Y),|);
end;
else if i=nrow(Y) then do;
X_i = X(|1:nrow(Y)-1,|);   Y_i = Y(|1:nrow(Y)-1,|);
end;
else do;
X_i = X(|1:(i-1),|)// X(|(i+1):nrow(Y),|);
Y_i = Y(|1:(i-1),|)// Y(|(i+1):nrow(Y),|);
end;

betahat_i = inv(X_i`*X_i)*X_i`*Y_i;
Yhat_i(|i,|) = X(|i,|)*betahat_i;

end; 

print Y Yhat_i;

PRESS = (Y-Yhat_i)`*(Y-Yhat_i);
MSEP=PRESS/nrow(Y);

Print PRESS MSEP;            

                                                                                                                                        
run;                                                                                                                                    
stop;