estm.m

Rice University的Robin C. Sickles professor开发的专门用于paneldata model test and estimation 的program

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

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,L1p,pp,R2k]=kss(y,x,n,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)
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

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;

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(9)==1|EST(5)==1;
disp('* Remember that the DEA and PSS3 estimators do not support the model with multiple outputs.')
end;

if EST(8)==1;
disp('* For BC, the log likelihood value is provided instead of R-squared.')
end;

diary off;


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