options nodate nonumber ps=55 ls=100;
title 'Retail Sales Regression Analysis';
title2 "U.S. Counties in SMSA's";

data retail;
infile 'C:\SAS\sta6208\smsa\retail.txt';
input county $ 22-55 estab 63-67 sales 71-77;
run;

proc sort; by county;
run;

data income;
infile 'C:\SAS\sta6208\smsa\income.txt';
input county $ 22-55 income 59-67;
run;

proc sort; by county;
run;

data pop;
infile 'C:\SAS\sta6208\smsa\pop.txt';
input county $ 22-55 pop mper100f;
run;

proc sort; by county;
run;

data fedexp;
infile 'C:\SAS\sta6208\smsa\fedexp.txt';
input county $ 22-55 fedexp 60-67;
run;

proc sort; by county;
run;

data all;
merge retail income pop fedexp;
by county;
pcsales=100*sales/pop;
pcestab=1000*estab/pop;
pcinc=10*income/pop;
pcfedexp=10*fedexp/pop;
mper100f=mper100f/10;
run;

proc sort; by pcsales;
run;


proc print;
var county pcsales pcestab pcinc pcfedexp mper100f;
run;



proc corr;
var pcsales pcestab pcinc pcfedexp mper100f;
run;

quit;



proc reg;
model pcsales = pcestab pcinc pcfedexp mper100f / ss1 ss2;
run;

quit;



proc iml;
use all;
read all;
Y=pcsales;
X0=j(845,1,1);
X1=pcestab;
X2=pcinc;
X3=pcfedexp;
X4=mper100f;

X01234=X0||X1||X2||X3||X4; P01234=X01234*inv(X01234`*X01234)*X01234`;
X0123=X0||X1||X2||X3; P0123=X0123*inv(X0123`*X0123)*X0123`;
X0124=X0||X1||X2||X4; P0124=X0124*inv(X0124`*X0124)*X0124`;
X0134=X0||X1||X3||X4; P0134=X0134*inv(X0134`*X0134)*X0134`;
X0234=X0||X2||X3||X4; P0234=X0234*inv(X0234`*X0234)*X0234`;
X1234=X1||X2||X3||X4; P1234=X1234*inv(X1234`*X1234)*X1234`;
X012=X0||X1||X2; P012=X012*inv(X012`*X012)*X012`;
X013=X0||X1||X3; P013=X013*inv(X013`*X013)*X013`;
X014=X0||X1||X4; P014=X014*inv(X014`*X014)*X014`;
X023=X0||X2||X3; P023=X023*inv(X023`*X023)*X023`;
X024=X0||X2||X4; P024=X024*inv(X024`*X024)*X024`;
X034=X0||X3||X4; P034=X034*inv(X034`*X034)*X034`;
X01=X0||X1; P01=X01*inv(X01`*X01)*X01`;
X02=X0||X2; P02=X02*inv(X02`*X02)*X02`;
X03=X0||X3; P03=X03*inv(X03`*X03)*X03`;
X04=X0||X4; P04=X04*inv(X04`*X04)*X04`;
X0=X0; P0=X0*inv(X0`*X0)*X0`;

i845=i(845);

r01234=y`*p01234*y;
r0123=y`*p0123*y;
r0124=y`*p0124*y;
r0134=y`*p0134*y;
r0234=y`*p0234*y;
r1234=y`*p1234*y;
r012=y`*p012*y;
r013=y`*p013*y;
r014=y`*p014*y;
r023=y`*p023*y;
r024=y`*p024*y;
r034=y`*p034*y;
r01=y`*p01*y;
r02=y`*p02*y;
r03=y`*p03*y;
r04=y`*p04*y;
r0=y`*p0*y;

r0_1234=r01234-r1234;
r1_0=r01-r0; r1_0234=r01234-r0234;
r2_01=r012-r01; r2_0134=r01234-r0134;
r3_012=r0123-r012; r3_0124=r01234-r0124;
r4_0123=r01234-r0123;




ssreg=y`*(p01234-p0)*y; 
dfreg=round(trace(((X01234`*X01234)*inv(X01234`*X01234)))-trace(((X0`*X0)*inv(X0`*X0))));
msreg=ssreg/dfreg;

ssres=y`*(i845-p01234)*y;
dfres=round(trace(i845)-trace(((X01234`*X01234)*inv(X01234`*X01234))));
msres=ssres/dfres;


f_obs=msreg/msres;
pval_f=1-probf(f_obs,dfreg,dfres);

sstotc=y`*(i845-p0)*y;
dftotc=round(trace(i845)-trace(((X0`*X0)*inv(X0`*X0))));

s=sqrt(msres);

r_square=ssreg/sstotc;

ybar=j(1,845,1/845)*Y;

Adj_rsq=1-((ssres/dfres)/(sstotc/dftotc));

CV=100*s/ybar;

betahat=inv(X01234`*X01234)*X01234`*Y;

s2_betahat=msres*inv(X01234`*X01234);

betahat0=betahat(|1,|); sbetahat0=sqrt(s2_betahat(|1,1|));
t0=betahat0/sbetahat0; p0=2*(1-probt(abs(t0),dfres));

betahat1=betahat(|2,|); sbetahat1=sqrt(s2_betahat(|2,2|));
t1=betahat1/sbetahat1; p1=2*(1-probt(abs(t1),dfres));

betahat2=betahat(|3,|); sbetahat2=sqrt(s2_betahat(|3,3|));
t2=betahat2/sbetahat2; p2=2*(1-probt(abs(t2),dfres));

betahat3=betahat(|4,|); sbetahat3=sqrt(s2_betahat(|4,4|));
t3=betahat3/sbetahat3; p3=2*(1-probt(abs(t3),dfres));

betahat4=betahat(|5,|); sbetahat4=sqrt(s2_betahat(|5,5|));
t4=betahat4/sbetahat4; p4=2*(1-probt(abs(t4),dfres));

print dfreg ssreg msreg f_obs pval_f;
print dfres ssres msres;
print dftotc sstotc;

print s r_square;

print ybar adj_rsq;

print cv;

print betahat0 sbetahat0 t0 p0 r0 r0_1234;
print betahat1 sbetahat1 t1 p1 r1_0 r1_0234;
print betahat2 sbetahat2 t2 p2 r2_01 r2_0134; 
print betahat3 sbetahat3 t3 p3 r3_012 r3_0124;
print betahat4 sbetahat4 t4 p4 r4_0123 r4_0123;

run;

stop iml;


quit;