*	See Gibbs estimation of microstructure models: Teaching notes
	March, 2010
	Joel Hasbrouck
____________________________________________________________________________________________________

	RollGibbsBeta01.sas

____________________________________________________________________________________________________;
options nodate nocenter nonumber ps=70 ls=120  sasautos=('./sas macros','.') mprint; 
libname this '.';

options orientation=landscape;
ods pdf file='RollGibbsBeta01.pdf';

%let Infinity=1e30;
%let eps=1e-30;

*___________________________________________________________________________________________________

	Generate some price data
____________________________________________________________________________________________________;
%let nObs=100;
%let sdu = .01;
%let c = 0.01;
%let beta = 1.1;
%let sdm = 0.01;
data sim;
	m = 0;
	pm = 0;	*	pm is the market index;
	do t=1 to &nObs;
		m = m + &sdu*normal(345699);
		pm = pm + &sdm*normal(12312);
		q = 1-2*(uniform(246)>.5);
		p = m + &beta*pm + &c*q;
		keep q p pm;
		output;
	end;
	run;

%let nSweeps=10000;


*___________________________________________________________________________________________________

	Simulate parameters and qs
____________________________________________________________________________________________________;
%RollGibbsBeta(sim, nSweeps=&nSweeps, cStart=&c, varuStart=%sysevalf(&sdu*&sdu),
	betaStart=&beta, qDraw=1);

data olsIn;
	set sim;
	r = dif(p);
	rm = dif(pm);
	run;
proc means data=olsin;
proc reg data=olsIn;
	model r = rm / noint;
	quit;
	run;

*___________________________________________________________________________________________________

	Analyze the simulated estimates
____________________________________________________________________________________________________;
%let pDrop=0.2;							*	Proportion of sweeps to drop (burn in);
%let nDrop = %sysevalf(&pDrop*&nSweeps);	*	Number of sweeps to drop;
proc means data=parmOut (where=(sweep>&nDrop)) fw=6;
	run;
goptions device=sasemf htext=1.8 ftext=complex;
proc univariate data=parmOut (where=(sweep>&nDrop));
	var c sdu beta;
	histogram c sdu beta / kernel;
run;
proc kde data=parmOut (where=(sweep>&nDrop));
	bivar  sdu c / levels out=myOut;
	run;
proc gcontour data=myOut;
	plot value1*value2=density;
	label value1='sdu' value2='c';
	title "Joint posterior with sdu=&sdu and c=&c";
	run;
	quit;