% Fixed and Random Effects Estimator
% Originally written in Gauss by Olvar Bergland 
% Modified and Translated into MATLAB by Wonho Song, October 2002
% E-mail: whsong@kiep.go.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);

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 = kron(in,dt)*y;
%xi = kron(in,dt)*x;

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);
R2w=1-(ew'*ew)/(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 = kron(in,jt')*y/t;          
%xm = kron(in,jt')*x/t;          

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       

xp  = [ones(nt,1) x]; 
sgu = sw;
sga = sm - (sw/t);
tht = sgu/(sgu + t*sga);
ah = 1 - sqrt(sw/(t*sm));


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-k-1);

sga1=s_pool -sw;

ah1 = 1-sqrt(sw/(sw+t*sga1));

%sga1
%sw
%ah1

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

%save rn ah t sw sm;

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));

er2  = y - xp*br;  

% Calculation of R2 
ey=y-mean(y);
R2r=1-(er2'*er2)/(ey'*ey);
R2r=[R2r;1-(1-R2r)*(nt-1)/(nt-k-n)];



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

vre=sgu*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);