% Fixed 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,FIX,EFIX,R2w]=panelf(y,x,n);

[nt,k]=size(x);


iota=ones(nt,1);
Piota=iota*inv(iota'*iota)*iota';
%Qiota=eye(nt)-Piota;                               % 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));