ods html close;

%let stocks =mmm aa axp t bac ba cat cvx csco ko dd xom ge hpq hd intc ibm jnj jpm mcd mrk msft pfe pg trv utx vz wmt dis;

libname dow30 "C:\temp\9.3\dow30";

proc datasets lib=work nolist kill;
quit;

/*Download stocks and prices*/
%get_stocks(&stocks,01JAN1998,30JUL2012,keepPrice=1);

/*output prices and returns to permanent location*/
data dow30.returns;
set returns;
run;

/*Put Week, Quarter, and Year on returns for aggregation*/
data returns;
format week qtr year 4.;
set returns;
week = week(date);
year = year(date);
qtr = qtr(date);
run;

data dow30.prices;
format week qtr year 4.;
set prices;
week = week(date);
year = year(date);
qtr = qtr(date);
run;

data dow30.prices_wk;
set dow30.prices;
by year week;
if last.week;
run;

/*Sum returns over each week to create weekly returns*/
proc summary data=returns;
var &stocks;
by year week;
output out=weekly(drop=_type_ _freq_) sum=;
run;

/*Copy weekly returns to perm library*/
proc datasets lib=work nolist;
copy out=dow30;
select weekly ;
run;
quit;

/*Add a record index*/
data weekly(index=(indx));
set weekly;
indx = _N_;
run;

/*Select the index values for each quarter end, 
  starting with year end 1999, and ending with
  year end 2011 */
proc sql noprint;
select indx, year, week
into :starts separated by ' '
from weekly
where (2012 > year > 1999
  		and week in (13,26,39,52)
	   )
  or (year = 1999 and week = 52)
order by year, week;
quit;

%put &starts;

/*Macro will loop over the starting index values 
  starts = period starting index values
  obs = number of observations for copula fitting
  fwd = number of periods forward to simulate
  draws = number of draws for simulation
  max = number of times to loop, m<1 means loop over all
        values in starts */
%macro loop(starts,obs=100,fwd=1,draws=100,max=0);
/*Total draws needed in sumulation*/
%let nDraws = %eval(&fwd*&draws);
/*Number of items in starts*/
%let nSims = %sysfunc(countw(&starts));

/*update max with nSims if max < 1*/
%if %sysevalf(&max < 1) %then 
	%let max=&nSims;
;

/*Delete optES if it exists*/
proc datasets lib=work nolist;
delete optES;
run;
quit;

/*Delete cov if it exists*/
proc datasets lib=dow30 nolist;
delete cov;
run;
quit;


/*Start main loop*/
%do i=1 %to &max;
	/*start= current starting index*/
	%let start = %scan(&starts,&i);
	/*e= ending index*/
	%let e = %eval(&start - &obs);

	/*create subset of weekly returns
	  get the as-of week and year*/
	data test;
	set weekly(where=(&e <= indx <=&start)) end=last;
	if _n_ = 1 then do;
       call symput("foy",put(year,4.));
	   call symput("fow",put(week,2.));
	end;
    if last then do;
	   call symput("aoy", put(year,4.));
	   call symput("aow", put(week,2.));
	end;
	run;

	/*de-mean the weekly returns over the subset*/
	proc standard mean=0 data=test out=test;
	var &stocks;
	run;

	/*Create the covariance matrix for the subset, 
	  store in perm location */
	proc corr data=test 
		out=cov(where=(_type_="COV")) cov noprint;
	var &stocks;
	run;

	data cov;
	format year week 4.;
	set cov;
	year=&aoy;
	week=&aow;
	run;

	proc append base=dow30.cov data=cov force;
	run;

	%put AS OF = &start;
	%put From &foy, Week &fow: To &aoy, Week &aow;

	ods select NONE ;

	/*Fit a T copula and simulate from it*/
	proc copula data=test;
	var &stocks;
	fit t /
		marginals=empirical
		method=MLE;

	simulate /
		ndraws=&nDraws
		seed=54321
		out=sim;
	run;

	/*Check the sim data.  If no observations then the 
	  fitting failed.  Fail over to a Normal/Gaussian 
	  Copula */
	proc sql noprint;
	select count(*) format=5. into :nsim from sim;
	quit;

	%if &nsim = 0 %then %do;
		%put NSIM=&nsim, switching to NORMAL Model;
		proc copula data=test ;
		var &stocks;
		fit normal / 
			marginals=empirical;

		simulate /
			ndraws=&nDraws
			seed=54321
			out=sim;
		run;
	%end;
	ods select default;

	/*Create an iteration number in the simulations*/
	data sim_&aoy._%left(&aow);
	format iter 8.;
	set sim;
	iter = mod(_N_,&draws);
	run;

	proc sort data=sim_&aoy._%left(&aow);
	by iter;
	run;

	/*Aggregate the iterations into a final number for 
	  the forward period*/
	proc summary data=sim_&aoy._%left(&aow);
	var &stocks;
	by iter;
	output out=sim_&aoy._%left(&aow)(drop=_type_ _freq_ iter) sum=;
	run;

	/*Use IML for the NL Optimization */
	proc iml noprint;

	/*Read the simulations into a matrix*/
	use sim_&aoy._%left(&aow);
	read all var _num_ into sim[colname=cnames];
	close sim_&aoy._%left(&aow);

	s = ncol(sim);
	w0 = repeat(1/s,s);

	/*Function calculating expected shortfall*/
	start F_ES(x) global(sim);
	   	v = sim*x`;
	   	call sort(v);
		cutoff = floor(nrow(sim)*.05);
		es = v[1:cutoff][:];
	   return (-es);
	finish F_ES;

	/*Setup constraints*/
	lb = repeat(0,1,s);
	ub = repeat(.1,1,s);
	ones = repeat(1,1,s);
	addc = ones || {0 1};

	con = (lb || {. .}) // 
	      (ub || {. .}) // 
	      addc;
	optn = {0 0}; 
	x = w0`;

	/*Call the nl quasi-Newton optimizer*/
	call nlpqn(rc,w,"F_ES",x,optn,con);

	/*Put the returns into a data set*/
	create results_&aoy._%left(&aow) from w[colname=cnames];
	append from w;
	close results_&aoy._%left(&aow);

	quit;

	/*Add the week and year to the results*/
	data results_&aoy._%left(&aow);
	year = &aoy;
	week = &aow;
	set results_&aoy._%left(&aow);
	run;

	/*Append the results into the final results*/
	proc append data=results_&aoy._%left(&aow) base=optES force;
	run;

%end;


/*Delete the simulation and result files*/
proc datasets lib=work nolist;
delete sim: results:;
run;
quit;

%mend;


/*Call the simulations and optimizations
  Use 104 weeks (2 years) in each subset
  simulate 13 weeks (1 qtr) forward
  draw 5000 iterations   */
%let st = %sysfunc(datetime());
options nonotes nomprint;
*loop(starts,obs=100,fwd=1,draws=100,max=0);
%loop(&starts,obs=104,fwd=13,draws=5000,max=0);
options notes nomprint;
%let el = %sysevalf(%sysfunc(datetime()) - &st);
%put Took: &el;

/*Copy the final results and covariance matrices to 
  output locationthe */
proc datasets lib=work nolist;
copy out=dow30;
select optES ;
run;
quit;

/*Download DIA and SPY returns*/
%get_stocks(dia spy,03JAN2000,30JUL2012,keepPrice=0);

/*Finally, export the optimal portfolio weights
  into R and save the DataFrame*/
proc iml;
call ExportDataSetToR("dow30.optES", "optES" );
call ExportDataSetToR("dow30.cov","cov");
call ExportDataSetToR("dow30.prices","prices");
call ExportDataSetToR("dow30.prices_wk","prices_wk");
call ExportDataSetToR("returns","dia_spy");
submit / R;
	save(optES, file="C:\\temp\\9.3\\dow30\\optES");
	save(cov, file="C:\\temp\\9.3\\dow30\\cov");
	save(prices, file="C:\\temp\\9.3\\dow30\\prices");
	save(prices_wk, file="C:\\temp\\9.3\\dow30\\prices_wk");
	save(dia_spy, file="C:\\temp\\9.3\\dow30\\dia");
endsubmit;
quit;