options nocenter ls=80;

*put your directories here;
libname temp 'd:\empirical methods\';
filename tempdir 'd:\empirical methods\';

*read ibbotson data;
data ibbotson;
  informat date mmddyy10.;
  format date mmddyy10.;
  infile tempdir(stock_bond_data.csv) firstobs=2 dlm=',' missover;
  input date sp500 small tbill t1yr eqprem tint clong tlong inf;
  proc print data=ibbotson (obs=100);
run;

*compute annualized mean and standard deviation from monthly data;
proc means data=ibbotson noprint;
  var sp500 small tbill t1yr eqprem tint clong tlong inf;
  output out=sumstats mean=sp500 small tbill t1yr eqprem tint clong tlong inf
  std=vsp500 vsmall vtbill vt1yr veqprem vtint vclong vtlong vinf;
run;

data sumstats;
set sumstats;
drop _type_ _freq_;
 fa=12;
 fas=sqrt(fa);
*annualize means;
 sp500=fa*sp500;
 small=fa*small;
 tbill=fa*tbill;
 t1yr=fa*t1yr;
 eqprem=fa*eqprem;
 tint=fa*tint;
 clong=fa*clong;
 tlong=fa*tlong;
 inf=fa*inf;
* annualize std;
 vsp500=fas*vsp500;
 vsmall=fas*vsmall;
 vtbill=fas*vtbill;
 vt1yr=fas*vt1yr;
 veqprem=fas*veqprem;
 vtint=fas*vtint;
 vclong=fas*vclong;
 vtlong=fas*vtlong;
 vinf=fas*vinf;
 proc print;
 run;

 *compute compounded returns over various horizons;
data compound (keep=date csp500 ctbill nrets);
  set ibbotson;

  *input desired horizon in months (e.g., 12 months for 1 year);
  horizon=12;
  *count observations;
  i=_n_;
  nrets=horizon-mod(i,horizon);
 
  if mod(i,horizon)=1 then do;
    csp500=1+sp500/100;
	ctbill=1+tbill/100;
  end;
  else do;
    csp500=csp500*(1+sp500/100);
	ctbill=ctbill*(1+tbill/100);
  end;
  if mod(i,horizon)=0 then output compound;
  retain csp500 ctbill;
proc means;
*proc print;
run;

data compound;
  set compound;
*get rate of return for the horizon;
  csp500=csp500-1;
  ctbill=ctbill-1;
  cprem=csp500-ctbill;

*produce confidence intervals for variables;
proc means data=compound n mean std p5 p95;
      var csp500 ctbill cprem;
run;