*	This SAS program performs the simulations and estimations reported in 
	"Stalking the efficient price in market microstructure specifications: An overview";

x 'cd c:\auser\docs\Cointegration JFM\Sas Programs';
options nocenter nosymbolgen mprint pagesize=54 ls=124;

%let nObs=50000;	*	number of observations to use in simulations;
%let np=20;			*	Order of estimated VECM;
%let nq=60;			*	Order of the cumulative impulse computation (VMA);

*	Description

	This code below consists of an estimation program, followed by the
	generation and estimation of the paper's models I, II and III.

	The model estimated in this program is the VECM:

	Dp(t) = A(1) Dp(t-1) + A(2) Dp(t-2) + ... + A(np) Dp(t-np) 
				+ g ([1 -1] x p(t-1)) + u(t)							(1)
	where:
	p(t)		2x1 vector of prices [p1(t) p2(t)]`
	Dp(t)		2x1 first difference
	A(1) ... A(np) are (2x2) coefficient matrices
	g			2x1 vector of adjustment coefficients
	u(t)		2x1 vector of disturbances

	The estimation procedure is implemented as a macro "PriceDiscoveryAnalysis".
	This macro assumes that a work dataset named "t1" has two price series ("p1" and "p2").
	Estimation is performed using the MODEL procedure.  Additional necessary calculations
	are performed using the IML ("Interactive Matrix Language") procedure.  The program 
	makes extensive use of SAS's macro facility.

	Programming notes:

	1.  Because the MODEL procedure does not allow a first difference on the lhs of an equation,
	the specification actually estimated is:

	p(t) = p(t-1) + A(1) Dp(t-1) + A(2) Dp(t-2) + ... + A(np) Dp(t-np) 
				+ g ([1 -1] x p(t-1)) + u(t)							(2)

	2.	The random-walk / information share calculation relies on the VMA representation:
		Dp(t) = Psi(L) u(t) 
	where Psi(L) is a matrix lag polynomial.  Only the sum is actually needed, that is,
	Psi(1).  This is obtained from the long-run cumulative impact function (from (2)) in
	response to unit shocks in p1 and p2.  The cumulative impact function is computed
	using the "forecast" component of the MODEL procedure.  The impact is computed through
	lead time "nq".

	3.	The estimation procedure in MODEL generates:
				The disturbance covariance matrix, Cov(u)
				Psi(1), described above
				The estimated adjustment coefficients, g.
	These are input to an IML ("Interactive Matrix Language") procedure that actually
	computes the random walk decompositions, information shares, and
	Granger-Gonzalo factor weights.

	4.	The MODEL procedure does not allow arrays of parameters.  The %Lz2 and %A macros
	facilitate a work-around.  An array-like reference "%A(6,1,2)" is expanded as
	"A060102".  That is, A(i,j,k) is the coefficient in equation j, of variable k, at lag i.

	5.	The %vset macro generates a sum where each term has the form:
		A(i,j,k) * Dp[k](t-i).
	This enables the matrix model (equation (2), above) to be specified without writing
	out each term explicitly.  The "list" option on the proc model statement can be used
	to print out the full expanded model that results.

	6.	SAS's IML package does not have a routine to compute a basis for the null space
	of a matrix.  (This is needed for the GG factor computation.)  I implemented an
	equivalent procedure using a singular value decompositon.  (See "Numerical recipes in C".)
;
%macro Lz2(i);
	%if	&i<10 %then 0&i;
	%else &i;
%mend Lz2;

%macro A(i,j,k);
%unquote(A%Lz2(&i)%Lz2(&j)%Lz2(&k))
%mend A;

*	The %vset macro generates sums of terms of the form: A(i,j,k) * lag(i, dif(v));
%macro vset(ieq,iv,LagStart,LagEnd,v);
	%do i=&LagStart %to &LagEnd;
		+ %A(&i,&ieq,&iv) * lag&np(&i,dif(&v))
	%end;
%mend vset;

*	The PriceDiscoveryAnalysis macro estimates a VECM and performs IS and PT/GG analyses.
*	Dataset "t1" is assumed to contain p1 and p2;
%macro PriceDiscoveryAnalysis;
proc datasets lib=work; *	Delete any old work datasets that are lying around.;
	delete shock1 shock2 t2 t3 GGWts PsiSumRow acv;
	run;
*	
	The shock1 and shock2 datasets generated in the following step
	will be used to generate the impulse response functions
	(VMA coefficients) by using the FORECAST option in proc model.;
data shock1 shock2;
*	lagged values assumed be zero;
	p1 = 0;	p2 = 0;
	do i=1 to &np; output; end;
*	impulses;
	i=&np+1;
	p1 = 1;	p2 = 0;
	output shock1;
	p1 = 0;	p2 = 1;
	output shock2;
*	future values;
	p1 = .; p2 = .;
	do i=&np+2 to &np+&nq+1; output; end;
	run;
proc model data=t1;
	p1 = lag(p1) + %vset(1,1,1,&np, p1) + %vset(1,2,1,&np,p2) + g1*(lag(p1)-lag(p2));
	p2 = lag(p2) + %vset(2,1,1,&np, p1) + %vset(2,2,1,&np,p2) + g2*(lag(p1)-lag(p2));
	fit p1 p2 / sur
		outs=vcov (drop=_nused_) 
		outest=est (keep=g1 g2) 
		out=res (rename=(p1=e1 p2=e2)); 
	solve / forecast data=shock1 out=ma1;
	solve / forecast data=shock2 out=ma2;
	title2 "Estimation and forecasting of VECM";
proc iml;
	start;
	print "Random-walk / information share analysis ordered to place p1 innovations first.";
	use ma1; read all var {p1 p2} into ma1 [colname=xname];
	use ma2; read all var {p1 p2} into ma2 [colname=xname];
	PsiSum = (ma1[nrow(ma1),] // ma2[nrow(ma2),])`;	*	Take only the last rows;
	print "Note: All rows should be the same-> "PsiSum;
	use vcov; read all into vcov [colname=cvname];
	print "VECM residual covariance matrix: " vcov;
	vw1 = PsiSum * vcov * PsiSum`;
	PsiSumRow = PsiSum[1,];
	*	Write weights to a dataset;
	create PsiSumRow var {w1 w2}; append from PsiSumRow;
	g = (PsiSumRow*t(root(vcov)))##2;	*	"root" performs the cholesky factorization.;
	vw = sum(g);
	print "Variance of rw component:" vw;
	is1 = g/vw;
	print "[Information share of p1; share of p2]: " is1;

	*	permute the variables to reverse the order;
	print "Random-walk / information share analysis ordered to place p2 innovations first.";
	vcovP = (vcov[2,2]||vcov[2,1])//(vcov[1,2]||vcov[1,1]);
	PsiSumRowP = PsiSumRow[1,2] || PsiSumRow[1,1];
	g = (PsiSumRowP*t(root(vcovP)))##2;
	vw = sum(g);
	print "Variance of rw component:" vw;
	is2 = g/vw;
	print "[Information share of p2; share of p1]: " is2;

	*	Granger-Gonzalo factors;
	use est; read all into g;	*	estimated cointegrating vector coefficients;
	print "Cointegrating vector coefficients:" g;
	call svd(u,q,v,g);	*	called to generate null space of g;
	if abs(q[2,1])>1e-10 then do;
		print "Error in null space analysis. svd q=" q;
		stop;
	end;
	GGWts = v[,2]`;
	GGWts = GGWts/sum(GGWts);	*	Normalize;
	print "Granger-Gonzalo factor weights:" GGWts;
	*	Write the weights to a dataset;
	create GGWts var {w1 w2};	append from GGWts;
	finish;
	run;
	quit;
data t2;	*	Construct GG common factor;
	retain w1 w2;
	if _n_=1 then set GGWts;
	set t1;
	f = w1*p1 + w2*p2;
proc arima data=t2;	*	Analyze GG common factor;
	identify var=f(1) nlag=5 noprint outcov=acv (keep=lag var n cov corr);
proc print data=acv noobs;
	title2 "Autocovariances of changes in GG common factor.";
	run;
data t3;	*	Construct (changes in) implicit random walk;
	retain w1 w2;
	if _n_=1 then set PsiSumRow;
	set res;	*	"res" contains the VECM residuals;
	dm = w1*e1 + w2*e2;
proc arima data=t3;
	identify var=dm nlag=5 noprint outcov=acv (keep=lag var n cov corr);
proc print data=acv noobs;
	title2 "Autocovariances of implicit random-walk.";
	run;
%mend PriceDiscoveryAnalysis;
*
	This ends the definition of the macros.  The remaining code simulates
	data for each of the models I, II and II, and invokes for each simulated
	dataset the PriceDiscoveryAnalysis macro.
;

/*title "(Hasbrouck) Price discovery example I: A two-market ""Roll"" model.";*/
/**	Generate the data;*/
/*data t1;*/
/*	sdu = 1;*/
/*	c=1;*/
/*	m = 0;*/
/*	do t=1 to &nObs;*/
/*		u = sdu*normal(123456);*/
/*		m = m+u;*/
/*		q1 = 1-2*ranbin(654321,1,0.5);*/
/*		q2 = 1-2*ranbin(678910,1,0.5);*/
/*		p1 = m+q1*c;*/
/*		p2 = m+q2*c;*/
/*		output;*/
/*		keep p1 p2;*/
/*	end;*/
/*	run;*/
/*%PriceDiscoveryAnalysis;*/

/*title "(Hasbrouck) Price discovery example II: Two markets with private informations.";*/
/**	Generate the data;*/
/*data t1;*/
/*	sdu = 1;*/
/*	c1=1;*/
/*	c2=1;*/
/*	m = 0;*/
/*	lambda = 1;*/
/*	do t=1 to &nObs;*/
/*		q1 = 1-2*ranbin(654321,1,0.5);*/
/*		q2 = 1-2*ranbin(678910,1,0.5);*/
/*		p2 = m+q2*c2;*/
/*		m = m + lambda * q1;*/
/*		u = sdu*normal(123456);*/
/*		p1 = m+q1*c1;*/
/*		output;*/
/*		keep p1 p2;*/
/*	end;*/
/*	run;*/
/*%PriceDiscoveryAnalysis;*/

title "(Hasbrouck) Price discovery example III: Two markets with public and private information.";
proc datasets lib=work;
	delete t1;
	run;
*	Generate the data;
data t1;
	sdu = 1;
	c1=1;
	c2=0;
	m = 0;
	lambda = 1;
	do t=1 to &nObs;
		q1 = 1-2*ranbin(654321,1,0.5);
		q2 = 1-2*ranbin(678910,1,0.5);
		p2 = m+q2*c2;
		u = sdu*normal(1234656);
		m = m + lambda * q1 + u;
		p1 = m+q1*c1;
		output;
		keep p1 p2;
	end;
	run;

%PriceDiscoveryAnalysis;