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};

/*  XY1 is a matrix containing n=24 rows and 3 columns (X1,X2,Y) we must form the X and Y matrices from it
This is similar to when we input a raw data frame (dataset) and create a matrix  */   

X0 = j(nrow(XY1),1,1);   /* Create a Column of 1s for the intercept (same number of rows as XY1) */
X = X0||XY1(|,1|)||XY1(|,2|);  /* Bind the columns for Intercept (X0), X1 (XY[,1]), X2 (XY[,2]) into X */ 
Y = XY1(|,3|);            /* Create Y from the third column of XY1   */

P = X * inv(X`*X) * X` ;         /* Compute P = X((X'X)^(-1))X   */ 
I_N = i(nrow(X));                  /* Compute I(n)   */
J_N = j(nrow(X),nrow(X),1/nrow(X));   /* Compute (1/n)J  */
                                                                                                                                      
SS_TOT_U = Y` * I_N *Y;    /* Compute SS(Total Uncorrected)  */
DF_TOT_U = sum(diag(I_N)) ; /* Compute df(Total Uncorrected)  */

SS_MODEL  = Y` * P * Y;         /* Compute SS(Model)   */
DF_MODEL = sum(diag(P));         /* Compute df(Model)   */

SS_RESIDUAL  = Y` * (I_N-P) * Y;         /* Compute SS(Residual)   */
DF_RESIDUAL = sum(diag(I_N-P));    /* Compute df(Residual)  */

SS_MEAN  = Y` * J_N * Y;        /* Compute SS(Model)   */
DF_MEAN = sum(diag(J_N));         /* Compute df(Model)  */

SS_REGRESSION = Y` * (P-J_N) * Y;         /* Compute SS(Regression)  */
DF_REGRESSION = sum(diag(P-J_N)); /* Compute df(Regression)   */

print SS_TOT_U DF_TOT_U;
print SS_MODEL DF_MODEL;
print SS_RESIDUAL DF_RESIDUAL;
print SS_MEAN DF_MEAN;
print SS_REGRESSION DF_REGRESSION;

/*  F-test for H0: beta1=beta2=0  (Test Statistic, Critical Value (alpha=0.05,P-value)  */ 
F_OBS = (SS_REGRESSION/DF_REGRESSION)/(SS_RESIDUAL/DF_RESIDUAL);  
F_CRIT = finv(.95,DF_REGRESSION,DF_RESIDUAL);
F_P = 1-probf(F_OBS,DF_REGRESSION,DF_RESIDUAL);

print F_OBS F_CRIT F_P;
run;                                                                                                                                    
stop;