% Battese and Coelli Estimator
% Written by Wonho Song, October 2002
% Updated May 2014
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com

function [b_bc,se_bc,TE,TEi,Llk]=bc(y,x,n)

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

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; 



bini=inv(x'*x)*x'*y; ee=y-x*bini; b0=bini(1); mu=0.0; eta=0.0; m2=mean(ee.^2);

%********************************************************************************* 
%            COLS estimators for b0, sigs2, gamma 
%*********************************************************************************

temp=[];
for i=1:9
gamma=0.1*i;
sigs=sqrt(m2/(1-2*gamma/pi));
bini(1)=b0+sqrt(2*gamma*sigs/pi);
para=[bini;sigs;sqrt(gamma);mu;eta];
temp=[temp;bcfun(para,y,x,n)];
end;

[ms,I]=sort(temp);
gamma=0.1*I(1);
sigs=sqrt(m2/(1-2*gamma/pi));
bini(1)=b0+sqrt(2*gamma*sigs/pi);
para=[bini;sigs;sqrt(gamma);mu;eta];

clear temp;

%************************************************************************************* 
%                                Maximization 
%*************************************************************************************

options=optimset('Display','off','LargeScale','off','HessUpdate','bfgs','GradObj','off');

% select HessUpdate method among bfgs, dfp, steepdesc, gillmurray

[epara,fval,exitflag,output,grad,hessian]=fminunc(@bcfun,para,options,y,x,n);

se=sqrt(diag(inv(hessian)));  % don't have to multiply (-) to hessian 
                              % because we already multiplied (-) to minimize the function
                              % BFGS method approximates HESSIAN!

%se=sqrt(diag((hessian)));  % DFP method approximates INVERSE HESSIAN!

epara=real(epara);

b=epara(1:p);     b_bc=b(2:p); se_bc=se(2:p);

Llk=-fval;   % log likelihood value: 'fval' is the minimum of (-)log likelihood)

%**** Remember that we reparametrized sigs2 and gamma to ensure that it is positive ****
sigs2=epara(p+1)^2;
gamma=epara(p+2)^2;
%*****************************************************************************
mu=epara(p+3);
eta=epara(p+4);


if gamma>=0.99; gamma=0.99; end;
if gamma<=0.01; gamma=0.01; end;
if    mu<=-4;   mu=-3.9; end;

sig2=gamma*sigs2;


%************************************************************** 
%                         Mean TE 
%**************************************************************

t1=1:t; t1=t1';
eta_t=exp(-eta*(t1-t));

temp1=normcdf(-eta_t*sqrt(sig2)+mu/sqrt(sig2))/normcdf(mu/sqrt(sig2));

TE=temp1.*exp(-eta_t*mu+0.5*eta_t.^2*sig2);

%************************************************************* 
%                      Individual TE 
%*************************************************************

eit=y-x*b; 
sv2=sigs2-sig2;
sigi2=sv2*sig2/(sv2+eta_t'*eta_t*sig2);

TEi=[];
for i=1:n
ind=t*(i-1)+1:t*i; ind=ind';
mui=(mu*sv2-eta_t'*eit(ind)*sig2)/(sv2+eta_t'*eta_t*sig2);
temp2=normcdf(-eta_t*sqrt(sigi2)+mui/sqrt(sigi2))/normcdf(mui/sqrt(sigi2));
TEi=[TEi temp2.*exp(-eta_t*mui+0.5*eta_t.^2*sigi2) ];
end;

%************************************************************************* 
%                  log-likelihood function 
%*************************************************************************

function [logL]=bcfun(para,y,x,n)

[nt,p]=size(x); t=nt/n;
t1=1:t; t1=t1';

b=para(1:p); 
%**** Remember that we reparametrized sigs2 and gamma to ensure that it is positive ****
sigs2=para(p+1)^2;
gamma=para(p+2)^2; 
%*****************************************************************************
mu=para(p+3);
eta=para(p+4);


if gamma>=0.99; gamma=0.99; end;
if gamma<=0.01; gamma=0.01; end;
if    mu<=-4;   mu=-3.9; end;

etai=exp(-eta*(t1-t));
z=mu/sqrt(gamma*sigs2);
z=real(z);

if z<-4; z=-4; end;

zistar=[];
for i=1:n
ind=t*(i-1)+1:t*i;ind=ind';

temp1=(mu*(1-gamma)-gamma*etai'*(y(ind)-x(ind,:)*b));
temp2=sqrt(gamma*(1-gamma)*sigs2*(1+(etai'*etai-1)*gamma));
temp=temp1/temp2;
zistar=[zistar;temp];
end;

zistar=real(zistar);
zistar([find(zistar<-4)])=-4;   % find elements smaller than -4 and replace them with -4
%if min(zistar)<-5; min(zistar); end;



ll=(y-x*b)'*(y-x*b)/((1-gamma)*sigs2);

l1=-0.5*nt*(log(2*pi)+log(sigs2)) -0.5*n*(t-1)*log(1-gamma) ;
l2=-0.5*n*log(1+(etai'*etai-1)*gamma)-n*log(normcdf(z))-0.5*n*z^2;
l3=sum(log(normcdf(zistar)))+0.5*sum(zistar.^2);
l4=-0.5*ll;

logL=l1+l2+l3+l4;
logL=-logL;  % we minimize (-) likelihood !!!


if gamma>=0.99; logL=logL+.01*abs(logL); end;
if gamma<=0.01; logL=logL+.01*abs(logL); end;
if    mu<=-4;   logL=logL+.01*abs(logL); end;


%********** End of File **********