/****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: STATDEMO                                            */
 /*   TITLE: STATISTICAL EXAMPLE                                 */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: MATRIX                                              */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: RHD                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

proc iml;
reset print;

x = {1 1 1 , 1 2 4 , 1 3 9 , 1 4 16, 1 5 25};
y = {1,5,9,23,36};
b = inv(x`*x)*x`*y;
yhat = x*b;
r = y-yhat;
sse =ssq(r);

start reg(x,y);
  xpxi=inv(x`*x);
  beta = xpxi*(x`*y);
  yhat = x*beta;
  resid = y-yhat;
  sse = ssq(resid);
  n = nrow(x);
  dfe = n-ncol(x);
  mse = sse/dfe;
  cssy = ssq(y-y[+]/n);
  rsquare = (cssy-sse)/cssy;
  print , "Regression Results" , , sse dfe mse rsquare;
  stdb = sqrt(vecdiag(xpxi)*mse);
  tratio = beta/stdb;
  prob = 1-probf(tratio#tratio,1,dfe);
  print , "Parameter Estimates" , , beta stdb tratio prob;
  print , y yhat resid;
 finish;

reset noprint;
run reg(x,y);