PROC IMPORT OUT= WORK.one 
            DATAFILE= "C:\DATA\SAS\F5.1.csv" 
            DBMS=CSV REPLACE;
	 /* Year qtr realgdp realcons realinvs realgovt realdpi cpi_u M1 tbilrate unemp pop infl realint */
     GETNAMES=YES;
     DATAROW=2; 
run;
data two;
set one;
lagc=lag(realcons);
int=1;
proc model data = two;
parms a b c;
var realcons realgdp lagc;
realcons =a+b*realdpi+c*lagc;
fit realcons;
test (b/(1-c))=1;
run;

proc iml;
use two;
read all var{realcons} into y where(lagc ^= .);
read all var{int realdpi lagc} into X where(lagc ^= .);

/*  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 dpi c_1}`;
Print "Model A";
run ols(x,y,vnameA); 
run;