Combining fixed effects industry dummies with two way clustered standard errors

[Orange123]

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!