/*** MIXED FIXED AND RANDOM COEFFICIENT MODEL APPROACH TO TESTING CAUSALITY
     IN PANEL DATA:  GAUSS VERSION  *************/

/*** Program by Diana Weinhold. Comments/Questions welcome. ***/
    /****** Email: d.weinhold@lse.ac.uk *****/

new;   /** clear memory of previous stuff to max. work space **/

/* The following program generates a lag length of two (2) for all variables*/
/* and is designed for a multivariate analysis with one (1) causal */
/* variable and two (2) addtional controls variables on the RHS. */
/* Any modifications to this must be made manually, but this */
/* is not a complicated process.   */


/*************REPLACE **'s WITH THE FOLLOWING INFORMATION: *************/

/* Note: For a lag of 2 there will be 2 coefficients for each variable*/

r= 6;  /*number of random coefficients to be estimated: include controls*/
f= 3;  /*number of fixed coefficients to be estimated */
N= ** ; /*number of individual cross section units in panel data set */
T= ** ; /*number of time periods (must be a time-balanced panel)*/

/************************* Read In Data ***************************/
       /*** NOTE: Data MUST be in "Stacked Panel" Form ***/

load datall[N*T,**]= **; /*# of columns and name of ascii data file*/


depv={**};  /* indicate location (column #) of dependent variable */
caus={**};  /* indicate location (column #) of causal variable */
con1={**};  /* you must include two control variables*/
con2={**};

/**************************************************************************/
/*Modifications May Be Made To The "Results" Output File At End of Program*/
/**********Otherwise Do Not Change Anything Below This Line.***************/
/**************************************************************************/

yy=submat(datall,0,depv);
xx=submat(datall,0,caus);
zz1=submat(datall,0,con1);
zz2=submat(datall,0,con2);

xi=zeros(T,N);
z1=ones(T-2,1);
z1i=zeros(T,N);
z2i=zeros(T,N);

yi=zeros(T,N);
sigmau=zeros(1,N);
fixd=zeros(T-2,f*N);
ran=zeros(T-2,r*N);

yyall=zeros((T-2)*N,1);
fixall=zeros((T-2)*N,f);
fixhat=zeros((T-2)*N,1);
ranall=zeros((T-2)*N,r);
yydep=zeros(T-2,N);
sighat=0;
thetabar=zeros(f+r,1);
thetai=zeros(f+r,N);
sigmai=zeros(1,N);
gamzisum=zeros((T-2)*N,1);

i=1;
do until i==N+1;

ggi=seqa(T*i-(T-1),1,T);
yi[.,i]=submat(yy,ggi,0);
xi[.,i]=submat(xx,ggi,0);
z1i[.,i]=submat(zz1,ggi,0);
z2i[.,i]=submat(zz2,ggi,0);
sss=seqa(3,1,T-2);
yydep[.,i]=submat(yi[.,i],sss,0);

/* lag each variable twice: modify this for greater lag lengths*/

xlag1=zeros(T,1);
xlag2=zeros(T,1);
ylag1=zeros(T,1);
ylag2=zeros(T,1);
z1lag1=zeros(T,1);
z1lag2=zeros(T,1);
z2lag1=zeros(T,1);
z2lag2=zeros(T,1);


k=3;
do until k==T+1;
m=k-1;
ylag1[k,.]=yi[m,i];
ylag2[k,.]=yi[m-1,i];
xlag1[k,.]=xi[m,i];
xlag2[k,.]=xi[m-1,i];
z1lag1[k,.]=z1i[m,i];
z1lag2[k,.]=z1i[m-1,i];
z2lag1[k,.]=z2i[m,i];
z2lag2[k,.]=z2i[m-1,i];
k=k+1;
endo;

/*eliminate first two observations of sample size*/

ss=seqa(3,1,T-2);
yylag1=submat(ylag1,ss,0);
yylag2=submat(ylag2,ss,0);
xxlag1=submat(xlag1,ss,0);
xxlag2=submat(xlag2,ss,0);
zz1lag1=submat(z1lag1,ss,0);
zz1lag2=submat(z1lag2,ss,0);
zz2lag1=submat(z2lag1,ss,0);
zz2lag2=submat(z2lag2,ss,0);

gi=seqa((T-2)*i-(T-3),1,T-2);
fi=seqa(f*i-(f-1),1,f);

fixd[.,fi]=z1~yylag1~yylag2;
fixall[gi,.]=fixd[.,fi];
yyall[gi,.]=yydep[.,i];

/*************************************/
/**** create orthogonal projection to use instead of xx ***********/

rhs1=z1~yylag1~yylag2~zz1lag1~zz1lag2~zz2lag1~zz2lag2;
param1=inv(rhs1'rhs1)*rhs1'xxlag1;
orlag1=xxlag1-(rhs1*param1);

rhs2=z1~yylag1~yylag2~zz1lag1~zz1lag2~zz2lag1~zz2lag2~orlag1;
param2=inv(rhs2'rhs2)*rhs2'xxlag2;
orlag2=xxlag2-(rhs2*param2);


/**************************************/

/*  if you do not want to use orthogonal projections as the */
/*  indep. vars with random variables, replace the "orlag*" terms */
/*  with "xxlag*" terms.*/

ci=seqa(r*i-(r-1),1,r);
ran[.,ci]=orlag1~orlag2~zz1lag1~zz1lag2~zz2lag1~zz2lag2 ;

ranall[gi,.]=ran[.,ci];

wi=ran[.,ci]~fixd[.,fi];

/*** ESTIMATION of individual OLS regression and begin to construct components
for the model ***/

thetai[.,i] = inv(wi'wi)*wi'yydep[.,i];

eihat=yydep[.,i]-wi*thetai[.,i];
sigmai[1,i]=eihat'eihat/((T-2)-(3+r));

thetabar=thetabar+thetai[.,i]*(1/N);

gi=seqa((T-2)*i-(T-3),1,T-2);

i=i+1;
endo;

/** ESTIMATE covariance matrix of the random coefficients and divide the
matrix into separate values for each coefficient **/

delta=0;

i=1;
do until i==N+1;

delta=delta+((1/(N-1))*(thetai[.,i]-thetabar)*(thetai[.,i]-thetabar)');

i=i+1;
endo;

varseq=seqa(1,1,r);
delta11=delta[varseq,varseq];

varb1=delta11[1,1];
varb2=delta11[2,2];

parta=0;
partb=0;

/** construct components for estimation algorithm **/

i=1;
do until i==N+1;

ci=seqa(r*i-(r-1),1,r);
fi=seqa(f*i-(f-1),1,f);

shi=(ran[.,ci]*delta11*ran[.,ci]')+(sigmai[1,i]*EYE(T-2));
inshi=inv(shi);
inzsz=inv(fixd[.,fi]'*inshi*fixd[.,fi]);
part1i=ran[.,ci]'*inshi*ran[.,ci];
part2i=ran[.,ci]'*inshi*fixd[.,fi]*inzsz*fixd[.,fi]'*inshi*ran[.,ci];
part3i=ran[.,ci]'inshi*yydep[.,i];
part4i=ran[.,ci]'*inshi*fixd[.,fi]*inzsz*fixd[.,fi]'*inshi*yydep[.,i];
parta=parta+(part1i-part2i);
partb=partb+(part3i-part4i);

i=i+1;
endo;

betahat=inv(parta)*partb; /*ESTIMATE of the means of the random coefficients*/

/** recover individual estimates of the random and fixed coefficients **/

idelta11=inv(delta11);
theta1i=zeros(r,N);

i=1;
do until i==N+1;
seqthet=seqa(1,1,r);
theta1i[.,i]=thetai[seqthet,i];

ci=seqa(r*i-(r-1),1,r);
fi=seqa(f*i-(f-1),1,f);

Izzi=inv(fixd[.,fi]'fixd[.,fi]);
terma=ran[.,ci]'ran[.,ci]-ran[.,ci]'fixd[.,fi]*Izzi*fixd[.,fi]'*ran[.,ci];
termb=(1/sigmai[1,i])*terma;
termc=termb+idelta11;
termd=inv(termc);
terme=termb*theta1i[.,i]+idelta11*betahat;
ibetahat=termd*terme; /*estimate of the individual random coefficients*/

term2a=yydep[.,i]-ran[.,ci]*theta1i[.,i];
term2b=fixd[.,fi]'term2a;
gammahat=izzi*term2b; /*estimate of the individual fixed coefficients*/

gi=seqa((T-2)*i-(T-3),1,T-2);
gamzisum[gi,.]=fixd[.,fi]*gammahat;

i=i+1;
endo;

/** estimate error term from mixed fixed and random model **/

yybar=meanc(yyall);
i=1;
tss=0;
do until i==(T-2)*N+1;
tss=tss+(yyall[i,.]-yybar)^2;
i=i+1;
endo;

uhat=yyall-(gamzisum+(ranall*betahat));

/** estimate coefficient standard errors and R-bar squared **/

sigone=uhat'uhat;
sighat=(1/((T-2)*N-(f*N+r)))*sigone;
bottom=tss/((T-2)*N-1);
rbarsq=1-(sighat/bottom);
regrsq=1-(sigone/tss);


xzall=ranall~fixall;
varcoef=sighat*inv(xzall'xzall);
sigsq=DIAG(varcoef);
sigcoef=SQRT(sigsq);

/** estimated standard errors for the estimate of the mean coefficient **/

betahat1=betahat[1,.];
sigbet1=submat(sigcoef,1,1);
betahat2=betahat[2,.];
sigbet2=submat(sigcoef,2,1);

bound1a=sigbet1*(-2*sqrt(N));
bound1b=sigbet1*(2*sqrt(N));
hatsig1=sqrt(varb1);
zscore1a=(bound1a-betahat1)/hatsig1;
zscore1b=(bound1b-betahat1)/hatsig1;

bound2a=sigbet2*(-2*sqrt(N));
bound2b=sigbet2*(2*sqrt(N));
hatsig2=sqrt(varb2);
zscore2a=(bound2a-betahat2)/hatsig2;
zscore2b=(bound2b-betahat2)/hatsig2;

/******************RESULTS*******************************/

/*  The est. means of the random variables are betahat1-betahat2 */
/*  The estimated variances are varb1-varb2 */
/*  The standard errors of the est. means are sigbet1-sigbet2 */
/*  Unadjusted R-squared is "regrsq" */
/*  Adjusted R-squared is "rbarsq"   */
/*  Individual estimates of the fixed coefficients may be recovered */
/*  from the vector "gammahat"  */

/*  Results are displayed simply by naming the desired variable */
/*  Below all the results are displayed and saved in an output */
/*  file called "results.out"   */
/************************************************************************/

output file=zsc12.out on;
zscore1a~zscore1b;
zscore2a~zscore2b;
output off;

output file=res12.out on;
betahat1 betahat2 ;
varb1 varb2  ;
sigbet1 sigbet2 ;
regrsq rbarsq;
output off;

/*************************** END ********************************/