options nodate nonumber ps=55 ls=80;

title 'RPD -- Examples 3.1+ -- X = Ozone, Y = Yield';

proc iml;

X = {1 0.02, 1 0.07, 1 0.11, 1 0.15};

Y = {242, 237, 231, 201};

XTX = X`*X;

XTY = X`*Y;

XTXINV = inv(XTX);

BETAHAT = XTXINV*XTY;

P = X*XTXINV*X`;

YHAT = P*Y;

E = (i(4)-P)*Y;

print X Y XTX XTY XTXINV;

print BETAHAT P YHAT E;


/* Much more computationally efficient ways to obtain BETAHAT, YHAT, and E: */

BETAHAT = solve(XTX,XTY);   /* X = solve(A,B) solves A*X = B */

YHAT = X*BETAHAT;

E = Y - YHAT;

print BETAHAT YHAT E;