@ ******************************************************************@
@-----
      Monte Carlo Study on asymptotic equivalence of OLS and GLS 
	in Panel One factor random effects linear model
          updated 5/5/98                                   --------@
@**The data generation process for simulation is similar to the one*@
@**by Phillips and Hansen (1990)                                 ***@
@*******************************************************************@
output file=1a.out reset; /*Output file*/
N=25;
t=25;
"n = ";n;
"t = ";t;
matrho={0.95};
rep1=1;
Do while rep1<=rows(matrho);
"rho=";matrho[rep1,1];	
matk={0.25, 1};
rep2=1;
Do while rep2<=rows(matk);
"k=";matk[rep2,1];
mattri={0,  0.2, 0.4, 0.6, 0.8};
rep3=1;
Do while rep3<=rows(mattri);
"delta=";mattri[rep3,1];
@***********************************************************************@
@ The following section tells us how to construct a mutlivariate random @
@ variables given a covariance matrix.                                  @
@***********************************************************************@
#include kernels.src
#include base.src


@----Start loop--------@
s1=5^13;         /*seed*/
s2=7^9;
s3=2^21;
i=1;
rec=10000;                    /*number of replications*/
LSDV=zeros(rec,1);         /* Initiate some vectors*/
OLS=zeros(rec,1);
LSBG=zeros(rec,1);
GLS=zeros(rec,1);
FD=zeros(rec,1);
GLSCO=zeros(rec,1);
GLSPW=zeros(rec,1);
vit=zeros(N,T);
ttt=zeros(N,T);
dif=zeros(N,1);
c=zeros(T,T);
cc=zeros(T,T);

do while i<=rec;

/*------Create an autoregressive (1) error structure--------------------*/


sigmamu=sqrt(10*mattri[rep3,1]);	/*std. dev. of random effect*/
sigmaer=sqrt(10*(1-mattri[rep3,1]));	/*std. dev. of error*/
mu=sigmamu*rndns(N+1000, 1, s2);	/*random effect normal(0,sigmamu^2)*/
err=sigmaer*rndns(T+1000, N, s3);	/*error normal(0,sigmaer^2)*/     
mui=mu[1001:rows(mu),.].*.ones(1,T);    /*discard first 1000*/
/**"mu=";mu;
"mui=";mui;***/
eit=err[1001:rows(err),.]';             /*discard first 1000*/

g=int(matk[rep2,1]*T); /*integer part of k*T*/
vgj=zeros(g+1,1);

j=1;              /*create vi1*/
Do while j<=N;
gg=g;
Do while gg>=0;     
vgj[gg+1,1]=(matrho[rep1,1]^gg)*err[(1001-gg),j]';      
gg=gg-1;
Endo;
vit[j,1]=sumc(vgj);
j=j+1;
Endo;

tt=2;   /*tt is the time dimension number*/
Do while tt<=T;
jj=1;
Do while jj<=N;
vit[jj,tt]=matrho[rep1,1]*vit[jj,tt-1]+eit[jj,tt];  /*AR(1)*/
jj=jj+1;
Endo;
tt=tt+1;
Endo;

s=1;   /*Create time trend*/
Do while s<=T;
l=1;
Do while l<=N;
ttt[l,s]=s;
l=l+1;
Endo;
s=s+1;
Endo;

beta=1;
y=mui+ttt*beta+vit;	/*y is generated*/

/*within*/
sec1=ttt-reshape(meanc(ttt'), N, T);
sec2=y;
wxx=sumc(sumc(sec1.*sec1));
wxy=sumc(sumc(sec1.*y));
wyy=sumc(sumc(sec2.*sec2));
bw=wxy/wxx;            /*fixed effects estimators*/
LSDV[i,.]=bw[1,.];            /*record the estimated coefficient of x.*/
/* "lsdv = "; bw[1,.]; */


/*OLS*/
txx=sumc(sumc(ttt.*ttt));
txy=sumc(sumc(ttt.*y));
bo=txy/txx;                /*OLS estimator*/
OLS[i,.]=bo[1,.];


/*First Difference*/
j=1;
Do while j<=N;
dif[j,.]=y[j,T]-y[j,1];
j=j+1;
Endo;
fdn=sumc(dif);
fdd=N*(T-1);
bfd=fdn/fdd;             /*first difference estimator*/
FD[i,.]=bfd[1,.];


/*Work on estimating GLS according to Baltagi chapter 5*/

j=2;
Do while j<=T;
c[j,j-1]=-matrho[rep1,1];
j=j+1;
Endo;
dgnl=(sqrt(1-matrho[rep1,1]^2))|ones(T-1,1);
cc=diagrv(c,dgnl);   /*Prais-Winsten transformation matrix for AR(1) error*/

ystr=(eye(N).*.cc)*reshape(y,N*T,1);
xstr=(eye(N).*.cc)*reshape(ttt,N*T,1); /*PW transformed observations*/

alpha=sqrt((1+matrho[rep1,1])/(1-matrho[rep1,1]));
lta=alpha|ones(T-1,1);
dd=alpha^2+(T-1);
jtabar=(lta*lta')/dd;
sigsqa=dd*(sigmamu^2)*((1-matrho[rep1,1])^2)+(sigmaer^2);
thetaa=1-(sigmaer/sqrt(sigsqa));

/*this matrix is too big to work with!!*/
/**pp=eye(N).*.eye(T)-thetaa*(eye(N).*.jtabar); /*eqn 5.13 in Baltagi p83*/
   ystrstr=pp*ystr;
  xstrstr=pp*xstr;   /*now run ols on the transformed observations*/**/

yitstr=reshape(ystr,N,T);
xitstr=reshape(xstr,N,T);
/*need to program equation 5.14 from Baltagi p83*/
bi=zeros(N,1);   /*initialize vectors*/
bix=zeros(N,1);
yitstst=zeros(N,T);
xitstst=zeros(N,T);
j=1;
Do while j<=N;
bi[j,1]=((alpha-1)*yitstr[j,1]+sumc(yitstr[j,.]'))/dd;
bix[j,1]=((alpha-1)*xitstr[j,1]+sumc(xitstr[j,.]'))/dd;
j=j+1;
Endo;

j=1;
Do while j<=N;
yitstst[j,1]=(yitstr[j,1]-thetaa*alpha*bi[j,1]);
xitstst[j,1]=(xitstr[j,1]-thetaa*alpha*bix[j,1]);
tt=2;
Do while tt<=T;
yitstst[j,tt]=(yitstr[j,tt]-thetaa*bi[j,1]);
xitstst[j,tt]=(xitstr[j,tt]-thetaa*bix[j,1]);
tt=tt+1;
Endo;
j=j+1;
Endo;     /*eqn(5.14) Baltagi p83*/

/*gls*/
gxx=sumc(sumc(xitstst.*xitstst));
gxy=sumc(sumc(xitstst.*yitstst));
bg=gxy/gxx;
GLS[i,.]=bg[1,.];
/*"gls="; bg[1,.];*/

/*Feasible GLS estimators*/
uhat=zeros(N,T);
j=1;
Do while j<=N;
tt=1;
Do while tt<=T;
uhat[j,tt]=y[j,tt]-bw*ttt[j,tt]; /*estimated residual from FE*/
tt=tt+1;
Endo;
j=j+1;
Endo;

uhat1=uhat';
num2=sumc(sumc(uhat1[2:rows(uhat1),.]'.*uhat1[1:rows(uhat1)-1,.]'));
dnm2=sumc(sumc(uhat1[1:rows(uhat1)-1,.]'.*uhat1[1:rows(uhat1)-1,.]'));
rhohat=num2/dnm2;  /*estimate of rho*/
if rhohat>=1;    /*truncate rhohat*/
	 rhohat=.99;
endif;
/*"rhohat=";rhohat;*/


/*get consistent estimates of the variance components(Baltagi and Li,1991)*/
j=2;
Do while j<=T;
c[j,j-1]=-rhohat;
j=j+1;
Endo;

dgnl=(sqrt(1-rhohat^2))|ones(T-1,1);
cc=diagrv(c,dgnl);  /*Prais-Winsten transformation matrix for AR(1) error*/

ystrhat=(eye(N).*.cc)*reshape(y,N*T,1);
xstrhat=(eye(N).*.cc)*reshape(ttt,N*T,1); /*PW transformed observations*/
ystht=reshape(ystrhat,N,T);
xstht=reshape(xstrhat,N,T);

alphahat=sqrt((1+rhohat)/(1-rhohat));
eta=alphahat|ones(T-1,1);
ddhat=alphahat^2+(T-1);
uhatstr=(eye(N).*.cc)*reshape(uhat,N*T,1);
uhtst=reshape(uhatstr,N,T);

ci=zeros(N,1);
uistr1=zeros(N,1);
uistrt=zeros(T-1,N);
ut=uhtst';
j=1;
Do while j<=N;
ci[j,1]=(alphahat*uhtst[j,1]+sumc(ut[2:rows(ut),j]))/ddhat;
uistr1[j,1]=uhtst[j,1];
uistrt[.,j]=ut[2:rows(ut),j];
j=j+1;
Endo;

sigsqaht=sumc(ddhat*(ci.*ci))/N;  /*estimate for sigmaalpha^2*/
cit=reshape(ci,T-1,N);
nm1=sumc((uistr1-alphahat*ci).*(uistr1-alphahat*ci));
nm2=sumc(sumc((uistrt'-cit').*(uistrt'-cit')));
sigsqeht=(nm1+nm2)/(N*(T-1));  /*estimate for sigmaeit^2*/
thetahat=1-(sqrt(sigsqeht)/sqrt(sigsqaht));

j=1;
Do while j<=N;
bi[j,1]=((alphahat-1)*ystht[j,1]+sumc(ystht[j,.]'))/ddhat;
bix[j,1]=((alphahat-1)*xstht[j,1]+sumc(xstht[j,.]'))/ddhat;
yitstst[j,1]=(ystht[j,1]-thetahat*alphahat*bi[j,1]);
xitstst[j,1]=(xstht[j,1]-thetahat*alphahat*bix[j,1]);
tt=2;
Do while tt<=T;
yitstst[j,tt]=(ystht[j,tt]-thetahat*bi[j,1]);
xitstst[j,tt]=(xstht[j,tt]-thetahat*bix[j,1]);
tt=tt+1;
Endo;
j=j+1;
Endo;

/*Cochrane-Orcutt*/
xht=xitstst';
yht=yitstst';
coxx=sumc(sumc(xht[2:rows(xht),.]'.*xht[2:rows(xht),.]'));
coxy=sumc(sumc(xht[2:rows(xht),.]'.*yht[2:rows(yht),.]'));
bco=coxy/coxx;
GLSCO[i,.]=bco[1,.];

/*Prais-Winsten*/
pwxx=sumc(sumc(xitstst.*xitstst));
pwxy=sumc(sumc(xitstst.*yitstst));
bpw=pwxy/pwxx;
GLSPW[i,.]=bpw[1,.];

i=i+1;
endo;
@-----------------------------------------------------------------------@
"varmu= "; (sigmamu^2);
"vareit= "; (sigmaer^2);
"meanbiaslsdv= "; meanc(LSDV)-beta;
"meanbiasols= "; meanc(OLS)-beta;
"meanbiasfd= "; meanc(FD)-beta;
"meanbiasgls= "; meanc(GLS)-beta;
"meanbiasglsco="; meanc(GLSCO)-beta;
"meanbiasglspw="; meanc(GLSPW)-beta;
/*"efflsdv= "; (stdc(GLS)^2)/(stdc(LSDV)^2);
"effols= "; (stdc(GLS)^2)/(stdc(OLS)^2);
"efffd= "; (stdc(GLS)^2)/(stdc(FD)^2);
"vargls= "; (stdc(GLS)^2);
"effglsco="; (stdc(GLS)^2)/(stdc(GLSCO)^2);
"effglspw="; (stdc(GLS)^2)/(stdc(GLSPW)^2);*/
"msegls="; ((stdc(GLS))^2)+((meanc(GLS)-beta)^2);

/*compute relative efficiencies using mse*/
"mseefflsdv="; 
(((stdc(GLS))^2)+((meanc(GLS)-beta)^2))/(((stdc(LSDV))^2)+((meanc(LSDV)-beta)^2));
"mseeffols=";
(((stdc(GLS))^2)+((meanc(GLS)-beta)^2))/(((stdc(OLS))^2)+((meanc(OLS)-beta)^2));
"mseefffd="; 
(((stdc(GLS))^2)+((meanc(GLS)-beta)^2))/(((stdc(FD))^2)+((meanc(FD)-beta)^2));
"mseeffglsco=";
(((stdc(GLS))^2)+((meanc(GLS)-beta)^2))/(((stdc(GLSCO))^2)+((meanc(GLSCO)-beta)^2));
"mseeffglspw=";
(((stdc(GLS))^2)+((meanc(GLS)-beta)^2))/(((stdc(GLSPW))^2)+((meanc(GLSPW)-beta)^2));
"***********************************************";
rep3=rep3+1;
Endo;
"***********************************************";
rep2=rep2+1;
Endo;
"***********************************************";
rep1=rep1+1;
Endo;
"***********************************************";
"replications=";rec;
"beta=";beta;

output off;