/*** MIXED FIXED AND RANDOM COEFFICIENT MODEL FOR PANEL DATA:  SAS VERSION  
(For either balanced or unbalanced panel data, but no provision here for missing values)

NOTE: This version does not use orthogonalized variables as mentioned in the papers 
below.  Please make appropriate adjustments or email author (below) for newer version 
if required.

Program by Diana Weinhold, November 2001, correspondence to:  d.weinhold@lse.ac.uk  

Source and discussion: 
1. Weinhold,D "A dynamic 'fixed effects' model for heterogeneous panel data,"  
London School of Economics, unpublished manuscript April 1999 
2. Nair-Reichert,U and D. Weinhold, “Causality Tests for 
Cross Country Panels: A look at FDI and Economic Growth in Less Developed Countries”,  
Oxford Bulletin of Economics and Statistics vol. 63 no. 2, 2001.  ***********************/
/******************************************************************************************/

libname now '.';

%global rplace nxs nz nr xlist zlist yvar rhs nrc;


/********  FILL IN NECESSARY VALUES BELOW ******************************/


%let rhs=%str(intercept x1 x2 x3);  /* full list of rhs variables, in order with fixed first and random last */
%let xlist=%str(x2 x3);             /* list of variables with random coefficients, in order as above */
%let zlist=%str(intercept x1);      /* list of variables with fixed coefficients, in order as above */
%let yvar=%str(y);                  /* dependent variable */
%let nxs= 20;                       /* number of cross section units */
%let nz=2;                          /* number of variables with fixed coefficients */
%let nr=2;                           /* number of variables with random coefficients */

%let rplace=%str(3 4); 
/* list where the random variables are in order of all variables (they must be at end).  
For example, if there are 9 variables and the last 4 are to 
have random coefficients, you would put:  rplace=%str(6 7 8 9); **/ 


data one;
infile 'test.dat';    /* ascii input file with data in stacked panel format */
input cc y x1 x2 x3;  /* Must include indicator variable cc=1..N where N is number of cross section units */


/*************************************************************************************/
/********** NO NEED TO CHANGE ANYTHING BETWEEN THIS LINE AND OUTPUT PART BELOW *********/
/**************************************************************************************/


intercept=1;

data two;
set one;

%macro here(ind);
if cc=&ind then do;
file "sepdel.&ind";
put &yvar &rhs;
end;
%mend here;

%macro loop1;
%do i=1 %to &nxs;
%here(&i);
%end;
%mend loop1;

%loop1;

proc datasets;
delete one two ;


%macro collect(ind);
data getdat&ind;
infile "sepdel.&ind";
input &yvar &rhs;
%mend collect;

%macro loop;
%do i=1 %to &nxs;
%collect(&i);
%end;
%mend loop;

%loop;

%macro matrix1;

proc iml;


ntobs=0;

%do ind=1 %to &nxs;

use getdat&ind;
setin getdat&ind;

read all var {&xlist} into xmat&ind ;
read all var {&yvar} into yvar&ind;
read all var {&zlist} into zmat&ind;
read all var {&rhs} into rhsmat&ind;

%let nrc=ncol(xmat&ind);

nt=nrow(yvar&ind);
t1=(rhsmat&ind`*rhsmat&ind);
t2=inv(t1);
theta&ind=t2*rhsmat&ind`*yvar&ind;
uhat&ind=yvar&ind-(rhsmat&ind*theta&ind);
I&ind=I(nt);
k=ncol(rhsmat&ind);
ss1=SSQ(uhat&ind);
usq&ind=ss1/(nt-k);

nt&ind=nrow(yvar&ind);
ntobs=ntobs+nt&ind;

%end;

/** CREATE TERMS ***/

thetabar=0;
delta=0;

%macro last;

%do ind=1 %to &nxs;
thetabar=thetabar+(1/&nxs)*theta&ind;
%end;

%do ind=1 %to &nxs;
delta=delta+((1/(&nxs-1))*(theta&ind-thetabar)*(theta&ind-thetabar)`);
%end;


rows={&rplace};
cols={&rplace};
deltax=delta[rows,cols];


%do ind=1 %to &nxs;
phi&ind=(xmat&ind*deltax*xmat&ind`)+usq&ind*I&ind;
%end;

aa=J(&nrc,&nrc,0);
bb=J(&nrc,&nrc,0);
cc=J(&nrc,1,0);
dd=J(&nrc,1,0);

%do ind=1 %to &nxs;
iphi&ind=inv(phi&ind);
aa=aa+(xmat&ind`*iphi&ind*xmat&ind);
b2=inv(zmat&ind`*iphi&ind*zmat&ind);
b1=(xmat&ind`*iphi&ind*zmat&ind);
b3=zmat&ind`*iphi&ind*xmat&ind;
bb=bb+(b1*b2*b3);
cc=cc+(xmat&ind`*iphi&ind*yvar&ind);
d3=zmat&ind`*iphi&ind*yvar&ind;
dd=dd+(b1*b2*d3);
%end;

%mend last;
%last;

betabar=inv(aa-bb)*(cc-dd);

tnvar=&nr+&nz;
zzv=(&nxs*&nz);

uu=J(ntobs,1,0);
zz=J(ntobs,zzv,0);
rr=J(ntobs,&nr,0);

locit=1;

%do ind=1 %to &nxs;  /*  calculate standard errors of random coefficient estimates*/
z1=inv(zmat&ind`*zmat&ind);
z2=zmat&ind`*(yvar&ind-xmat&ind*betabar);
bzhat&ind=z1*z2;
yhat&ind=zmat&ind*bzhat&ind+xmat&ind*betabar;

loc1=locit;
loc2=locit+nt&ind-1;
row=loc1:loc2;

loc3=(&ind*&nz);
loc4=loc3-&nz+1;
col=loc4:loc3;
colr=1:&nr;
uu[row,1]=yvar&ind-yhat&ind;
zz[row,col]=zmat&ind;
rr[row,colr]=xmat&ind;

locit=locit+nt&ind;
%end;

xx=zz||rr;

ssqt=SSQ(uu);
cte=(&nxs*&nz+&nr);
sighat=(1/(ntobs-cte))*ssqt;
varcoef=sighat*inv(xx`*xx);
sigsq1=DIAG(varcoef);

dum1=J(cte,1,1);
dum2=J(&nr,1,1);

sigsq=sigsq1*dum1;
sigcoef=SQRT(sigsq);

fir=cte-&nr+1;
row=fir:cte;
sigbet=sigcoef[row,];


varb1=DIAG(DELTAX);
varb=varb1*dum2;


/************************ OUTPUT *************************/
/*** output file is called test.out and contains the estimated mean betahats of the 
random coefficients, their standard error (for calc. t-stats) and the estimated 
variance from the MFR ************************************************************/


head='mean(beta)    sterr(beta)    var(beta)';
sp='  ';

filename out 'test.out';

file out;
put head;
do i=1 to nrow(betabar);
put (betabar[i,1]) sp @ ;
put (sigbet[i,1]) sp @ ;
put (varb[i,1]);
end;

closefile out;

quit;

%mend matrix1;

%matrix1;