options nocenter ls=80;

libname temp 'd:\empiricalmethods\programs\';
filename tempdir 'd:\empiricalmethods\programs\';

data daily;
  informat date yymmdd8.;
  format date mmddyy8.;

  infile tempdir(dailyindexreturns.txt) firstobs=6 missover;
  input permno date ret vwret ewret;
  drop permno;
* calculate log returns;
  vw=log(1+vwret);
  ew=log(1+ewret);
  lvw=lag(vw);
  lew=lag(ew);
if date<mdy(7,31,1962) or year(date)>1994 then delete;
proc print data=daily (obs=100);
proc univariate;
  var vw ew;
run;

proc autoreg;
  model ew=/all nlag=4;
run;

proc iml;
  use daily;
  read all var {ew} into ew;
  q=4;
  t=nrow(ew);
  print t;
  one=j(t,1,1);
  r1=j(q-1,1,0);
  d1=j(q-1,1,0);

*compute mean;
  mu=one`*ew/t;
print mu;

* compute variance;
  dev1=ew-mu;
  dev0=dev1#dev1;
  var=one`*dev0/(t);
  print var;

*compute k order autocorrelation, variance, and het consistent standard error;
  do k=1 to q-1;
  do z=k+1 to t by 1;
    r1[k]=r1[k]+(ew[z]-mu)*(ew[z-k]-mu);
	d1[k]=d1[k]+(ew[z]-mu)**2*(ew[z-k]-mu)**2;
  end;
  end;
  r1=r1/t;
  print r1;
  rho=r1/var;
  d1=t*d1/(t*var)**2;

  vrq=0;
  thetaq=0;
  do k=1 to q-1;
    vrq=vrq+(1-k/q)*rho[k];
	thetaq=thetaq+(1-k/q)**2*d1[k];
  end;
  vrq=1+2*vrq;
  thetaq=4*thetaq;
  psi=sqrt(t)*(vrq-1)/sqrt(thetaq);
  print vrq rho thetaq psi;

quit;