Hello community,    Currently I am working on a project with economical penal data. I am trying to replicate a paper in which the author runs a regression with industry dummies and standard errors adjusted by a two-dimensional cluster at the firm and year levels. I think I know how to include the industry dummies (fixed effect) into my code. I am using the following:    proc glm data=&dataset; class ff12; model &y = &x ff12 /solution ; run; quit;   Online I have found one source for two-way clustered standard errors which can be found here: https://sites.google.com/site/markshuaima/home/two-way-clustered-standard-errors-and-sas-code And it looks as follows:  %MACRO REG2DSE(y, x, firm, time, multi, dataset, output); proc surveyreg data=&dataset; cluster &firm; model &Y = &X /covb ; ods output covb=firm; ods output FitStatistics=fit; run;quit; proc surveyreg data=&dataset; cluster &time; model &Y = &X /covb ; ods output covb=time; run;quit; %if &multi=1 %then %do; proc surveyreg data=&dataset; cluster &time &firm; model &y = &x /adjrsq covb; ods output covb=both ; ods output parameterestimates=parm; run;quit; data parm; set parm;keep parameter estimate;run; %end; %else %if &multi=0 %then %do; /* >> my attempt proc glm data=&dataset; class ff12; model &y = &x ff12 /solution ; ods output acovest=both ; ods output parameterestimates=parm; run; quit; /*>>original code proc reg data=&dataset; model &y = &x / adjrsq vif tol hcc acov covb ; ods output acovest=both ; ods output parameterestimates=parm; run;quit; */ data both; set both; parameter=Variable; run; data both; set both;drop variable Dependent Model;run; data parm; set parm;parameter=Variable;Estimates=Estimate;keep parameter estimates;run; %end; data parm1; set parm; n=_n_;m=1;keep m n;run; data parm1;set parm1; by m;if last.m;keep n;run; data both; set both; keep intercept &x; run; data firm; set firm; keep intercept &x; run; data time; set time; keep intercept &x; run; data fit1; set fit; parameter=Label1; Estimates=nValue1; if parameter="R-square" then output; run; data fit1; set fit1; n=1; keep parameter Estimates n; run; proc iml;use both;read all var _num_ into Z;print Z;use firm;read all var _num_ into X;print X; use time;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; data results; merge parm B ; tstat=estimates/stderr;n=0;run; data resultsfit; merge results fit1;by n; run; data &output; set resultsfit; drop n; run; %MEND REG2DSE;   I have marked my attempt as well as the original code. Could someone please explain to me how I can combine these two issues?   Any help is highly appreciated!  ... View more