PROC IMPORT OUT= WORK.one 
            DATAFILE= "C:\DATA\SAS\clothing.csv" 
            DBMS=CSV REPLACE;
	 /* tsales,sales,margin,nown,nfull,npart,naux,hoursw,hourspw,inv1,inv2,ssize,start */
     GETNAMES=YES;
     DATAROW=2; 
     GUESSINGROWS=20;
run;
data two;
set one;
int=1;
proc reg data = two;
var obs;
model sales = hoursw ssize ;
model sales = ssize nfull npart ;
model sales = hoursw ssize nfull npart/ selection=adjrsq aic sbc;
run;
proc autoreg data = two;
model sales = hoursw ssize;
model sales = hoursw ssize nown;
test nown;
model sales = hoursw ssize nown npart;
test npart;
run;
/* Reset for model A, Non nested F test */
proc autoreg;
model sales = hoursw ssize /reset ;
model sales = hoursw ssize nfull npart;
test nfull,npart;
model sales = ssize nfull npart hoursw;
test hoursw;
run;

proc iml;
use two;
read all var{sales} into sales;
read all var{int hoursw ssize } into X1 ;
read all var{int ssize nfull npart} into X2;
read all var{int hoursw ssize nfull npart} into X3;
read all var{nown} into nown;
read all var{npart} into npart;

/*  Wald Test */
start Wald(y,x,R,lr);
  T=nrow(X);
  K=ncol(X);
  J=nrow(R);
  t_k=t-k;
  b=inv(x`*x)*x`*y;
  sig2=(y-x*b)`*(y-x*b)/(T_K);
  W=(R*b-lr)`*inv(R*inv(x`*x)*R`)*(R*b-lr)/(J*sig2);
  pval=1-probf(w,j,t_k);
  print w J T_K pval;
finish;

/* The RESET Test */
start RESET(y,x);
  k=ncol(x);
  b=inv(x`*x)*x`*y;
  yhat=x*b;
  e2=yhat#yhat;
  e3=yhat#yhat#yhat;
  z1=x||e2;
    r=j(1,k,0)||1;
    lr=0;
    Print "Reset Power 2";
    run wald(y,z1,r,lr);
  z1=x||e2||e3;
    r=j(2,k,0)||i(2);
    lr=0;
    Print "Reset Power 3";
    run wald(y,z1,r,lr);
finish;

/* R-square */
start rsquare(x,y);
  n=nrow(y);
  k=ncol(x);
  b=inv(x`*x)*x`*y;
  yhat=x*b;
  ybar=sum(y)/n;
  SSR=ssq(yhat-ybar);
  SST=ssq(y-ybar);
  SSE=ssq(y-x*b);
  r2=ssr/sst;
  rbar2 = 1-((n-1)/(n-k))*(1-r2);
  print SSR SST SSE r2 rbar2;
finish;

start modelsel(x,y);
  n=nrow(X);
  K=ncol(X);
  b=inv(x`*x)*x`*y;
  sse = (y-x*b)`*(y-x*b);
  AIC = n*log(sse/n) + 2*k/1;
  pc = (sse/(n-k))*(1 + k/n);
  sc = n*log(sse/n) + k*log(n)/1;
  print aic pc sc;
  /* Greene's versions of sc and aic are obtained by dividing by n */
finish;

start ols(x,y,vname);
  T=nrow(X);
  K=ncol(X);
  b=inv(x`*x)*x`*y;

  /*  Sigma squared */
  sse = (y-x*b)`*(y-x*b);
  sig2 = sse /(T-K);
  print sse sig2;

  /* Estimated Covariance Matrix */
  cov=sig2*inv(x`*x);
  *print "The variance covariance matrix", cov;

  /* Standard Errors */
  se=sqrt(vecdiag(cov));
  /* Tvalues */
  tvals=b/se;
  pvals=2*(1-probt(abs(tvals),t-k));

  print "Ordinary Least Squares";
  print vname b se tvals pvals;
finish;

vnameA = {int hoursw ssize }`;
Print "Model A";
run ols(x1,sales,vnameA); 
run reset(sales,x1);
vnameB = {int ssize nfull npart}`;
Print "Model B";
run ols(x2,sales,vnameB); 
vnameC = {int hoursw ssize nfull npart}`;
Print "Model C";
run ols(x3,sales,vnameC);
run reset(sales,x3); 

Print "Nonnested F-Test";
Bigr={0 0 0 1 0,0 0 0 0 1};
  lr={0,0};
  run wald(sales,x3,bigr,lr);
BigR={0 1 0 0 0};
  lr=0;
  run wald(sales,x3,bigr,lr);

/* J test */

px1=x1*inv(x1`*x1)*x1`;
px2=x2*inv(x2`*x2)*x2`;
yhat1=px1*sales;
yhat2=px2*sales;

Xj1=x1||yhat2;
xj2=x2||yhat1;
Print "J-test results";
vnameAj = {int hoursw ssize yhat}`;
run ols(xj1,sales,vnameAj);
vnameBj = {int ssize nfull npart yhat}`;
run ols(xj2,sales,vnameBj);

print "Model Selection Results";
run rsquare(x1,sales);
run rsquare(x2,sales);
run rsquare(x3,sales);
run modelsel(x1,sales);
run modelsel(x2,sales);
run modelsel(x3,sales);
run;