options nodate nonumber ps=55 ls=80;

title 'RPD -- Examples 4.2-4.5 -- 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;

SSTOTALU = Y`*Y;   /*  Alternative:  Y[##]  */

SSMODEL = BETAHAT`*XTY;

SSRESID = SSTOTALU - SSMODEL;

SSMU = Y`*((1/4)#j(4))*Y;   /* j(m) = m-by-m matrix of 1's */

   /*  More efficient alternative:  4 * Y[:]**2  */

SSTOTALC = SSTOTALU - SSMU;

SSREG = SSMODEL - SSMU;

print SSTOTALU SSMU SSTOTALC SSREG SSRESID;


S2 = SSRESID/2;

S2BETAHAT = S2#XTXINV;

print S2 S2BETAHAT;