function [b, se, BTE, BTEci, EF, TE, otherpara,R2k] = bie(y,x,n,dist)
% Originally written by Joe J. Qian (jhqian@sjtu.edu.cn)
% Modified by Liu Junrong (liujunrong_529@hotmail.com)
% Updated Sep 2015 by Wonho Song (whsong73@hotmail.com)

% if    dist == 0, truncated exponential 
%       dist == 1, truncated half normal (gamma)
%       dist == 2, doubly truncated normal

if (nargin < 4) || isempty(dist) % 
    dist = 0;
end

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

%%%%%%%%%% upper and lower bounds for optimset %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% updated by Song, Sep 2015
X_lb = ones(1,p)*(-1);
X_ub = ones(1,p)*(1);
sig_lb = 0;
sig_ub = Inf;
gam_lb = 0;
gam_ub = 1;
mu_lb = -Inf;
mu_ub = Inf; 
B_lb = ones(1,t)*0;
B_ub = ones(1,t)*Inf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if dist == 0;   % truncated exp:            sigma_u, sigma_v, B
       lb = [X_lb sig_lb sig_lb B_lb ]'; 
       ub = [X_ub sig_ub sig_ub B_ub]'; 
elseif dist == 1;   % truncated half normal:    sigma, gamma, B
       lb = [X_lb sig_lb gam_lb B_lb ]'; 
       ub = [X_ub sig_ub gam_ub B_ub]'; 
elseif dist == 2;   % doubly truncated normal:  sigma, gamma, mu, B
       lb = [X_lb sig_lb gam_lb mu_lb B_lb ]'; 
       ub = [X_ub sig_ub gam_ub mu_ub B_ub]';
end;

options = optimset('LargeScale','on',...
    'GradObj','on',...
    'HessUpdate','bfgs',...
    'DerivativeCheck','off',...
    'Diagnostics','off'); %{'bfgs'} | 'dfp' | 'gillmurray' | 'steepdesc'
warning('off','all');
if dist==0 % truncated exp
theta0 = initipanelexpt([y x],t);
[theta,fval,exitflag,output,L,G,H] = fmincon(@ipanelexpt,theta0,[],[],[],[],lb,ub,[],options,[y x]);
setht = full(sqrt(diag(inv(H))));

b = theta(1:p);
se = setht(1:p);
sigma_u = theta(p+1); sigma_v = theta(p+2); B = theta(p+3:p+t+2); Bse = setht(p+3:p+t+2);
otherpara = [sigma_u sigma_v]';
BB = kron(ones(n,1),B);
e = y - x*b;
ustar = -e - sigma_v^2/sigma_u;
u = ustar + sigma_v * (normpdf(-ustar/sigma_v) - normpdf((BB-ustar)/sigma_v)) ./ ...
    (normcdf((BB-ustar)/sigma_v) - normcdf(-ustar/sigma_v)); % TE
elseif dist == 1 % truncated half normal
theta0 = initipanelhalftgam([y x],t);
[theta,fval,exitflag,output,L,G,H] = fmincon(@ipanelhalftgam,theta0,[],[],[],[],lb,ub,[],options,[y x]);
setht = full(sqrt(diag(inv(H))));
b = theta(1:p);
se = setht(1:p);
sigma = theta(p+1); gamma = theta(p+2); B = theta(p+3:p+t+2); Bse = setht(p+3:p+t+2);
lamda = sqrt(gamma/(1-gamma));
sigma_u = sigma/sqrt(1+1/lamda^2); 
sigma_v = sqrt(sigma^2-sigma_u^2);
otherpara = [sigma_u sigma_v]';
BB = kron(ones(n,1),B);
e = y - x*b;
ustar = -e*sigma_u^2/sigma^2;
sigmastar = sigma_u*sigma_v/sigma;
u = ustar + sigmastar*(normpdf(-ustar/sigmastar)-normpdf((BB-ustar)/sigmastar))...
    ./(normcdf((BB-ustar)/sigmastar) - normcdf(-ustar/sigmastar));
elseif dist == 2 % doubly truncated normal
theta0 = initipaneldoublytgam([y x],t);
[theta,fval,exitflag,output,L,G,H] = fmincon(@ipaneldoublytgam,theta0,[],[],[],[],lb,ub,[],options,[y x]);
setht = full(sqrt(diag(inv(H))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cov_BIE = full(inv(H));
Cov = Cov_BIE(1:p,1:p);
%disp('Covariant matrix of method BIE is:');
%disp(Cov);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b = theta(1:p);
se = setht(1:p);
sigma = theta(p+1); gamma = theta(p+2); mu = theta(p+3); B = theta(p+4:p+t+3); Bse = setht(p+4:p+t+3);
lamda = sqrt(gamma/(1-gamma));
sigma_u = sigma/sqrt(1+1/lamda^2); 
sigma_v = sqrt(sigma^2-sigma_u^2);
otherpara = [sigma_u sigma_v]';
BB = kron(ones(n,1),B);
e = y - x*b;
ustar = (mu*sigma_v^2-e*sigma_u^2)/sigma^2;
sigmastar = sigma_u*sigma_v/sigma;
u = ustar + sigmastar*(normpdf(-ustar/sigmastar)-normpdf((BB-ustar)/sigmastar))...
    ./(normcdf((BB-ustar)/sigmastar) - normcdf(-ustar/sigmastar));
end;    

%%%%%%%%%%%%%%Calculation of R2  %%%%%%%%%%%%%
ey = y - mean(y);
res = e - mean(e);
R2k=1-(res'*res)/(ey'*ey);
R2k=[R2k;1-(1-R2k)*(nt-1)/(nt-p-n)];
%disp('the R2 of BIE')
%disp(R2bie)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%Calculation of Efficiency  %%%%%%%%%%%%%
u = offinf(u,t);
u(u==Inf) = median(u);
EF = reshape(u,t,n);                       % modified Sep 2015 by Song
SF = min(EF');                             % modified Sep 2015 by Song
TE = exp(-EF+kron(SF',ones(1,n)));          % modified Sep 2015 by Song
BTE = exp(-(B-SF'));
BTEci = [exp(-(B+1.96*Bse-SF')) exp(-(B-1.96*Bse-SF'))];
EF = u;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function theta0 = initipanelexpt(X,T)
% estimate initial value for ipanelexpt
% theta0:
%   alpha_1~alpha_p
%   sigma_u
%   sigma_v
%   B(1)~B(T)
[n,m] = size(X);
y = X(:,1);
Z = X(:,2:m);
rr = panelre(y,Z,T);
b = rr.b; % alpha_1~alpha_p
e = y - Z*b;
sigma_u = rr.sigmau;
sigma_v = rr.sigmav;
et = (reshape(e,T,n/T))';
B = max(-et);
theta0 = [b' sigma_u sigma_v B]';

function theta0 = initipanelhalftgam(X,T)
% estimate initial value for ipanelexpt
% theta0:
%   alpha_1~alpha_p
%   sigma
%   gamma
%   B(1)~B(T)

[n,m] = size(X);
y = X(:,1);
Z = X(:,2:m);
rr = panelre(y,Z,T);
b = rr.b; % alpha_1~alpha_p
e = y - Z*b;
sigma_u = rr.sigmau;
sigma_v = rr.sigmav;
sigma = sqrt(sigma_u^2 + sigma_v^2);
gamma = sigma_u^2 / sigma^2;
et = (reshape(e,T,n/T))';
B = max(-et);
theta0 = [b' sigma gamma B]';

function theta0 = initipaneldoublytgam(X,T)
% estimate initial value for ipanelexpt
% theta0:
%   alpha_1~alpha_p
%   sigma
%   gamma
%   B(1)~B(T)

[n,m] = size(X);
y = X(:,1);
Z = X(:,2:m);
rr = panelre(y,Z,T);
b = rr.b; % alpha_1~alpha_p
e = y - Z*b;
sigma_u = rr.sigmau;
sigma_v = rr.sigmav;
sigma = sqrt(sigma_u^2 + sigma_v^2);
gamma = sigma_u^2 / sigma^2;
mu = mean(e);
et = (reshape(e,T,n/T))';
B = max(-et);
theta0 = [b' sigma gamma mu B]';

function rr = panelre(y,X,T) 
[m,k] = size(X); % m = N*T, k: number of regressors
N = m/T;
rr_po = ols(y,[ones(m,1) X]); % estimate pooled regression
rr_fe = panelfe(y,X,T); % fixed effect model
s2_e = rr_fe.s2;
s2_u = rr_po.s2 - s2_e;
theta = 1 - sqrt(s2_e/(s2_e + s2_u * T));
G = (eye(T) - theta * ones(T) / T) / sqrt(s2_e);
yg = zeros(m,1);
Xg = zeros(m,k);
for ii = 1:N
    yg((ii-1)*T+1:ii*T) = G*y((ii-1)*T+1:ii*T);
end
for ii = 1:N
    Xg((ii-1)*T+1:ii*T,:) = G*X((ii-1)*T+1:ii*T,:);
end
%yg = kron(eye(N),G)*y;
%Xg = kron(eye(N),G)*X;
r = ols(yg, Xg);
e = y - X*r.b;
u = mean(reshape(e,T,N));
v = e-kron(u',ones(T,1));
rr.b = r.b;
rr.se = r.se;
rr.bvar = r.bvar;
rr.sse = e'*e;
rr.s2 = rr.sse/(m - k);
rr.R2 = 1 - e'*e / (y'*y);
%rr.R2bar = 1-(e'*e/(m-k)) / (y'*(eye(m)-ones(m,m)/m)*y/(m-1));
rr.sigmau = std(u);
rr.sigmav = std(v);

function rr = panelfe(y,X,T);
% estimate fixed effect panel data model
% regressee: y = (y(11) y(12) ... y(1T) ... y(N1) y(N2) ... y(NT))'
% regressor: X

[m,k] = size(X); % m = N*T, k: number of regressors
N = m/T;
Z = [X kron(sparse(eye(N)), ones(T, 1))];
rr = ols(y, Z);

function rr = ols(y, X)
[n,m] = size(X);
b = inv(X'*X)*X'*y;
e = y - X*b;
sse = e'*e;
s2 = sse/(n - m);
bvar = inv(X'*X)*s2;
se = sqrt(diag(bvar));
R2 = 1 - e'*e / (y'*y);
y1 = y-mean(y);
R2bar = 1-(e'*e/(n-m)) / (y1'*y1/(n-1));

rr.b = b;
rr.se = se;
rr.R2 = R2;
rr.R2bar = R2bar;
rr.bvar = bvar;
rr.s2 = s2;

function y = offinf(x,T)
% get rid of "Inf" in x by performing a robust filtering.

m = length(x);
n = m/T;
y = [];
for i = 0:n-1
    z = x(i*T+1:(i+1)*T);
    z(z==Inf) = median(z);
    y = [y;z];
end

function [logl,g,gi] = ipanelhalftgam(theta,P)
% Log likelihood function for Normal-Truncated 
%   Half Normal panel model with independent time effects
% The model: 
%   y(it) = x(it)'alpha + v(it) - u(it), where
%       v(it) \sim N(0,sigma_v^2)
%       u(it) \sim Truncated Half Normal
% Inputs:
% theta: parameters to estimated
%   theta(1:p): model parameters
%   theta(p+1): sigma
%   theta(p+2): gamma
%   theta(p+3:p+T+2): B(1):B(T)
% P: data
%   P(:,1): y=(y(11) y(12) ... y(1T) ... y(N1) y(N2) ... y(NT))'
%   P(:,2:p+1): x
% 
% Joe J. Qian, 2007 June

[n,p] = size(P);
p = p - 1; % number of regressors
T = length(theta)-p-2;
N = n/T;
alpha = theta(1:p);
sigma = theta(p+1); gamma = theta(p+2); B = theta(p+3:p+T+2);

y = P(:,1); x = P(:,2:p+1); 
e = y - x * alpha;

lamda = sqrt(gamma/(1-gamma));
sigma_u = sigma/sqrt(1+1/lamda^2); 

BB = kron(ones(N,1),B);
A1 = BB/sigma_u;
A3 = ( (BB+e)*lamda + BB/lamda ) / sigma;
A4 = e*lamda / sigma;
PhiA1 = normcdf(A1);
PhiA3 = normcdf(A3);
PhiA4 = normcdf(A4);

phiA1 = normpdf(A1);
phiA3 = normpdf(A3);
phiA4 = normpdf(A4);

dA3A4 = PhiA3-PhiA4; 
if min(dA3A4)<=0
    dA3A4(dA3A4<=0) = 1e-8;
end
dA1A2 = PhiA1-1/2; 
if min(dA1A2)<=0
    dA1A2(dA1A2<=0) = 1e-8;
end
    
logl =  sum(-e.^2/(2*sigma^2) + log(dA3A4)) - sum(log(dA1A2)) - n*log(sigma) - n*log(2*pi)/2;
logl = -logl;

if nargout == 1     
    return;
end

me = reshape((reshape(e,T,N))',N,T);
mA3 = ( (kron(ones(N,1),B')+me)*lamda + kron(ones(N,1),B')/lamda ) / sigma;
mA4 = me*lamda / sigma;

d12 = 1 ./dA1A2;
d34 = 1 ./ dA3A4;
Dlam = sigma/(sqrt(1+1/lamda^2)*lamda)^3;
Dsig = 1/sqrt(1+1/lamda^2);

mphiA3 = normpdf(mA3);
mPhiA3 = normcdf(mA3);
mphiA4 = normpdf(mA4);
mPhiA4 = normcdf(mA4);
dmPhiA3mPhiA4 = mPhiA3 - mPhiA4;
if length(dmPhiA3mPhiA4 == 0)>0
    dmPhiA3mPhiA4(dmPhiA3mPhiA4==0) = 1e-8;
end
md34 = 1 ./ (dmPhiA3mPhiA4);

gi_a = kron( e/sigma - d34.*(phiA3 - phiA4)*lamda,ones(1,p)) .* x / sigma;
gi_sig = e.^2/sigma^3 - d34 .* (phiA3.*A3 - phiA4.*e*lamda/sigma)/sigma ...
    + d12.*phiA1.*BB*Dsig/sigma_u^2 - 1/sigma;
gi_lam = d34 .* (phiA3.*((BB+e)/sigma - BB/(lamda^2*sigma)) - phiA4.*e/sigma)...
    + d12.*phiA1.*BB*Dlam/sigma_u^2;
gi_B = md34 .* mphiA3 * (lamda/sigma + 1/(lamda*sigma))...
    - reshape(d12.*phiA1/sigma_u,T,N)';

lamdagamma = 1/(2*lamda) * (1/(1-gamma)^2);
gi_gam = gi_lam*lamdagamma;

g_a = sum(gi_a);
g_sig = sum(gi_sig);
g_gam = sum(gi_gam);
g_B = sum(gi_B);

g = -[g_a g_sig g_gam g_B]';
if (nargout > 2)
gi = -[gi_a gi_sig gi_gam gi_B]';
end

function [logl,g] = ipanelexpt(theta,P)
% Log likelihood function for Normal-Truncated 
%   Exponetial panel model with independent time effects
% The model: 
%   y(it) = x(it)'alpha + v(it) - u(it), where
%       v(it) \sim N(0,sigma_v^2)
%       u(it) \sim Truncated exponential, f(x) = 1/sigma_u*exp(-x/sigma_u)
%           with upper-bound B(t)
% Inputs:
% theta: parameters to estimated
%   theta(1:p): model parameters
%   theta(p+1): sigma_u
%   theta(p+2): sigma_v
%   theta(p+3:p+T+2): B(1):B(T)
% P: data
%   P(:,1): y=(y(11) y(12) ... y(1T) ... y(N1) y(N2) ... y(NT))'
%   P(:,2:p+1): x
% 
% Joe J. Qian, 2006 June

[n,p] = size(P);
p = p - 1; % number of regressors
T = length(theta)-p-2;
N = n/T;
alpha = theta(1:p);
sigu = theta(p+1); sigv = theta(p+2); B = theta(p+3:p+T+2);

y = P(:,1); x = P(:,2:p+1);
e = y - x * alpha;

BB = kron(ones(N,1),B);
A1 = (BB+e)/sigv + sigv/sigu;
A2 = e/sigv + sigv/sigu;
phiA1 = normpdf(A1);
PhiA1 = normcdf(A1);
phiA2 = normpdf(A2);
PhiA2 = normcdf(A2);
dA1A2 = PhiA1 - PhiA2;
if length(dA1A2<=0) > 0
    dA1A2(dA1A2<=0) = 1e-8;
end
logl = -n*log(sigu) - N*sum(log(1-exp(-B/sigu))) + (n/2) * (sigv/sigu)^2 ...
    + sum(e)/sigu + sum(log(dA1A2));
logl = -logl;

if nargout == 1     
    return;
end

% calculate gradient
me = reshape((reshape(e,T,N))',N,T);
mA1 = (kron(ones(N,1),B')+me)/sigv + sigv/sigu;
mA2 = me/sigv + sigv/sigu;
mphiA1 = normpdf(mA1);
mPhiA1 = normcdf(mA1);
mPhiA2 = normcdf(mA2);
dmA1A2 = mPhiA1 - mPhiA2;
if length(dmA1A2==0) > 0
    dmA1A2(dmA1A2==0) = 1e-8;
end
g_a = -sum(x)/sigu + sum(kron((-phiA1+phiA2)./(dA1A2),ones(1,p)) .* x/sigv);
g_sigu = -n/sigu + N*sum(B/sigu^2 .* exp(-B/sigu) ./ (1-exp(-B/sigu))) - 1/sigu^2*sum(e) - n*sigv^2/sigu^3 ...
    - sum((phiA1 - phiA2) ./ (dA1A2))*sigv/sigu^2;
g_sigv = n*sigv/sigu^2 + sum((phiA1 .* (-(BB+e)/sigv^2 + 1/sigu) - phiA2 .* (-e/sigv^2 + 1/sigu)) ./ (dA1A2));
g_B = (-N/sigu*exp(-B/sigu)./(1-exp(-B/sigu)))' + sum((mphiA1/sigv) ./ dmA1A2);
g = -[g_a g_sigu g_sigv g_B]';

function [logl,g,gi] = ipaneldoublytgam(theta,P)
% Log likelihood function for Normal-Doubly Truncated 
%   Normal panel model with independent time effects
% The model: 
%   y(it) = x(it)'alpha + v(it) - u(it), where
%       v(it) \sim N(0,sigma_v^2)
%       u(it) \sim Doubly truncated normal
% Inputs:
% theta: parameters to estimated
%   theta(1:p): model parameters
%   theta(p+1): sigma
%   theta(p+2): gamma
%   theta(p+3): mu
%   theta(p+4:p+T+3): B(1):B(T)
% P: data
%   P(:,1): y=(y(11) y(12) ... y(1T) ... y(N1) y(N2) ... y(NT))'
%   P(:,2:p+1): x
% 
% Joe J. Qian, 2007 June

[n,p] = size(P);
p = p - 1; % number of regressors
T = length(theta)-p-3;
N = n/T;
alpha = theta(1:p);
sigma = theta(p+1); gamma = theta(p+2); mu = theta(p+3); B = theta(p+4:p+T+3);
mu=0;
y = P(:,1); x = P(:,2:p+1); 
e = y - x * alpha;

lamda = sqrt(gamma/(1-gamma));
sigma_u = sigma/sqrt(1+1/lamda^2); 

BB = kron(ones(N,1),B);
A1 = (BB-mu)/sigma_u;
A2 = -mu/sigma_u;
A3 = ( (BB+e)*lamda + (BB-mu)/lamda ) / sigma;
A4 = (e*lamda - mu/lamda) / sigma;
PhiA1 = normcdf(A1);
PhiA2 = normcdf(A2);
PhiA3 = normcdf(A3);
PhiA4 = normcdf(A4);

phiA1 = normpdf(A1);
phiA2 = normpdf(A2);
phiA3 = normpdf(A3);
phiA4 = normpdf(A4);

dA3A4 = PhiA3-PhiA4; 
if min(dA3A4)<=0
    dA3A4(dA3A4<=0) = 1e-8;
end

dA1A2 = PhiA1-PhiA2; 
if min(dA1A2)<=0
    dA1A2(dA1A2<=0) = 1e-8;
end

logl =  sum(-(e+mu).^2/(2*sigma^2) + log(dA3A4)) - sum(log(dA1A2)) - n*log(sigma) - n*log(2*pi)/2;
logl = -logl;
%disp('loglikehood of BIE is:')
%disp(logl)

if nargout == 1     
    return;
end

me = reshape((reshape(e,T,N))',N,T);
mA3 = ( (kron(ones(N,1),B')+me)*lamda + (kron(ones(N,1),B')-mu)/lamda ) / sigma;
mA4 = (me*lamda - mu/lamda) / sigma;

d12 = 1 ./ dA1A2;
d34 = 1 ./ dA3A4;
Dlam = sigma/(sqrt(1+1/lamda^2)*lamda)^3;
Dsig = 1/sqrt(1+1/lamda^2);

mphiA3 = normpdf(mA3);
mPhiA3 = normcdf(mA3);
mphiA4 = normpdf(mA4);
mPhiA4 = normcdf(mA4);
dmPhiA3mPhiA4 = mPhiA3 - mPhiA4;
if length(dmPhiA3mPhiA4 == 0)>0
    dmPhiA3mPhiA4(dmPhiA3mPhiA4==0) = 1e-8;
end
md34 = 1 ./ (dmPhiA3mPhiA4);

gi_a = kron( (e+mu)/sigma - d34.*(phiA3 - phiA4)*lamda,ones(1,p)) .* x / sigma;
gi_sig = (e+mu).^2/sigma^3 - d34 .* (phiA3.*A3/sigma - phiA4.*A4/sigma) ...
    + d12.*( phiA1.*(BB-mu)*Dsig/sigma_u^2 + phiA2*mu*Dsig/sigma_u^2) - 1/sigma;
gi_lam = d34 .* (phiA3.*((BB+e)/sigma - (BB-mu)/(lamda^2*sigma)) - phiA4.*(e/sigma + mu/(sigma*lamda^2)) )...
    + d12.*( phiA1.*(BB-mu)*Dlam/sigma_u^2 + phiA2*mu*Dlam/sigma_u^2 );
gi_mu = -((e+mu)/sigma^2 + d34.*(phiA3 - phiA4)/(sigma*lamda))...
    + d12.*(phiA1 - phiA2)/sigma_u;
gi_B = md34 .* mphiA3 * (lamda/sigma + 1/(lamda*sigma))...
    - reshape(d12.*phiA1/sigma_u,T,N)';

lamdagamma = 1/(2*lamda) * (1/(1-gamma)^2);
gi_gam = gi_lam*lamdagamma;

g_a = sum(gi_a);
g_sig = sum(gi_sig);
g_gam = sum(gi_gam);
g_mu = sum(gi_mu);
g_B = sum(gi_B);

g = -[g_a g_sig g_gam g_mu g_B]';
if (nargout > 2)
gi = -[gi_a gi_sig gi_gam gi_mu gi_B]';
end