options nocenter ls=90;

libname temp 'c:\empirical_methods_class\programs\';
filename tempdir 'c:\empirical_methods_class\programs\';

data all (keep=date vwretd ewretd r1-r10 year month);
  informat date yymmdd8.;
  format date date7.;

  infile tempdir(decilereturns_nyseamex.txt) missover firstobs=2 dlm=',';
  input date vwretd vwretx ewretd ewretx r1-r10;
  if vwretd=-99 then vwretd=.;
  if r1=-99 then r1=.;
  if r2=-99 then r2=.;
  if r3=-99 then r3=.;
  if r4=-99 then r4=.;
  if r5=-99 then r5=.;
  if r6=-99 then r6=.;
  if r7=-99 then r7=.;
  if r8=-99 then r8=.;
  if r9=-99 then r9=.;
  if r10=-99 then r10=.; 
  if vwretd=. then delete;
  if year(date)<1926 or year(date)>1999 then delete;
  year=year(date);
  month=month(date);
proc means;
proc sort;
  by year month;
run;
 
data tbill (keep=date t30ret t90ret year month);
  informat date yymmdd8.;
  format date date7.;

  infile tempdir(riskfree.txt) missover firstobs=2 dlm=',';
  input date b30ret b30ind b20ret b20ind b10ret b10ind b7ret b7ind b5ret b5ind
  b2ret b2ind b1ret b1ind t90ret t90ind t30ret t30ind cpiret cpiind;
  if t30ret=-99 then t30ret=.;
  if t90ret=-99 then t90ret=.;
  if year(date)<1926 or year(date)>1999 then delete;
  year=year(date);
  month=month(date);
proc means;
proc sort;
  by year month;
run;

data all;
  merge all (in=A) tbill;
  by year month;
  if A;
  r1=r1-t30ret;
  r2=r2-t30ret;
  r3=r3-t30ret;
  r4=r4-t30ret;
  r5=r5-t30ret;
  r6=r6-t30ret;
  r7=r7-t30ret;
  r8=r8-t30ret;
  r9=r9-t30ret;
  r10=r10-t30ret;

  rvw=vwretd-t30ret;
  if year<1965 then delete;
  if year>1994 then delete;
proc sort;
  by date;
*proc print;
proc means;
run;

data market (keep=year month date rvw);
  set all;
run;

proc means data=market noprint;
  var rvw;
  output out=mkt mean=mu std=sigmam n=t;
run;

data mkt;
  set mkt;
  sigmam=sqrt((t-1)/t)*sigmam;
proc print data=mkt;
run;
  
proc model data=all;
  var r1 r10 rvw;
  parms a1-a10 b1-b10;
  r1=a1+b1*rvw;
  r2=a2+b2*rvw;
  r3=a3+b3*rvw;
  r4=a4+b4*rvw;
  r5=a5+b5*rvw;
  r6=a6+b6*rvw;
  r7=a7+b7*rvw;
  r8=a8+b8*rvw;
  r9=a9+b9*rvw;
  r10=a10+b10*rvw;
  fit r1 r2 r3 r4 r5 r6 r7 r8 r9 r10/sur outs=sigma outest=est1;
  test a1, a2, a3, a4, a5, a6, a7, a8, a9, a10,/wald;
run;

proc print data=sigma;
run;
proc print data=est1;
run;
proc print data=cov;
run;

proc iml;
  use sigma;
  read all var {r1 r2 r3 r4 r5 r6 r7 r8 r9 r10} into sigmae;
  print sigmae;

  use est1;
  read all var {a1 a2 a3 a4 a5 a6 a7 a8 a9 a10} into alpha;
  print alpha;

  use mkt;
  read all var {mu} into mu;
  read all var {sigmam} into sigma;
  read all var {t} into t;
  print mu;
  print sigma;
  print t;

  isigmae=inv(sigmae);
  j0=t*alpha*isigmae*alpha`/(1+(mu*mu)/(sigma*sigma));
  j1=(t-10-1)/10*((1+(mu**2)/(sigma**2))**-1)*alpha*isigmae*alpha`;
  print j0 j1;
quit;
  
proc model data=all;
  parms a1-a10 b1-b10;
  r1=a1+b1*rvw;
  r2=a2+b2*rvw;
  r3=a3+b3*rvw;
  r4=a4+b4*rvw;
  r5=a5+b5*rvw;
  r6=a6+b6*rvw;
  r7=a7+b7*rvw;
  r8=a8+b8*rvw;
  r9=a9+b9*rvw;
  r10=a10+b10*rvw;
  fit r1 r2 r3 r4 r5 r6 r7 r8 r9 r10/gmm outest=est1 outcov covb;
  test a1, a2, a3, a4, a5, a6, a7, a8, a9, a10,/wald;
run;

data est;
  set est1;
  if _n_=1;
proc print data=est;
run;

data vcov;
  set est1;
  if _n_=1 then delete;
  if _n_>11 then delete;
proc print data=vcov;
run;

proc iml;
  use vcov;
  read all var {a1 a2 a3 a4 a5 a6 a7 a8 a9 a10} into vcov;
  print vcov;

  ivcov=inv(vcov);
print ivcov;

  use est;
  read all var {a1 a2 a3 a4 a5 a6 a7 a8 a9 a10} into alpha;
  print alpha;

  j7=alpha*ivcov*alpha`;
  print j7;
run;