% PSS2 Estimator
% Written by Park, Sickles, and Simar
% Updated by Wonho Song, May 2014
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com


function [bhat,se_pss2w,bglshat,se_pss2g,PSS2W,PSS2G,EPSS2W,EPSS2G,R2p2w,R2p2g]=pss2(y,x,n,s)

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


%*****************************************
%*****************************************
%         START THE COMPUTATIONS
%*****************************************
%*****************************************
% calculate \bar x_i and \bar y_i
% to get the WITHIN estimator of \beta
%
xbar=[];
ybar=[];
for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2]';
xbar=[xbar;mean(x(ind,:))];
ybar=[ybar;mean(y(ind))];

%[i mean(y(ind))]

end;
%*****************************************
% get \tilde \beta= btilde
%
xwithin=x-kron(xbar,ones(t,1));
ywithin=y-kron(ybar,ones(t,1));
cpw=inv(xwithin'*xwithin);
btilde=cpw*xwithin'*ywithin;


%*****************************************
% get \tilde \rho = rtilde
% the within residuals are given by uw
%
restilde = y-x*btilde; % the OLS residuals
%***********************************************************************
% build the vector of OLS WITHIN residuals and lagged values by firms
%***********************************************************************
uw=ywithin-xwithin*btilde; % the OLS WITHIN residuals
%
% estimator of Var of beta tilde, within OLS
% not necessary in the Monte-Carlo experiment
%
% swols2=uw'*uw/(nt-p);
% Vbwols=swols2*cpw; 
%
%************************************************
% build the vector of residuals(with btilde) by firms: 
% uwt:   lag=0, no observations lost
% uwt0:  lag=0, but first period lost for each firm.
% uwt1:  lag=1
% uwt02:  lag=0, but 2 first periods lost for each firm.
% uwt2:  lag=2
%
uwt=[];
uwt01=[];
uwt1=[];
uwt02=[];
uwt2=[];
for i=1:n
   j1=t*(i-1)+1;j2=j1+t-1; % current index (it) vectorized, by firms
   ind=[j1:j2];
   ind01=[j1+1:j2]';ind1=[j1:j2-1]';
   ind02=[j1+2:j2];ind2=[j1:j2-2];
   uwt=[uwt restilde(ind)];% uwt will be (t x n) at the end of the loop
   uwt01=[uwt01  restilde(ind01)];% uwt01 will be (t-1 x n) at the end of the loop
   uwt1=[uwt1  restilde(ind1)];% uwt1 will be (t-1 x n) at the end of the loop
   uwt02=[uwt02  restilde(ind02)];% uwt02 will be (t-2 x n) at the end of the loop
   uwt2=[uwt2  restilde(ind2)];% uwt2 will be (t-2 x n) at the end of the loop
end;
C_i0=diag(uwt'*uwt)/t; % this is a (n x 1) matrix
C_i1=diag(uwt01'*uwt1)/(t-1);% this is a (n x 1) matrix
C_i2=diag(uwt02'*uwt2)/(t-2);% this is a (n x 1) matrix
rtilde=sum(C_i1-C_i2)/sum(C_i0-C_i1);
%
%************************************
%************************************
% Get the weights c_t and e_t
%
trho=(1-rtilde^2)+(t-1)*(1-rtilde)^2;
ct=ones(t,1)*((1-rtilde)^2)/trho;
ct(1)=(1-rtilde)/trho;ct(t)=(1-rtilde)/trho;
%
et=ct*trho*(1+rtilde)/((t-1)*(1-rtilde));
ebig=kron(ones(n,1),et);
%
%*************************************************************************
% calculate \tilde w_{i}, \tilde x_{i}, \tilde z_{it}and \tilde\sigma^2
%*************************************************************************
%
wtilde=[];
xtilde=[];
ytilde=[];
for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2]';
temp1=ct'*restilde(ind);
temp2=ct'*x(ind,:);
temp3=ct'*y(ind);
wtilde=[wtilde;temp1];% at the end wtilde is (n x 1)
xtilde=[xtilde;temp2];% at the end xtilde is (n x p)
ytilde=[ytilde;temp3];% at the end ytilde is (n x 1)
end;
ztilde = restilde-kron(wtilde,ones(t,1));
stilde2 =(sqrt(ebig).*ztilde)'*(sqrt(ebig).*ztilde)/n;
%
%
xstar1=[];
ystar1=[];
for i=1:n
   j1=t*(i-1)+1;
   temp=x(j1,:)-xtilde(i,:);
   tempy=y(j1)-ytilde(i);
   xstar1=[xstar1;temp];
   ystar1=[ystar1;tempy];
end;
xstar1=sqrt(1-rtilde^2)*xstar1; % is a (n x p) matrix
ystar1=sqrt(1-rtilde^2)*ystar1; % is a (n x 1) vector
%
xstar=[];
ystar=[];
for i=1:n
   j1=t*(i-1)+1;j2=j1+t-1;
	ind=[j1+1:j2]';ind1=[j1:j2-1]';
   temp=x(ind,:)-rtilde*x(ind1,:)-(1-rtilde)*ones(t-1,1)*xtilde(i,:);
   xstar=[xstar; xstar1(i,:); temp];
   tempy=y(ind)-rtilde*y(ind1)-(1-rtilde)*ones(t-1,1)*ytilde(i);
   ystar=[ystar; ystar1(i); tempy];
end;
% xstar is nt x p stored by firms, and ystar is nt x 1
%
%  GLS within estimator
%
btildegls=inv(xstar'*xstar)*xstar'*ystar;

%****************************************
% calculate \tilde x_{\cdot} = xtilbar
%
xtilbar=mean(xtilde);
xtilcent=xtilde-ones(n,1)*xtilbar; % xtilcent is (n x p) matrix
%
%*********************************************************************
% calculate I_f, \Sigma_1, \Sigma_2, \hat I and  \hat \beta
%*********************************************************************
%
ztilde1=[];
x1=[];
for i=1:n
   j1=t*(i-1)+1;
   ztilde1=[ztilde1;ztilde(j1)];
   x1=[x1;x(j1,:)];
end;
term1=((1-rtilde^2)/stilde2)*(ztilde1'*x1);
%*******************************************
zxrho=[];
for i=1:n
	j1=t*(i-1)+1;j2=j1+t-1;
	ind=[j1+1:j2]';ind1=[j1:j2-1]';
	temp1=ztilde(ind)-rtilde*ztilde(ind1);
   temp2=x(ind,:)-rtilde*x(ind1,:);
   zxrho=[zxrho;temp1'*temp2];
end;
term2=sum(zxrho)/stilde2;
%******************************************************
% ********** density estimator ************************
% with current bandwith s chosen above by the s loop
%******************************************************
%
temp= kron(ones(1,n),exp(wtilde/s));
temp=temp./temp';
%
% now temp_{ij}=exp((\tilde w_i -\tilde w_j)/s)
%
K =(temp./((1+temp).^2))/s;
Kprime = (temp.*(1-temp)./((1+temp).^3))/(s^2);
%
% K is  a (n x n) matrix, so is Kprime (derivatives)
%
% sum over j
K =sum(K'); % now K is a (1 x n) vector: \hat f(w_i),i=1,...,n
%
Kprime =sum(Kprime');% now Kprime is a (1 x n) vector: \hat f^{{1)}(w_i),i=1,...,n
%
wpdw = Kprime./ K; % wpdw is a (1 x n) vector of {\hat f{{1)}\over \hat f}
%
term3=-(wpdw*xtilcent);
%
I_f=mean(wpdw.^2);
%
% Sigma_1 et Sigma_2
%
%
S1 = (xstar'*xstar)/n;
S2 = (xtilcent'*xtilcent)/n;
%
% calculate \hat I and \hat \beta
%
Ihat=S1/stilde2 + S2*I_f;
correc=inv(Ihat)*(term1+term2+term3)'/n;
bhat =btilde +correc;


% Variance of beta gls
Vbwgls=stilde2*inv(S1)/n;
% Variance of beta hat, efficient estimator
Vbhat=inv(Ihat)/n;

% standard error of pss2w
se_pss2w=sqrt(diag(Vbhat));


ep2w=y-x*bhat; ep2w=ep2w-mean(ep2w);

% Calculation of R2 
ey=y-mean(y);
R2p2w=1-(ep2w'*ep2w)/(ey'*ey);
R2p2w=[R2p2w;1-(1-R2p2w)*(nt-1)/(nt-p-n)];


%  
wtilde1=wtilde;
%*******************************************************************************
%*******************************************************************************
% redo the calculations with \hat \beta 
% for better estimating \hat rho, \hat sigma and \hat alpha
%
%*******************************************************************************
%
% get \hat rho= rhat
%---------------------------------------
% the within residuals are given by uw
%
reshat =y-x*bhat ; % the efficient residuals
%
%************************************************
% build the vector of residuals(with bhat) by firms: 
% uwt:   lag=0, no observations lost
% uwt0:  lag=0, but first period lost for each firm.
% uwt1:  lag=1
% uwt02:  lag=0, but 2 first periods lost for each firm.
% uwt2:  lag=2
%
uwt=[];
uwt01=[];
uwt1=[];
uwt02=[];
uwt2=[];
for i=1:n
   j1=t*(i-1)+1;j2=j1+t-1; % current index (it) vectorized, by firms
   ind=[j1:j2];
   ind01=[j1+1:j2]';ind1=[j1:j2-1]';
   ind02=[j1+2:j2];ind2=[j1:j2-2];
   uwt=[uwt reshat(ind)];% uwt will be (t x n) at the end of the loop
   uwt01=[uwt01  reshat(ind01)];% uwt01 will be (t-1 x n) at the end of the loop
   uwt1=[uwt1  reshat(ind1)];% uwt1 will be (t-1 x n) at the end of the loop
   uwt02=[uwt02  reshat(ind02)];% uwt02 will be (t-2 x n) at the end of the loop
   uwt2=[uwt2  reshat(ind2)];% uwt2 will be (t-2 x n) at the end of the loop
end;
C_i0=diag(uwt'*uwt)/t; % this is a (n x 1) matrix
C_i1=diag(uwt01'*uwt1)/(t-1);% this is a (n x 1) matrix
C_i2=diag(uwt02'*uwt2)/(t-2);% this is a (n x 1) matrix
rhat=sum(C_i1-C_i2)/sum(C_i0-C_i1);

if rhat>=1; rhat=0.99; end;    % revised on Sep. 2003
if rhat<=-1; rhat=-0.99; end;

%
%************************************
%************************************
% Get the NEW weights c_t and e_t
%
trho=(1-rhat^2)+(t-1)*(1-rhat)^2;
ct=ones(t,1)*((1-rhat)^2)/trho;
ct(1)=(1-rhat)/trho;
ct(t)=(1-rhat)/trho;
%
et=ct*trho*(1+rhat)/((t-1)*(1-rhat));
ebig=kron(ones(n,1),et);
%
% recalculate \tilde w_{i}, \tilde x_{i}, \tilde z_{it}and \hat \sigma^2
%
wtilde=[];
xtilde=[];
for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2]';
temp1=ct'*reshat(ind);
temp2=ct'*x(ind,:);
wtilde=[wtilde;temp1];
xtilde=[xtilde;temp2];
end;
ztilde = reshat-kron(wtilde,ones(t,1));
shat2 =(sqrt(ebig).*ztilde)'*(sqrt(ebig).*ztilde)/n;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     INDIVIDUAL EFFECTS PSS2W for WITHIN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alphapss2w=wtilde;

PSS2W=wtilde;
PSS2W=PSS2W-mean(PSS2W);

EPSS2W=exp(PSS2W-max(PSS2W));


%***********************************************************
%***********************************************************
% redo the same thing with the GLS within as first step
%***********************************************************
%***********************************************************
% save OLS values of \tilde beta, \tilde rho and \tilde sigma^2
btilde1=btilde;rtilde1=rtilde;stilde21=stilde2;
%
%wtilde1=wtilde;K1=K;restilde1=restilde;Ihat1=Ihat;correc1=correc;
%
btilde=btildegls; % NOW the GLS estimator !!!!!!
%
restilde = y-x*btilde; % NOW the GLS residuals !!!!!!!
%
%************************************************
% build the vector of residuals(with btilde) by firms: 
% uwt:   lag=0, no observations lost
% uwt0:  lag=0, but first period lost for each firm.
% uwt1:  lag=1
% uwt02:  lag=0, but 2 first periods lost for each firm.
% uwt2:  lag=2
%
uwt=[];
uwt01=[];
uwt1=[];
uwt02=[];
uwt2=[];
for i=1:n
   j1=t*(i-1)+1;j2=j1+t-1; % current index (it) vectorized, by firms
   ind=[j1:j2];
   ind01=[j1+1:j2]';ind1=[j1:j2-1]';
   ind02=[j1+2:j2];ind2=[j1:j2-2];
   uwt=[uwt restilde(ind)];% uwt will be (t x n) at the end of the loop
   uwt01=[uwt01  restilde(ind01)];% uwt01 will be (t-1 x n) at the end of the loop
   uwt1=[uwt1  restilde(ind1)];% uwt1 will be (t-1 x n) at the end of the loop
   uwt02=[uwt02  restilde(ind02)];% uwt02 will be (t-2 x n) at the end of the loop
   uwt2=[uwt2  restilde(ind2)];% uwt2 will be (t-2 x n) at the end of the loop
end;
C_i0=diag(uwt'*uwt)/t; % this is a (n x 1) matrix
C_i1=diag(uwt01'*uwt1)/(t-1);% this is a (n x 1) matrix
C_i2=diag(uwt02'*uwt2)/(t-2);% this is a (n x 1) matrix
rtilde=sum(C_i1-C_i2)/sum(C_i0-C_i1);
%
%************************************
%************************************
% Get the weights c_t and e_t
%
trho=(1-rtilde^2)+(t-1)*(1-rtilde)^2;
ct=ones(t,1)*((1-rtilde)^2)/trho;
ct(1)=(1-rtilde)/trho;ct(t)=(1-rtilde)/trho;
%
et=ct*trho*(1+rtilde)/((t-1)*(1-rtilde));
ebig=kron(ones(n,1),et);
%
%*************************************************************************
% calculate \tilde w_{i}, \tilde x_{i}, \tilde z_{it}and \tilde\sigma^2
%*************************************************************************
%
wtilde=[];
xtilde=[];
ytilde=[];
for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2]';
temp1=ct'*restilde(ind);
temp2=ct'*x(ind,:);
temp3=ct'*y(ind);
wtilde=[wtilde;temp1];% at the end wtilde is (n x 1)
xtilde=[xtilde;temp2];% at the end xtilde is (n x p)
ytilde=[ytilde;temp3];% at the end ytilde is (n x 1)
end;
ztilde = restilde-kron(wtilde,ones(t,1));
stilde2 =(sqrt(ebig).*ztilde)'*(sqrt(ebig).*ztilde)/n;
%
%
xstar1=[];
ystar1=[];
for i=1:n
   j1=t*(i-1)+1;
   temp=x(j1,:)-xtilde(i,:);
   tempy=y(j1)-ytilde(i);
   xstar1=[xstar1;temp];
   ystar1=[ystar1;tempy];
end;
xstar1=sqrt(1-rtilde^2)*xstar1; % is a (n x p) matrix
ystar1=sqrt(1-rtilde^2)*ystar1; % is a (n x 1) vector
%
xstar=[];
ystar=[];
for i=1:n
   j1=t*(i-1)+1;j2=j1+t-1;
	ind=[j1+1:j2]';ind1=[j1:j2-1]';
   temp=x(ind,:)-rtilde*x(ind1,:)-(1-rtilde)*ones(t-1,1)*xtilde(i,:);
   xstar=[xstar; xstar1(i,:); temp];
   tempy=y(ind)-rtilde*y(ind1)-(1-rtilde)*ones(t-1,1)*ytilde(i);
   ystar=[ystar; ystar1(i); tempy];
end;
% xstar is nt x p stored by firms, and ystar is nt x 1
%
%****************************************
% calculate \tilde x_{\cdot} = xtilbar
%
xtilbar=mean(xtilde);
xtilcent=xtilde-ones(n,1)*xtilbar; % xtilcent is (n x p) matrix
%
%*********************************************************************
% calculate I_f, \Sigma_1, \Sigma_2, \hat I and  \hat \beta
%*********************************************************************
%
ztilde1=[];
x1=[];
for i=1:n
   j1=t*(i-1)+1;
   ztilde1=[ztilde1;ztilde(j1)];
   x1=[x1;x(j1,:)];
end;
term1=((1-rtilde^2)/stilde2)*(ztilde1'*x1);
%*******************************************
zxrho=[];
for i=1:n
	j1=t*(i-1)+1;j2=j1+t-1;
	ind=[j1+1:j2]';ind1=[j1:j2-1]';
	temp1=ztilde(ind)-rtilde*ztilde(ind1);
   temp2=x(ind,:)-rtilde*x(ind1,:);
   zxrho=[zxrho;temp1'*temp2];
end;
term2=sum(zxrho)/stilde2;
%*****************************************
% ********** density estimator ***********
% with current bandwith chosen above
%
%
temp= kron(ones(1,n),exp(wtilde/s));
temp=temp./temp';
%
% now temp_{ij}=exp((\tilde w_i -\tilde w_j)/s)
%
K =(temp./((1+temp).^2))/s;
Kprime = (temp.*(1-temp)./((1+temp).^3))/(s^2);
%
% K is  a (n x n) matrix, so is Kprime (derivatives)
%
% sum over j
K =sum(K'); % now K is a (1 x n) vector: \hat f(w_i),i=1,...,n
%
Kprime =sum(Kprime');% now Kprime is a (1 x n) vector: \hat f^{{1)}(w_i),i=1,...,n
%
wpdw = Kprime./ K; % wpdw is a (1 x n) vector of {\hat f^{{1)}\over \hat f}
%
term3=-(wpdw*xtilcent);
%
I_f=mean(wpdw.^2);
%
% Sigma_1 et Sigma_2
%
%
S1 = (xstar'*xstar)/n;
S2 = (xtilcent'*xtilcent)/n;
%
% calculate \hat I and \hat \beta
%
Ihat=S1/stilde2 + S2*I_f;
correc=inv(Ihat)*(term1+term2+term3)'/n;
bglshat =btilde +correc;
% Variance of beta gls hat, efficient estimator
Vbglshat=inv(Ihat)/n;   

% standard error of pss2g

se_pss2g=sqrt(diag(Vbglshat));



ep2g=y-x*bglshat; ep2g=ep2g-mean(ep2g);

% Calculation of R2 
ey=y-mean(y);
R2p2g=1-(ep2g'*ep2g)/(ey'*ey);
R2p2g=[R2p2g;1-(1-R2p2g)*(nt-1)/(nt-p-n)];



wtilde1=wtilde;
%*******************************************************************************
%*******************************************************************************
% redo the calculations with \hat \beta 
% for better estimating \hat rho, \hat sigma and \hat alpha
%
%*******************************************************************************
%
% get \hat rho= rhat
%---------------------------------------
% the within residuals are given by uw
%
reshat =y-x*bglshat ; % the efficient residuals
%
%************************************************
% build the vector of residuals(with bhat) by firms: 
% uwt:   lag=0, no observations lost
% uwt0:  lag=0, but first period lost for each firm.
% uwt1:  lag=1
% uwt02:  lag=0, but 2 first periods lost for each firm.
% uwt2:  lag=2
%
uwt=[];
uwt01=[];
uwt1=[];
uwt02=[];
uwt2=[];
for i=1:n
   j1=t*(i-1)+1;j2=j1+t-1; % current index (it) vectorized, by firms
   ind=[j1:j2];
   ind01=[j1+1:j2]';ind1=[j1:j2-1]';
   ind02=[j1+2:j2];ind2=[j1:j2-2];
   uwt=[uwt reshat(ind)];% uwt will be (t x n) at the end of the loop
   uwt01=[uwt01  reshat(ind01)];% uwt01 will be (t-1 x n) at the end of the loop
   uwt1=[uwt1  reshat(ind1)];% uwt1 will be (t-1 x n) at the end of the loop
   uwt02=[uwt02  reshat(ind02)];% uwt02 will be (t-2 x n) at the end of the loop
   uwt2=[uwt2  reshat(ind2)];% uwt2 will be (t-2 x n) at the end of the loop
end;
C_i0=diag(uwt'*uwt)/t; % this is a (n x 1) matrix
C_i1=diag(uwt01'*uwt1)/(t-1);% this is a (n x 1) matrix
C_i2=diag(uwt02'*uwt2)/(t-2);% this is a (n x 1) matrix
rhat=sum(C_i1-C_i2)/sum(C_i0-C_i1);

if rhat>=1; rhat=0.99; end;
if rhat<=-1; rhat=-0.99; end;


%
%************************************
%************************************
% Get the NEW weights c_t and e_t
%
trho=(1-rhat^2)+(t-1)*(1-rhat)^2;
ct=ones(t,1)*((1-rhat)^2)/trho;
ct(1)=(1-rhat)/trho;
ct(t)=(1-rhat)/trho;
%
et=ct*trho*(1+rhat)/((t-1)*(1-rhat));
ebig=kron(ones(n,1),et);
%
% recalculate \tilde w_{i}, \tilde x_{i}, \tilde z_{it}and \hat \sigma^2
%
wtilde=[];
xtilde=[];
for i=1:n
j1=t*(i-1)+1;j2=j1+t-1;
ind=[j1:j2]';
temp1=ct'*reshat(ind);
temp2=ct'*x(ind,:);
wtilde=[wtilde;temp1];
xtilde=[xtilde;temp2];
end;
ztilde = reshat-kron(wtilde,ones(t,1));
shat2 =(sqrt(ebig).*ztilde)'*(sqrt(ebig).*ztilde)/n;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     INDIVIDUAL EFFECTS PSS2G for GLS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alphapss2g=wtilde;

PSS2G=wtilde;
PSS2G=PSS2G-mean(PSS2G);

EPSS2G=exp(PSS2G-max(PSS2G));

save pss2 shat2 rhat alphapss2w alphapss2g;