x 'cd c:\Active\CointProgs02';
title 'CointTest01: Bivariate price discovery analysis.';
options spool nocenter sasautos=(sasautos,"./SasMacros") mautosource;

*******************************************************************************
	Simulate some data;
data prices;
	m=0;
	do t=1 to 5000;	*	number of observations to simulate;
		m = m + normal(123456);
		p1 = m + normal(6789);
		output;
	end;
data prices;
	set prices;
	p2 = 4+lag2(m) + .1*normal(2468);
*	Uncomment the following lines to produce a price plot;
/*symbol1 v=none i=smooth;*/
/*goptions htext=1 ftext=simplex;*/
/*proc gplot data=prices;*/
/*	plot p1*t p2*t / overlay;*/
/*	title "Price plot";*/
/*	format t time.;*/
/*	run;*/
/*	quit;*/

*******************************************************************************
	Estimate the VECM;
*	The following is a SAS ODS (Output Delivery System) statement that
	will cause the LagLength to be saved for later use;
ods output LagLength=dsLag;
proc model data=prices outmodel=model list;
	dp1 = dif(p1);
	dp1L = lag(dp1);
	dp2 = dif(p2);
	dp2L = lag(dp2);	
	%PDL(PDL11,9,2);
	%PDL(PDL12,9,2);
	%PDL(PDL21,9,2);
	%PDL(PDL22,9,2);
	z2 = lag(p1) - lag(p2) - zmean2;
	p1 = lag(p1) + %PDL(PDL11,dp1L) + %PDL(PDL12,dp2L) + g12*z2;
	p2 = lag(p2) + %PDL(PDL21,dp1L) + %PDL(PDL22,dp2L) + g22*z2;
	fit p1-p2 / sur
		outs=cov (drop=_nused_) 
		outest=parm (keep=zmean2);
	title2 "Estimating VECM";
	quit;
	run;
*******************************************************************************
	Construct the VMA representation of the VECM by forecasting the system
	subsequent to a unit shock in each of the price variables;
%let nq=120; 				*	Number of periods ahead to forecast;
*	Build a dataset that has unit shocks in each of the price variables;
data shock;
	array p(*) p1-p2;
	set parm; 	*	Read in estimates of the zmeans.;
	set dsLag;	*	Read in the maximum lag used in the estimation;
	MaxLag = nValue1;
	do k=1 to 2;	*	loop over prices;
		p1 = 0;
		p2 = -zmean2;
		do i=1 to MaxLag; output; end;

		*	Initial shock;
		p2 = 0;
		p(k) = 1;
		p2 = p2 - zmean2;
		output;

		*	values to forecast;
		p1 = .; p2=.;
		do i=1 to &nq; output; end;
		keep k p1 p2;
	end;
proc model model=model;
	solve p1-p2 / forecast data=shock out=impResp (rename=(_lag_=t));
	title2 "Forecasting of VECM";
	by k;
	run;
*	Normalize forecast of p2 by adding back the mean difference from p1;
data ImpResp;
	if _n_=1 then set Parm; *	This dataset contains the estimates of zmean2;
	set ImpResp;
	p2 = p2 + zmean2;
	keep t p1 p2 k;
	run;

*******************************************************************************
	Graph the impulse response functions;
proc transpose data=ImpResp out=ImpRespT (rename=(col1=p));
	var p1-p2;
	by k t;
	run;
proc format;
	value k 1="Initial unit shock to p1" 2="Initial unit shock to p2";
symbol1 v=none i=smooth;
goptions htext=1 ftext=simplex;
legend1 label=none;
proc gplot data=ImpRespT;
	plot p*t=_name_ / legend=legend1;
	by k;
	format k k.;
	label k='00'x;
	title2 "Impulse response functions";
	format t time.;
	run;
	quit;

*******************************************************************************
	Perform the random walk analysis.  The actual calculations are done in the
	macro %RandomWalkAnalysis.  The following code sets up the call.
	See the %RandomWalkAnalysis macro for more details.
;
data ImpResp2;	*	Pull off the last forecasts;
	set ImpResp;
	by k;
	if last.k then output;
proc transpose data=ImpResp2 out=Coeff prefix=p;
	var p1 p2;
	run;
title2 'Random walk analysis and information shares';
%RowDifference(Coeff);
%RandomWalkAnalysis(RwOut,Cov,Coeff);
