options ps=55 ls=76 nodate nonumber;
title 'Ad Expenitures and Retained Impressions';

data one;
input brand $1-15 expend impress;
expend2=expend**2;
expend3=expend**3;
cards;
MILLER LITE	      50.1	32.1
PEPSI                            74.1	99.6
STROH'S	                      19.3	11.7
FED'L EXPRESS           22.9	21.9
BURGER KING              82.4	60.8
COCO-COLA                 40.1	78.6
MC DONALD'S	     185.9	92.4
MCI	                        26.9	50.7
DIET COLA	       20.4	21.4
FORD	                      166.2	40.1
LEVI'S	                        27.0	40.8
BUD LITE	                       45.6  10.4
ATT/BELL	     154.9	88.9
CALVIN KLEIN	         5.0	12.0
WENDY'S	        49.7	29.2
POLAROID	        26.9	38.0
SHASTA	                         5.7	10.0
MEOW MIX	          7.6	12.3
OSCAR MEYER                9.2	23.4
CREST	                         32.4	 71.1
KIBBLES 'N BITS      6.1   4.4
;
symbol1 c=black v=dot;
symbol2 c=black v=dot i=rl l=1;
symbol3 c=black v=dot i=rq l=1;
symbol4 c=black v=dot i=rc l=1;

proc gplot;
plot impress*expend=1 / frame;
title2 'Raw Data';
run;
quit;

proc gplot;
plot impress*expend=2 / frame;
title2 'Linear fit';
run;
quit;

proc gplot;
plot impress*expend=3 /frame;
title2 'Quadratic fit';
run;
quit;

proc gplot;
plot impress*expend=4 /frame;
title2 'Cubic fit';
run;
quit;

proc reg;
title2;
model impress=expend/ss1 ss2;
model impress=expend expend2/ss1 ss2;
model impress=expend expend2 expend3/ss1 ss2;
run;

quit;



data two; 
set one;
keep impress expend expend2 expend3;
run;

proc iml;
use two;
read all;
y=impress;
x=expend||expend2;
intcpt=j(nrow(x),1);
x=intcpt||x;

print x y;

orthx=orpol(expend,2);

print x orthx;

x0=intcpt; x1=expend; x2=expend2;
X01=x0||x1; X02=X0||X2; X12=X1||X2;
PX0=X0*inv(X0`*X0)*X0`;
PX01=X01*inv(X01`*X01)*X01`;
PX02=X02*inv(X02`*X02)*X02`;
PX012=X*inv(X`*X)*X`;

RX1_0=Y`*(PX01-PX0)*Y;
RX2_0=Y`*(PX02-PX0)*Y;
RX2_01=Y`*(PX012-PX01)*Y;
RX1_02=Y`*(PX012-PX02)*Y;

ox0=orthx(|,1|); ox1=orthx(|,2|); ox2=orthx(|,3|);
ox01=ox0||ox1; ox02=ox0||ox2; ox12=ox1||ox2;
Pox0=ox0*inv(ox0`*ox0)*ox0`;
Pox01=ox01*inv(ox01`*ox01)*ox01`;
Pox02=ox02*inv(ox02`*ox02)*ox02`;
Pox012=orthx*inv(orthx`*orthx)*orthx`;

Rox1_0=Y`*(Pox01-Pox0)*Y;
Rox2_0=Y`*(Pox02-Pox0)*Y;
Rox2_01=Y`*(Pox012-Pox01)*Y;
Rox1_02=Y`*(Pox012-Pox02)*Y;

print rx1_0 rx1_02 rx2_0 rx2_01;

print rox1_0 rox1_02 rox2_0 rox2_01;

run;

stop;