PROC IMPORT OUT= WORK.one 
            DATAFILE= "C:\DATA\SAS\wage1.csv" 
            DBMS=CSV REPLACE;
	 /* wage	educ	exper	expersq	female	married	lwage	nonwhite */
     GETNAMES=YES;
     DATAROW=2; 
run;
data two;
set one;
int=1;
proc model data = two;
parms b1-b5;
var lwage educ exper expersq female;
lwage =b1+b2*educ+b3*exper+b4*expersq+b5*female;
fit ;
test b2,b3,b4,b5;
fit lwage/hccme=0;
fit lwage/hccme=3;
test b2,b3,b4,b5;
run;

proc iml;
use two;
read all var{lwage} into lwage;
read all var{int educ exper expersq female } into X1 ;

/*  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;

/*  Wald Test */
start Wald2(y,x,R,lr,cov);
  J=nrow(R);
  b=inv(x`*x)*x`*y;
  W=(R*b-lr)`*inv(R*cov*R`)*(R*b-lr);
  pval=1-probchi(w,j);
  print w J pval;
finish;

/*  Lambda 2 form of F-test */
start lambda(y,x,bigr,lr);
  T=nrow(X);
  K=ncol(X);
  J=nrow(bigr);
  b=inv(x`*x)*x`*y;
  sseu=(y-x*b)`*(y-x*b);
  sig2=(sseu)/(T-K);
  RLS=b+inv(x`*x)*BigR`*inv(BigR*inv(x`*x)*Bigr`)*(Lr-BigR*b);
  sser=(y-x*rls)`*(y-x*rls);
  lambda2=(sser-sseu)/(J*sig2);
  pval2=probf(lambda2,J,T-K);

  Print "Other forms of the test statistics";
  print lambda2 pval2;

finish;

/* HCCME */
start HCCME0(y,x);
  k=ncol(x);
  invxtx=inv(x`*x);
  b=invxtx*x`*y;
  e=y-x*b;
  e2=e#e;
  cov=invxtx*x`*diag(e2)*x*invxtx;
  se=sqrt(vecdiag(cov));
  print b se;
  return(cov);
finish;

start HCCME3(y,x);
  k=ncol(x);
  invxtx=inv(x`*x);
  b=invxtx*x`*y;
  e=y-x*b;
  e2=e#e;
  Hat=X*invxtx*x`;
  et=e2/(1-vecdiag(hat))##2;
  cov=invxtx*x`*diag(et)*x*invxtx;
  se=sqrt(vecdiag(cov));
  print b se;
  return(cov);
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 educ exper exper2 female }`;
Print "Model A";
run ols(x1,lwage,vnameA); 
Print "OLS with HCCME0 std errors ";
cov1=hccme0(lwage,x1);
Print "OLS with HCCME3 std errors ";
cov2=hccme3(lwage,x1);
bigr = j(4,1,0)||I(4);
lr=j(4,1,0);
Print "Overall Regression Significance (Chisquare form) test based on Hccme3";
run wald2(lwage,x1,bigr,lr,cov2);
run;