*****************************************************************************************************************************************************************; 
*	Last Update: March 2014																	*;
* 	This SAS code creates a SAS dataset called 'Results' with OLS parameter estimates and standard errors that are robust to clustering along two dimensions*;
* 	as described in Cameron, Gelbach and Miller: Robust inference with multiway clustering, Journal of Business and Economic Statistics, 2011.		*;
*																				*;
* 	This code was written by John McInnis at U Texas - Austin (https://webspace.utexas.edu/johnmac/www/).							*;																	
*	I made some very minor modifications, which I feel will facilitate the use of his code. 								*;
*																				*;
* 	You only need to change:																*;																
* 	- name of your dataset ('yourdata'), 															*;														
*	- your cluster identifiers ('firmid' and 'time'), 													*;													
* 	- dependent and independent variables ('y' and 'x'),													*;
*	in the "CUSTOMIZE CODE HERE"-section.															*;
*****************************************************************************************************************************************************************; 


*CUSTOMIZE CODE HERE;
%let yourdata=z;

%let firmid=permno;
%let time=yyyymm;

%let y=sue;
%let x=rec ln_size ln_bm cum_return AAC_int y2-y26;


*cluster by first dimension (e.g., firm);
proc surveyreg data=&yourdata;
	cluster &firmid;
	model &y = &x /covb;
	ods output covb=firm;
run;
quit;

*cluster by second dimension (e.g., year);
proc surveyreg data=&yourdata;
	cluster &time;
	model &y = &x /covb;
	ods output covb=year;
run;
quit;

*cluster by intersection of the two dimensions (e.g, firm-year);
proc surveyreg data=&yourdata;
	cluster &firmid &time;
	model &y = &x /covb;
	ods output covb=both; 
	ods output parameterestimates=parm;
run;
quit;

*keeps original parameter estimates;
data parm; set parm;
	keep parameter estimate;
run;

*returns a dataset with a scalar for the dimensions of the var/cov matrix. This is needed to extract the square roots of the diagonals later on;
data parm1; set parm;
	n=_n_;
	m=1;
	keep m n;
run;

data parm1; set parm1; by m;
	if last.m;
	keep n;
run;

*uses matrix algebra interface to construct Var-cov matrix and extract the standard errors;
proc iml;
	use both; 
	read all var _num_ into Z;
	print Z;
	use firm; 
	read all var _num_ into X;
	print X;
	use year; 
	read all var _num_ into Y;
	print Y;
	use parm1; 
	read all var _num_ into n;
	print n;

	B=X+Y-Z; 
	C=I(n);
	D=J(n,1);
	E=C#B;
	F=E*D;
	G=F##.5;

	print B;
	print G;
	create b from G [colname='stderr']; ;
	append from G;
quit;

*creates a dataset called 'results' that contains the parameter estimates, the SE's, and the t-stats;
data results; merge parm B;
	tstat=estimate/stderr;
run;