% Main Procedures for Efficiency.m
% Written by Wonho Song
% Updated Oct 2005
% Updated May 2014
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com


function []=estm(y,x,n,g_para,pss_para,zp,kss_para,EST);  

[nt,p]=size(x);  t=nt/n;

 trunc1 = g_para(1);
 trunc2 = g_para(2);
   tmtr = g_para(3);
out_fig = g_para(4:6);
out_tab = g_para(7:9);
firmeff = g_para(10);

     NB = pss_para(1);
 gr_pss = pss_para(2:4);
pss_out = pss_para(5);
     p1 = pss_para(6);

     LI = kss_para(1);
     LS = kss_para(2);
 gr_kss = kss_para(3:5);



tt=(1:t)';
tm=kron(ones(n,1),tt);    % linear time trend

if tmtr==1; x_tr=[tm x]; p1=p1+1; [nt,p]=size(x_tr); 
else; x_tr=x;
end;


if p1>=p; 
disp(' ')
disp('p1 must be less than the number of variables p.')
disp(' ')
return;
end;

if sum(EST)==0;
disp(' ')
disp('No estimator is chosen. Please select at least one estimator.')
disp(' ')
return;
end;

if trunc1>0&trunc1<1/n;
disp(' ')
disp('trunc1 is too small. Please choose larger value. ')
disp(' ')
return;
end;

if trunc2>0&trunc2<1/n;
disp(' ')
disp('trunc2 is too small. Please choose larger value. ')
disp(' ')
return;
end;


disp(' ')
disp('Please be patient!')
disp('It may take long time to compute estimates if the data size is large.')
disp(' ')


%************************************************************
%    Finding bandwidth for PSS estimators using bootstrap   *
%************************************************************

[s1,s2w,s2g,s3]=bootpss(y,x_tr,n,NB,gr_pss,pss_out,p1,EST);
 

Label=[]; Label1=[]; Label2=[]; COEF=[]; SE=[]; EFFECTS=[]; EFFICIENCY=[]; 
AVGEFF=[]; metacurve=[];  allcurves=[]; lgd=[]; lgdtv=[]; R2all=[]; INDEFF=[];

disp(' ')
disp('Computing... ')
disp(' ')


%************************************************
%       TRADITIONAL PANEL DATA ESTIMATION       *
%************************************************

if EST(1)==1;

Label=[Label '   FIX    ' ];
Label1=[Label1 '   FIX    ' ];
Label2=[Label2 '   FIX    ' ];

[b_with,se_with,sigf,ef,FIX,EFIX,R2w]=panelf(y,x_tr,n);

INDEFF=[INDEFF FIX];

if trunc1>0;
[temp2,J]=sort(EFIX);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EFIX(mmm)=[];
[n2,temp]=size(EFIX);
else; n2=n;
end;


FIX=kron(FIX,ones(t,1));
EFIX=kron(EFIX,ones(t,1));

COEF=[COEF b_with];
SE=[SE se_with];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS FIX];
EFFICIENCY=[EFFICIENCY EFIX];
end;

AVGEFF=[AVGEFF mean(EFIX)];

VFIX=reshape(EFIX,t,n2);      VFIX2=mean(VFIX')';
if firmeff==1; save VFIX; end;


R2all=[R2all R2w];

allcurves=[allcurves VFIX2];
lgd=[lgd; 'FIX  '];

disp('FIX done')

elseif EST(1)==2;

Label=[Label '   FIX   ' '    RND    '];
Label1=[Label1 '   FIX   ' '    RND    '];
Label2=[Label2 '   FIX   ' '    RND    '];

[b_with,se_with,sigf,ef,b_gls,se_gls,sigr,er,hw,FIX,RND,EFIX,ERND,R2w,R2r]=panel(y,x_tr,n);

INDEFF=[INDEFF FIX RND];

if trunc1>0;
[temp2,J]=sort(EFIX);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EFIX(mmm)=[];
[nw2,temp]=size(EFIX);

[temp2,J]=sort(ERND);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
ERND(mmm)=[];
[ng2,temp]=size(ERND);
else; nw2=n; ng2=n;
end;


FIX=kron(FIX,ones(t,1));
RND=kron(RND,ones(t,1));

EFIX=kron(EFIX,ones(t,1));
ERND=kron(ERND,ones(t,1));

COEF=[COEF b_with b_gls];
SE=[SE se_with se_gls];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS FIX RND];
EFFICIENCY=[EFFICIENCY EFIX ERND];
end;

AVGEFF=[AVGEFF mean(EFIX) mean(ERND)];

VFIX=reshape(EFIX,t,nw2);      VFIX2=mean(VFIX')';
VRND=reshape(ERND,t,ng2);      VRND2=mean(VRND')';

if firmeff==1; save VFIX; save VRND; end;

R2all=[R2all R2w R2r];

allcurves=[allcurves VFIX2 VRND2];
lgd=[lgd; 'FIX  ';'RND  '];

disp('FIX and RND done')

end;



%************************************************
%              Hausman-Taylor                   *
%************************************************

if EST(2) ==1;

Label=[Label '   HT     '];
Label1=[Label1 '   HT     '];
Label2=[Label2 '   HT     '];

[b_ht,se_ht,HT,EHT,R2h]=ht(y,x_tr,p1,n);

INDEFF=[INDEFF HT];

if trunc1>0;
[temp2,J]=sort(EHT);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EHT(mmm)=[];
[n2,temp]=size(EHT);
else; n2=n;
end;


HT=kron(HT,ones(t,1));
EHT=kron(EHT,ones(t,1));

COEF=[COEF b_ht];
SE=[SE se_ht];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS HT];
EFFICIENCY=[EFFICIENCY EHT];
end;

AVGEFF=[AVGEFF mean(EHT)];

VHT=reshape(EHT,t,n2);        VHT2=mean(VHT')';

if firmeff==1; save VHT; end;

R2all=[R2all R2h];

allcurves=[allcurves  VHT2];
lgd=[lgd; 'HT   '];

disp('HT done')

end;


%************************************************
%                  PSS1                         *
%************************************************

if EST(3)==1;

Label=[Label '  PSS1    '];
Label1=[Label1 '  PSS1    '];
Label2=[Label2 '  PSS1    '];

[b_pss1,se_pss1,PSS1,EPSS1,R2p1]=pss1(y,x_tr,p1,n,s1);

INDEFF=[INDEFF PSS1];

if trunc1>0;
[temp2,J]=sort(EPSS1);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EPSS1(mmm)=[];
[n2,temp]=size(EPSS1);
else; n2=n;
end;


PSS1=kron(PSS1,ones(t,1));
EPSS1=kron(EPSS1,ones(t,1));

COEF=[COEF b_pss1];
SE=[SE se_pss1];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS PSS1];
EFFICIENCY=[EFFICIENCY EPSS1];
end;

AVGEFF=[AVGEFF mean(EPSS1)];

VPSS1=reshape(EPSS1,t,n2);    VPSS12=mean(VPSS1')';

if firmeff==1; save VPSS1; end;

R2all=[R2all R2p1];

allcurves=[allcurves  VPSS12];
lgd=[lgd; 'PSS1 '];

disp('PSS1 done')

end;


%************************************************
%                   PSS2                        *
%************************************************

if EST(4)==1;

Label=[Label '  PSS2W  ' '   PSS2G  '];
Label1=[Label1 '  PSS2W  ' '   PSS2G  '];
Label2=[Label2 '  PSS2W  ' '   PSS2G  '];

[b_pss2w,se_pss2w,temp1,temp2,PSS2W,temp4,EPSS2W,temp5,R2p2w,temp6]=pss2(y,x_tr,n,s2w);

[temp3,temp,b_pss2g,se_pss2g,temp4,PSS2G,temp5,EPSS2G,temp6,R2p2g]=pss2(y,x_tr,n,s2g);

INDEFF=[INDEFF PSS2W PSS2G];

if trunc1>0;
[temp2,J]=sort(EPSS2W);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EPSS2W(mmm)=[];
[nw2,temp]=size(EPSS2W);

[temp2,J]=sort(EPSS2G);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EPSS2G(mmm)=[];
[ng2,temp]=size(EPSS2G);
else; nw2=n; ng2=n;
end;


PSS2W=kron(PSS2W,ones(t,1));
EPSS2W=kron(EPSS2W,ones(t,1));

PSS2G=kron(PSS2G,ones(t,1));
EPSS2G=kron(EPSS2G,ones(t,1));

COEF=[COEF b_pss2w b_pss2g];
SE=[SE se_pss2w se_pss2g];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS PSS2W PSS2G];
EFFICIENCY=[EFFICIENCY EPSS2W EPSS2G];
end;

AVGEFF=[AVGEFF mean(EPSS2W) mean(EPSS2G)];

VPSS2W=reshape(EPSS2W,t,nw2);  VPSS2W2=mean(VPSS2W')';
VPSS2G=reshape(EPSS2G,t,ng2);  VPSS2G2=mean(VPSS2G')';

if firmeff==1; save VPSS2W; save VPSS2G; end;

R2all=[R2all R2p2w R2p2g];

allcurves=[allcurves  VPSS2W2 VPSS2G2];
lgd=[lgd; 'PSS2W';'PSS2G'];

disp('PSS2 done')

end;


%************************************************
%                     PSS3                      *
%************************************************

if EST(5) ==1;

Label=[Label '   PSS3   '];
Label1=[Label1 '   PSS3   '];
Label2=[Label2 '   PSS3   '];

[b_pss3,b_pss3l,se_pss3,se_pss3l,PSS3,EPSS3,betai,gammai,se_beta,se_gamma,R2p3]=pss3(y,x_tr,n,s3);

INDEFF=[INDEFF PSS3];

if trunc1>0;
[temp2,J]=sort(EPSS3);
mmm=J(ceil([1:trunc1*n n-trunc1*n+1:n]));
EPSS3(mmm)=[];
[n2,temp]=size(EPSS3);
else; n2=n;
end;


PSS3=kron(PSS3,ones(t,1));
EPSS3=kron(EPSS3,ones(t,1));

COEF=[COEF b_pss3];
SE=[SE se_pss3];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS PSS3];
EFFICIENCY=[EFFICIENCY EPSS3];
end;

AVGEFF=[AVGEFF mean(EPSS3)];

VPSS3=reshape(EPSS3,t,n2);    VPSS32=mean(VPSS3')';

if firmeff==1; save VPSS3; end;


R2all=[R2all R2p3];

allcurves=[allcurves  VPSS32];
lgd=[lgd; 'PSS3 '];

disp('PSS3 done')

end;


%***********************************************
%                  CSS                         *
%***********************************************

if EST(6) ==1;

Label=[Label '   CSSW   '];
Label1=[Label1 '   CSSW   '];
Label2=[Label2 '   CSSW   '];

[b_cssw,se_cssw,CSSW,R2cw]=cssw(y,x,n);

INDEFF=[INDEFF mean(CSSW)'];

if trunc1>0;

mmm=[];
for i=1:t
[temp2,J]=sort(CSSW(i,:));
mmm=[mmm J(ceil([1:trunc1*n n-trunc1*n+1:n]))];
end;
I=unique(mmm);  
CSSW(:,I)=[];
[t2,n2]=size(CSSW);
nt2=t2*n2;

else; nt2=nt; n2=n;
end;


ECSSW=[]; 
for i=1:t
ECSSW=[ECSSW;exp(CSSW(i,:)-max(CSSW(i,:)))];
end;

CSSW=reshape(CSSW,nt2,1);
ECSSW=reshape(ECSSW,nt2,1);

if tmtr==1; b_cssw=[0;b_cssw]; se_cssw=[0;se_cssw]; end;

COEF=[COEF b_cssw ];
SE=[SE se_cssw ];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS CSSW];
EFFICIENCY=[EFFICIENCY ECSSW];
end;

AVGEFF=[AVGEFF mean(ECSSW)];

VCSSW=reshape(ECSSW,t,n2);  VCSSW2=mean(VCSSW')';

if firmeff==1; save VCSSW; end;

R2all=[R2all R2cw];


metacurve=[metacurve VCSSW2];
allcurves=[allcurves VCSSW2];
lgd=[lgd; 'CSSW '];
lgdtv=[lgdtv; 'CSSW '];

disp('CSSW done')


elseif EST(6)==2;

Label=[Label '   CSSW   ' '   CSSG   ' ];
Label1=[Label1 '   CSSW   ' '   CSSG   ' ];
Label2=[Label2 '   CSSW   ' '   CSSG   ' ];

[b_cssw,se_cssw,b_cssg,se_cssg,CSSW,CSSG,R2cw,R2cg]=css(y,x,n,zp);   % updated May 2014

INDEFF=[INDEFF mean(CSSW)' mean(CSSG)'];


if trunc1>0;

mmm=[];
for i=1:t
[temp2,J]=sort(CSSW(i,:));
mmm=[mmm J(ceil([1:trunc1*n n-trunc1*n+1:n]))];
end;
I=unique(mmm);  
CSSW(:,I)=[];
[tw2,nw2]=size(CSSW);
ntw2=tw2*nw2;


mmm=[];
for i=1:t
[temp2,J]=sort(CSSG(i,:));
mmm=[mmm J(ceil([1:trunc1*n n-trunc1*n+1:n]))];
end;
I=unique(mmm);  
CSSG(:,I)=[];
[tg2,ng2]=size(CSSG);
ntg2=tg2*ng2;

else; ntw2=nt; ntg2=nt; nw2=n; ng2=n;
end;


ECSSW=[]; ECSSG=[];
for i=1:t
ECSSW=[ECSSW;exp(CSSW(i,:)-max(CSSW(i,:)))];
ECSSG=[ECSSG;exp(CSSG(i,:)-max(CSSG(i,:)))];
end;

CSSW=reshape(CSSW,ntw2,1);
ECSSW=reshape(ECSSW,ntw2,1);

CSSG=reshape(CSSG,ntg2,1);
ECSSG=reshape(ECSSG,ntg2,1);

if tmtr==1; b_cssw=[0;b_cssw]; b_cssg=[0;b_cssg];
            se_cssw=[0;se_cssw]; se_cssg=[0;se_cssg];
end;

COEF=[COEF b_cssw b_cssg];
SE=[SE se_cssw se_cssg];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS CSSW CSSG];
EFFICIENCY=[EFFICIENCY ECSSW ECSSG];
end;

AVGEFF=[AVGEFF mean(ECSSW) mean(ECSSG)];

VCSSW=reshape(ECSSW,t,nw2);  VCSSW2=mean(VCSSW')';
VCSSG=reshape(ECSSG,t,ng2);  VCSSG2=mean(VCSSG')';

if firmeff==1; save VCSSW; save VCSSG; end;

R2all=[R2all R2cw R2cg];

metacurve=[metacurve VCSSW2 VCSSG2];
allcurves=[allcurves VCSSW2 VCSSG2];
lgd=[lgd; 'CSSW ';'CSSG '];
lgdtv=[lgdtv; 'CSSW '; 'CSSG '];

disp('CSSW and CSSG done')

end;


%************************************************
%                     KSS                       *
%************************************************

if EST(7)==1;

Label=[Label '   KSS    '];
Label1=[Label1 '   KSS    '];
Label2=[Label2 '   KSS    '];

[b_kss,se_kss,KSS,L,pp,R2k,zz]=kss(y,x,n,LI,LS,gr_kss);

INDEFF=[INDEFF mean(KSS)'];


if trunc1>0;

mmm=[];
for i=1:t
[temp2,J]=sort(KSS(i,:));
mmm=[mmm J(ceil([1:trunc1*n n-trunc1*n+1:n]))];
end;
I=unique(mmm);  
KSS(:,I)=[];
[t2,n2]=size(KSS);
nt2=t2*n2;

else; nt2=nt; n2=n;
end;

EKSS=[];
for i=1:t
EKSS=[EKSS;exp(KSS(i,:)-max(KSS(i,:)))];
end;

KSS=reshape(KSS,nt2,1); EKSS=reshape(EKSS,nt2,1);

if tmtr==1; b_kss=[0;b_kss]; se_kss=[0;se_kss]; end;

COEF=[COEF b_kss]; 
SE=[SE se_kss];

if trunc1==0&trunc2==0;
EFFECTS=[EFFECTS KSS];
EFFICIENCY=[EFFICIENCY EKSS];
end;

AVGEFF=[AVGEFF mean(EKSS)];

VKSS=reshape(EKSS,t,n2);    VKSS2=mean(VKSS')';

if firmeff==1; save VKSS; end;

R2all=[R2all R2k];


metacurve=[metacurve VKSS2];
allcurves=[allcurves VKSS2];
lgd=[lgd; 'KSS  '];
lgdtv=[lgdtv; 'KSS  '];

disp('KSS done')

end;


%************************************************
%              Battese and Coelli               *
%************************************************

if EST(8)==1;

Label1=[Label1 '    BC    '];
Label2=[Label2 '    BC    '];

[b_bc,se_bc,TE,BCi,Llk]=bc(y,x_tr,n);

if trunc2>0;

mmm=[];
for i=1:t
[temp2,J]=sort(BCi(i,:));
mmm=[mmm J(ceil([1:trunc2*n n-trunc2*n+1:n]))];
end;
I=unique(mmm);  
BCi(:,I)=[];
[t2,n2]=size(BCi);
nt2=t2*n2;

else; nt2=nt; n2=n;
end;


EBC=reshape(BCi,nt2,1);


COEF=[COEF b_bc];
SE=[SE se_bc];

if trunc1==0&trunc2==0;
EFFICIENCY=[EFFICIENCY EBC];
end;

AVGEFF=[AVGEFF mean(EBC)];

VBC=reshape(EBC,t,n2);      VBC2=mean(VBC')';

if firmeff==1; save VBC; end;

Llk=[Llk;0];

R2all=[R2all Llk];

metacurve=[metacurve VBC2];
allcurves=[allcurves VBC2];
lgd=[lgd; 'BC   '];
lgdtv=[lgdtv; 'BC   '];

disp('BC done')

end;


%************************************************
%                    DEA                        *
%************************************************

% x must be level values (positive)
% use single y, not multiple outputs

if EST(9)==1;

Label2=[Label2 '   DEA    '];

[DEA]=dea(exp(y),exp(x),n);   %  DEA uses level values 


if trunc2>0;

mmm=[];
for i=1:t
[temp2,J]=sort(DEA(i,:));
mmm=[mmm J(ceil([1:trunc2*n n-trunc2*n+1:n]))];
end;
I=unique(mmm);  
DEA(:,I)=[];
[t2,n2]=size(DEA);
nt2=t2*n2;

else; nt2=nt; n2=n;
end;


EDEA=reshape(DEA,nt2,1); 

if trunc1==0&trunc2==0;
EFFICIENCY=[EFFICIENCY EDEA];
end;

AVGEFF=[AVGEFF mean(EDEA)];

VDEA=reshape(EDEA,t,n2);    VDEA2=mean(VDEA')';

if firmeff==1; save VDEA; end;


metacurve=[metacurve VDEA2];
allcurves=[allcurves VDEA2];
lgd=[lgd; 'DEA  '];
lgdtv=[lgdtv; 'DEA  '];

disp('DEA done')

end;


%************************************************
%                   Output                      *
%************************************************




diary on;


disp(' ')
disp('    **** PANEL DATA ESTIMATION ****')
disp(' ')
disp('       Parameters Estimates')
disp(' ')
disp(['   ' Label1 ])
disp(COEF)
disp(' ')
disp('       Standard Errors of Parameters Estimates')
disp(' ')
disp(['   ' Label1])
disp(SE)
disp(' ')
disp('       R2 and adR2')
disp(' ')
disp(['   ' Label1])
disp(R2all)
if EST(8)==1;
disp('* For BC, the log likelihood value is provided instead of R-squared.')
end;
disp(' ')
disp('       Average Efficiencies')
disp(' ')
disp(['   ' Label2])
disp(AVGEFF)


if sum(EST)>1&trunc1==0&trunc2==0;

disp('Standard Correlation of Effects')
disp(['   ' Label ])
disp(corrcoef(EFFECTS))

disp('Spearman Rank Correlation of Effects')
disp(['   ' Label ])
disp(spearman(EFFECTS))

disp('Standard Correlation of Efficiencies')
disp(['   ' Label2 ])
disp(corrcoef(EFFICIENCY))

disp('Spearman Rank Correlation of Efficiencies')
disp(['   ' Label2 ])
disp(spearman(EFFICIENCY))

end;

if EST(1)==2;

disp('Hausman-Wu test p-value')
disp(hw)

end;


if EST(5) ==1;

disp('PSS3 estimator for lagged dependent variable ')
disp([b_pss3l se_pss3l])
disp(' ')

disp('Initial estimator for beta (s.e.) used in PSS3')
disp([betai se_beta])

disp('Initial estimator for gamma (s.e.) used in PSS3')
disp([gammai se_gamma])

disp(' ')
end;



xis=(1:1:t)';

%plot(xis,VCSSW); 
%xlabel('time')
%title('Results from CSSW (all firms)')
%print -dmfile figCSSW                          

%plot(xis,VCSSG); 
%xlabel('time')
%title('Results from CSSG (all firms)')
%print -dmfile figCSSG

%plot(xis,VBC); 
%xlabel('time')
%title('Results from BC (all firms)')
%print -dmfile figBC

%plot(xis,VKSS); 
%xlabel('time')
%title('Results from KSS (all firms)')
%print -dmfile figKSS

%plot(xis,VDEA); 
%xlabel('time')
%title('Results from DEA (all firms)')
%print -dmfile figDEA


save metacurve metacurve;


ESTVAR=sum(EST(6:9));

if ESTVAR > 0;

if ESTVAR==1; avgcurve=metacurve; 
else;
avgcurve=mean(metacurve')';
end;

if out_fig(1)==1;
plot(xis,avgcurve); 
xlabel('Time')
ylabel('Efficiency')
xlim([1 t])
%ylim([0 1])
title('Average of time-variant estimators')
end;


if out_fig(2)==1;
figure;
plot(xis,metacurve); 
xlabel('Time')
ylabel('Efficiency')
xlim([1 t])
%ylim([0 1])
legend(lgdtv);
title('Time-variant estimators')
end;


if out_fig(3)==1;
figure;
plot(xis,allcurves); 
xlabel('Time')
xlim([1 t])
%ylim([0 1])
ylabel('Efficiency')
legend(lgd);
title('Averaged efficiencies from each estimator')
end;


% Print Tables

if out_tab(1)==1;
disp('Average of time-variant estimators')
disp([xis avgcurve]); 
end;



if out_tab(2)==1;
disp('Efficiencies from the Time-variant estimators')
temp1= sum(EST(6:9));
[tp1,tp2]=size(Label2);
temp2= Label2(tp2-temp1*10+1:tp2);
disp(['             ' temp2])
disp([xis metacurve]); 
end;


if out_tab(3)==1;
disp('Individual Effects from each estimator')
disp(['             ' Label])
disp([(1:n)' INDEFF])
end;

end;   % end of ESTVAR>0


disp('Notes:')

if tmtr == 1; 
disp('* The first coefficient of FR, HT, PSS123, BC is the one for the time trend.')
end;

disp('* Correlations among the various estimators are calculated only when both trunc1=0 and trunc2=0.')

if firmeff ==1;
disp('* To load efficiencies for each firm, simply type in "load V**" (**=the name of the estimator) in the command line.')
end;

if EST(7)==1;
disp('     L    zz   pp')
disp([ L zz pp] )
end;



if EST(9)==1|EST(5)==1;
disp('* Remember that the DEA and PSS3 estimators do not support the model with multiple outputs.')
end;


diary off;


delete pss1.mat;
delete pss2.mat;
delete pss3.mat;