x 'cd c:\Active\CointProgs02';	*	Set current directory;
libname this ".\Sas Datasets";
libname TempDir "c:/TempComp";	*	This will be used as a work directory;
options spool nocenter sasautos=(sasautos,"./SasMacros") mautosource nonotes nomprint;
*options notes mprint;

%let TempLib = TempDir;
%let dsPrices = Prices;		*	Price dataset (in current directory);
%let dsPriceMap = PriceMap;	*	Index dataset for the price data (in current directory);
%let nPrice=3;				*	Number of prices in the model;
%let PriceVars = spy fes fsp;
%let PriceVarsQuoted = 'SPY' 'FES' 'FSP';
%let PriceDefs = spy=p1 fsp=p2 fes=p3;	*	The price discovery routines require
											that the prices be called p1, p2, etc.;
%let nq=600; 				*	Number of periods ahead to compute impulse response
								functions.  This must be set sufficiently high that
								(for a given innovation) the forecasted price have
								converged.;


title "Price discovery analysis: &PriceDefs";

*******************************************************************************
	Build a data file (&TempLib..PriceDay) that contains all the times and
	prices for a given day.;
%macro DayBuild(iDate,PrintLevel=0);
%global pDate;
%global FirstPriceRec;
%global LastPriceRec;
%global CountPriceRec;
proc datasets lib=&TempLib nolist;
	delete temp1 temp2 PriceDay;
	quit;
data _null_;
	%put %str(iDate:) &iDate;
	set this.&dsPriceMap (firstobs=&iDate);
	call symput("FirstPriceRec",compress(put(First,12.)));
	call symput("LastPriceRec",compress(put(Last,12.)));
	call symput("CountPriceRec",compress(put(Count,12.)));
	call symput("pDate",compress(put(Date,date.)));
	stop;
	run;
%put DayBuild Date=&pDate First=&FirstPriceRec Last=&LastPriceRec Count=&CountPriceRec;
data &TempLib..temp1;
	set this.&dsPrices (firstobs=&FirstPriceRec obs=&LastPriceRec);
	if time>"9:30:00"t and time<"16:00:00"t 
		and symbol in(&PriceVarsQuoted) then output;
proc sort data=&TempLib..temp1;
	by time;
*	Align the prices by time;
proc transpose data=&TempLib..temp1 out=&TempLib..temp2 (drop=_name_);
	var price;
	id symbol;
	by time;
*	Fill in missing prices with the most recent values.
	Rename the prices to p1, p2, etc.;
proc expand data=&TempLib..temp2 out=&TempLib..PriceDay 
	(rename=(&PriceDefs))
	to=second method=step;
	id time;
data &TempLib..PriceDay;
	set &TempLib..PriceDay;
	time = timepart(time);
	format time time.;
	array p(*) p1-p&nPrice;
	do i=1 to &nPrice;
		p(i) = 10000*log(p(i));
	end;
	run;
%if &PrintLevel>0 %then %do;
	%let TimeCondition= time>"10:00:00"t and time<"10:01:00"t;
	proc print data=&TempLib..temp2;	*	Printout used in documentation;
		where &TimeCondition;
		title2 "DayBuild iDate=&iDate input data.";
		run;
	proc print data=&TempLib..PriceDay;	*	Printout used in documentation;
		where &TimeCondition;
		title2 "DayBuild iDate=&iDate expanded data.";
		run;
%end;
%mend DayBuild;

/**	Test of DayBuild macro;*/
/*%DayBuild(1,PrintLevel=1);*/

*	Generate a price graph for a single day;
%macro GraphPrices(dsIn);
proc transpose data=&dsIn out=&TempLib..PriceDayT prefix=p;
	var p1-p&nPrice;
	by time;
symbol1 v=none i=smooth;
goptions htext=1 ftext=simplex;
proc gplot data=&TempLib..PriceDayT;
	plot p1*time=_name_;
	title2 "Daily price plot";
	run;
%mend GraphPrices;
/*%GraphPrices(&TempLib..PriceDay);*/

*******************************************************************************
	Estimate the VECM;
%let ModelVECMFirstCall=1;	*	Flag to indicate that this is the first call;
%macro ModelVECM(PrintLevel=1);

%*	The SAS ODS (output delivery system) is used to save the lag length in a
	dataset.  For ODS to work, the print option must be used on the first call
	of ModelVECM.  To force this, ModelVECMFirstCall is initially set to one
	(above) and then reset.;
	%if &ModelVECMFirstCall %then ods output LagLength=&TempLib..dsLag %str(;);

%*	If you want to check the looping logic of the program and don''t want to 
	wait forever, stick an ""(obs=2000)"" parameter in the proc model statement.;
proc model data=&TempLib..PriceDay outmodel=&TempLib..Model

	%if &ModelVECMFirstCall %then %let ModelVECMFirstCall=0;

%else %if &PrintLevel<1 %then noprint;
	;

	*	Set up dp's, lagged dp's and moving averages;
	%do i=1 %to &nPrice;
		dp&i = dif(p&i);
		dp&i.L = lag(dp&i);
		%do j=1 %to &nPrice;
			%PDL(PDL&i._&j,9,2);
		%end;
		%let nMA = 2;	*	Number of moving average terms in specification;
		dp&i._L11 = lag11(dp&i);
		ma&i._1 = movavg10(dp&i._L11);
		ma&i._2 = movavg20(dp&i._L11);
	%end;
	
	*	Build deviation variables;
	%do i=2 %to &nPrice;
		z&i = lag(p1) - lag(p&i) - zmean&i;
	%end;

	*	Build equations for each price;
	%do i=1 %to &nPrice;
		p&i = lag(p&i)
			%do j=1 %to &nPrice;
				%*	Error correction term;
				%if &j>=2 %then %do;
					+ g&i._&j * z&j
				%end;
				+ %PDL(PDL&i._&j,dp&j.L)
				%*	Moving average terms;
				%do k=1 %to &nMA;
					+ A&i._&j._&k * ma&j._&k
				%end;
			%end;
		;
	%end;

	fit p1-p&nPrice /sur
		outs=&TempLib..Cov (drop=_nused_) 
		outest=&TempLib..Parm (keep=zmean2-zmean&nPrice);
	title2 "Estimating VECM";
	run;
	quit;
	ods output clear;
%mend ModelVECM;

/*%ModelVECM;*/

*******************************************************************************
	Construct the VMA representation of the VECM by forecasting the system
	subsequent to a unit shock in each of the price variables;
%macro ImpulseResponse(PrintLevel=1);
*	Build a dataset that has unit shocks in each of the price variables;
data &TempLib..Shock;
	set &TempLib..Parm; 	*	Read in estimates of the zmeans.;
	set &TempLib..dsLag;	*	Read in the maximum lag used in the estimation;
	MaxLag = nValue1;
	array p(1:&nPrice) 		p1-p&nPrice; 
	array zmean(2:&nPrice) 	zmean2-zmean&nPrice; 

	do kInitialPriceShock=1 to &nPrice; 
		p(1) = 0;
		do i=2 to &nPrice;
			p(i) = -zmean(i);
		end;
		do i=1 to MaxLag; output; end;
		if kInitialPriceShock=1 then p(kInitialPriceShock) = 1;
		else p(kInitialPriceShock) = 1-zmean(kInitialPriceShock);
		output;
		do i=1 to &nPrice; p(i) = .; end;
		do i=1 to &nq; output; end;
	end;
proc model model=&TempLib..Model
	%if &PrintLevel<1 %then noprint;
	;
	solve p1-p&nPrice / forecast data=&TempLib..Shock out=&TempLib..ImpResp (rename=(_lag_=t));
	title2 "Forecasting VECM to get impulse response functions";
	by kInitialPriceShock;
	run;
*	Normalize forecasts of p2...p&nPrice by adding back the mean difference from p1;
data &TempLib..ImpResp;
	array p(1:&nPrice) 		p1-p&nPrice; 
	array zmean(2:&nPrice) 	zmean2-zmean&nPrice; 
	set &TempLib..ImpResp;
	if _n_=1 then set &TempLib..Parm; *	This dataset contains the estimates of zmean2;
	do i=2 to &nPrice;
		p(i) = p(i) + zmean(i);
	end;
	keep t p1-p&nPrice kInitialPriceShock;
	run;
%if &PrintLevel>0 %then %do;
	proc print data=&TempLib..ImpResp;
		by kInitialPriceShock;
		title3 "Impulse response output";
		run;
%end;
%mend ImpulseResponse;

/*%ImpulseResponse;*/

*******************************************************************************
	Graph the impulse response functions;
%macro GraphImpulseResponse(dsIn);
proc transpose data=&dsIn out=&TempLib..ImpRespT (rename=(col1=p));
	var p1-p&nPrice;
	by kInitialPriceShock t;
	run;
proc format;
	value kShock 
		%do i=1 %to &nPrice;
			&i="Initial unit shock to p&i"
		%end;
		;
symbol1 v=none i=smooth;
goptions htext=1 ftext=simplex;
legend1 label=none;
proc gplot data=&TempLib..ImpRespT;
	plot p*t=_name_ / legend=legend1;
	by kInitialPriceShock;
	format kInitialPriceShock kShock.;
	label kInitialPriceShock='00'x;
	title2 "Impulse response functions";
	format t time.;
	run;
	quit;
%mend GraphImpulseResponse;

/*%GraphImpulseResponse;*/

*******************************************************************************
	Perform the random walk / information share analysis;
%macro RandomWalkAnalysis2(PrintLevel=1);
data ImpResp2;	*	Pull off the last forecasts;
	set &TempLib..ImpResp;
	by kInitialPriceShock;
	if last.kInitialPriceShock then output;
proc transpose data=ImpResp2 out=Coeff prefix=p;
	var p1-p&nPrice;
*	Make sure that the rows are approximately identical;
%RowDifference(Coeff,PrintLevel=0);
data _null_;
	%if %eval(&MaxRelDiff<0.001) %then %do;
		%put Warning.  RandomWalkAnalysis2;
		%put Maximum relative row difference in coefficient matrix is &MaxRelDiff;
	%end;
	run;
%RandomWalkAnalysis(&TempLib..RwOut,&TempLib..Cov, Coeff, PrintLevel=&PrintLevel, AllRows=0);
%if &PrintLevel>0 %then %do;
proc print data=&TempLib..RwOut;
	title2 "Output from RandomWalkAnalysis";
	run;
%end;
%mend RandomWalkAnalysis2;

/*%RandomWalkAnalysis2;*/

*******************************************************************************
	Perform the analysis (loop) over all days in the sample;
%macro DoAllDates;
	proc datasets lib=this nolist;
		delete RWSummary IRSummary;
		quit;
	proc datasets lib=&TempLib nolist;
		delete PriceDay Model Cov Parm Shock ImpResp dsLag;
		quit;
	%NumObs(this.PriceMap,nDates);
	%put %str(nDates:) &nDates in &TempLib..PriceMap;
	%let PrintLevel=1;		*	This will print out details for first iteration (iDate=1);
	%*let nDates=2;			*	Set to 1 or 2 for debugging.;
	%do iDate=1 %to &nDates;*	Loop over all dates;
		%CurrentDateTime(cd);
		%put Starting iDate loop for iDate= &iDate of &nDates. Starting at &cd;

		*	The following macros build the dataset and do the estimations for the day.;
		%DayBuild(&iDate);
		%ModelVECM(PrintLevel=&PrintLevel);
		%ImpulseResponse(PrintLevel=&PrintLevel);
		%RandomWalkAnalysis2(PrintLevel=&PrintLevel);
		%ComputeCorrelations(&TempLib..Cov, &TempLib..Corr);

		*	Build a summary dataset of estimates for day;
		data &TempLib..DaySummary;
			set	&TempLib..RwOut
				&TempLib..Corr (rename=(_type_=type));
			rename _name_=name;
			run;
		*	Add iDate to the summary dataset;
		data &TempLib..DaySummary;
			set &TempLib..DaySummary;
			iDate = &iDate;
		*	Add iDate to the impulse response dataset;
		data &TempLib..ImpResp;
			set &TempLib..ImpResp;
			iDate = &iDate;
		proc append base=this.RWSummary data=&TempLib..DaySummary;
		proc append base=this.IRSummary data=&TempLib..ImpResp;
		%let PrintLevel=0;	%*	Don''t print out any more detail;
		run;
	%end;
	%CurrentDateTime(cdEnd);
	%put iDate loop completed at &cdEnd;
%mend DoAllDates;

*******************************************************************************
	Run DoAllDates;
%DoAllDates;

*******************************************************************************
	Summarize the numeric results;

%CurrentDateTime(d);
title2 "Summarizing &d";
proc sort data=this.RWSummary;
	by type name idate;
proc means data=this.RWSummary (drop=iDate) noprint;
	var p1-p&nPrice VarRandWalk;
	output out=DescStat;
	by type name;
proc sort data=DescStat;
	by type _stat_ name;
proc print data=DescStat (drop=_type_ _freq_) noobs;
	by type;
	run;

*******************************************************************************
	Average and graph impulse response functions;
proc sort data=this.IRSummary;
	by kInitialPriceShock t iDate;
proc means data=this.IRSummary noprint;
	var p1-p&nPrice;
	output out=IRMeans mean=;
	by kInitialPriceShock t;
	run;
%GraphImpulseResponse(IRMeans);