% Fixed and Random Effects Estimator
% Originally written in Gauss by Olvar Bergland 
% Modified and Translated into MATLAB by Wonho Song, October 2002
% Updated May 2014
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com

function [bw,sef,sw,ew,br_1,ser_1,sr,er2,hw,FIX,RND,EFIX,ERND,R2w,R2r]=panel(y,x,n);

[nt,k]=size(x);

iota=ones(nt,1);
Piota=iota*inv(iota'*iota)*iota';
%Qiota=eye(nt)-Piota;                          % eye(nt): out of memory problem encountered 
y=y-Piota*y;
x=x-Piota*x; 


t  =  nt/n;          
k1 = k+1;
nk = n*k1;
n1 = n-1;

in = eye(n);                 
it = eye(t);                 
jt = ones(t,1);              
jn = ones(n,1);              

%    fixed-effects model 

dt = it - jt*jt'/t;


yi=[];xi=[];

for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2];
yi=[yi;dt*y(ind,:)];
xi=[xi;dt*x(ind,:)];
end;

xxi = inv(xi'*xi);      
bw  = xxi*xi'*yi;                
ew  = yi - xi*bw;               
sw  = ew'*ew/(nt-n-k);           
sef = sqrt(sw*(diag(xxi)));     

% bw = coefficient
% sef= standard error
% sw = variance of residual


% Calculation of R2 
ey=y-mean(y);                           
er=y-x*bw; 
R2w=1-(er'*er)/(ey'*ey);
R2w=[R2w;1-(1-R2w)*(nt-1)/(nt-k-n)];


%ey=y-mean(y);                           % Eviews uses this method
%er=yi-xi*bw; 
%R2w=1-(er'*er)/(ey'*ey);
%R2w=[R2w;1-(1-R2w)*(nt-1)/(nt-k-n)];


%ey=yi;                                  % STATA uses this method
%er=yi-xi*bw; 
%R2w=1-(er'*er)/(ey'*ey);
%R2w=[R2w;1-(1-R2w)*(nt-1)/(nt-k-n)];



temp=y-x*bw;
temp=temp-mean(temp);

FIX=mean(reshape(temp,t,n))';   % Fixed Effects
FIX=FIX-mean(FIX);

EFIX=exp(FIX-max(FIX));

vfe=sw*xxi;

  

ym=[];xm=[];

for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2];
ym=[ym;jt'*y(ind,:)/t];
xm=[xm;jt'*x(ind,:)/t];
end;

% between regression:   

xm  = [ones(n,1) xm];             
xxi = inv(xm'*xm);      
bm  = xxi*xm'*ym;                
em  = ym - xm*bm;               
sm  = em'*em/(n-k1);             
%se  = sqrt(sm*(diag(xxi)));     


%%%%%%%%%%%% random-effects model %%%%%%%%%%%%%%%%5

xp  = [ones(nt,1) x]; 

vc=3;

if vc==1;

%%%% Greene's method %%%%%%

swg  = ew'*ew/(nt);        
x_tmp=[ones(nt,1) x];
b_tmp=inv(x_tmp'*x_tmp)*x_tmp'*y;
e_tmp=y-x_tmp*b_tmp;
s_pool=e_tmp'*e_tmp/(nt);
sga1=s_pool -swg;
ah1 = 1-sqrt(swg/(swg+t*sga1));

%%%%%%%% STATA method %%%%%%%%

elseif vc==2; 

sga1 = max(0,sm - sw/t);    % sm = between RSS
ah1 = 1-sqrt(sw/(sw+t*sga1)); 

%%%%%%%% Baltagi(1995) method %%%%%%%%

elseif vc==3; 

ui=y-x*bw;

um=[];
for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2];
um=[um;jt'*ui(ind,:)/t];
end;

sg2 = t*um'*um/(n);
sw2 = (ui-kron(um,ones(t,1)))'*(ui-kron(um,ones(t,1)))/(nt-n-k);
sga1 = (sg2-sw2)/t;
ah1 = 1-sqrt(sw2/sg2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end; 
 
%sqrt(sga1)
%sqrt(sw2)
%sqrt(sw)


yv = y  - ah1*kron(ym,jt);
xv = xp - ah1*kron(xm,jt);

xxi = inv(xv'*xv);      
br  = xxi*xv'*yv;                
er  = yv - xv*br;               
sr  = er'*er/(nt-k1);            
ser = sqrt(sr*(diag(xxi)));     


temp=y-xp*br;
temp=temp-mean(temp);

RND=mean(reshape(temp,t,n))';    % Random Effects
RND=RND-mean(RND);
ERND=exp(RND-max(RND));


% Calculation of R2 
er2=y-xp*br;                             % Eviews uses this method  
ey=y-mean(y);
R2r=1-(er2'*er2)/(ey'*ey);
R2r=[R2r;1-(1-R2r)*(nt-1)/(nt-n-k-1)];

%er2 = yi - xi*br(2:k+1);                % STATA uses this method
%ey = yi;
%R2rr=1-(er2'*er2)/(ey'*ey)


br_1 = br(2:k+1);
ser_1= ser(2:k+1);
 
% br = coefficients: X_1=constant: 
% ser= standard error
% sr = Var residuals

vre=sw*xxi;


w=(bw-br(2:k1))'*inv(vfe-vre(2:k1,2:k1))*(bw-br(2:k1)); 
hw=[w 1-chi2cdf(w,k)];

%   k = rows(b);
%   t = b ./ s;                    
%   p = 2*cdftc(abs(t ),n-k);