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

function [b_pss3,b_pss3l,se_pss3,se_pss3l,PSS3,EPSS3,betai,gammai,se_beta,se_gamma,R2p3]=pss3(y,x,n,bn)

AB=1;
cn=0;


%%%%%%% INITIAL IV ESTIMATION  %%%%%%%

[nt,k]=size(x);
r=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; 



xit=[];yit=[];yil=[]; vxit=[]; vyit=[]; vyil=[]; ivy=[];
dy2=[]; dy3=[];y2=[]; y3=[]; x1=[];x2=[];dx1=[];dx2=[];dx3=[];


for i=1:n

j1=r*(i-1)+2;
j2=r*i;

xit=[xit;x(j1:j2,:)];
yit=[yit;y(j1:j2,:)];
yil=[yil;y(j1-1:j2-1,:)];

end;

ivy=[];

for i=1:n

jj1=r*(i-1)+1;
jj2=r*i;

vxit=[vxit;x(jj1+4:jj2,:)];
vyit=[vyit;y(jj1+4:jj2,:)];
vyil=[vyil;y(jj1+3:jj2-1,:)];

temp=y(jj1+3:jj2-1,:);

dy2=y(jj1+2:jj2-2,:)-y(jj1+1:jj2-3,:);
dy3=y(jj1+1:jj2-3,:)-y(jj1:jj2-4,:);
y2=y(jj1+2:jj2-2,:);
y3=y(jj1+1:jj2-3,:);
x1=x(jj1+3:jj2-1,:);
x2=x(jj1+2:jj2-2,:);
dx1=x(jj1+3:jj2-1,:)-x(jj1+2:jj2-2,:);
dx2=x(jj1+2:jj2-2,:)-x(jj1+1:jj2-3,:);
dx3=x(jj1+1:jj2-3,:)-x(jj1:jj2-4,:);

%iv=[dy2 dy3 y2 x1 dx1 dx2];
iv=[dy2 dy3 y2];

ivy=[ivy;iv*inv(iv'*iv)*iv'*temp];
end;

clear temp;


%if MC==200;
%disp('iv=[dy2 x1]')
%end;


%%%% Preliminary Definitions %%%%%%%

ybar=[]; xbar=[];
for i=1:n
j1=r*(i-1)+1;j2=j1+r-1;
ind=[j1:j2];
xbar=[xbar;mean(x(ind,:))];
ybar=[ybar;mean(y(ind,:))];
end;

vxbar=[]; vybar=[]; vylbar=[]; ivybar=[];
for i=1:n
j1=(r-4)*(i-1)+1;j2=j1+(r-4)-1; ind=[j1:j2];
ivybar=[ivybar;mean(ivy(ind,:))];
vxbar=[vxbar;mean(vxit(ind,:))];
vybar=[vybar;mean(vyit(ind,:))];
vylbar=[vylbar;mean(vyil(ind,:))];
end;

xitwith=xit-kron(xbar,ones(r-1,1));
yitwith=yit-kron(ybar,ones(r-1,1));
yilwith=yil-kron(ybar,ones(r-1,1));

vxitwith=vxit-kron(vxbar,ones(r-4,1));
vyitwith=vyit-kron(vybar,ones(r-4,1));
vyilwith=vyil-kron(vylbar,ones(r-4,1));
ivywith=ivy-kron(ivybar,ones(r-4,1));

%%%%%% Within OLS / Arelleno and Bond IV Estimator  estimator %%%%%%%%%%%%%

if AB==1; 
vxxit=[vxit ivy];

[bf,sef,sw,ew,thetatilde,setilde,sighat2,ss2,hw,FIX,RND,EFIX,ERND]=panel(vyit,vxxit,n);


betai=thetatilde(1:k);   % Arelleno-Bond estimator for beta coefficient
gammai=thetatilde(k+1);  % Arelleno-Bond estimator for lagged dependent variable
se_beta=setilde(1:k);   % setilde(1:k) is the standard error for Arelleno-Bond estimator for beta coefficient
se_gamma=setilde(k+1);  % setilde(k+1) is the standard error for Arelleno-Bond estimator for lagged dependent variable


elseif AB==0;

%%%%  Arelleno and Bond IV using OLS %%%%%%%%%
%abx=[vxitwith ivywith];
%abx2=[vxitwith vyilwith];

%thetatilde=inv(abx'*abx)*abx'*vyitwith;   

%betai=thetatilde(1:k);
%gammai=thetatilde(k+1);

ss2=vyitwith-abx2*thetatilde;
sighat2=sum(ss2.^2)/(n*(r-4));

abxwith=inv(abx'*abx);
setilde=sqrt(diag(sighat2*abxwith));

end;

%%%%%%%% Definitions %%%%%%%%%%%%%%

zit = yit-yil*gammai-xit*betai;                    

zbar=[];
for i=1:n
j1=(r-1)*(i-1)+1;j2=j1+(r-1)-1;ind=[j1:j2];
zbar=[zbar;mean(zit(ind,:))];
end;

zitwith=zit-kron(zbar,ones(r-1,1));


lowgam=zeros(r-1,r-1);

gammaij=[];
for i=0:r-2
gammaij=[gammaij;gammai^i];
end;

for i=1:r-1
lowgam(i:r-1,i)=gammaij(1:r-i);
end;

xiw_tilde=[];
xitw=[]; 

for i=1:n
xxl=x(r*(i-1)+1:i*r-1,:);    
xxitw=lowgam*xxl;   
xxiw=sum(xxitw)/r;
xiw_tilde=[xiw_tilde;xxiw ];   
xitw=[xitw;xxitw];
end;

xitw_with=xitw-kron(xiw_tilde,ones(r-1,1));

ziw_tilde=[];

for i=1:n
j1=(r-1)*(i-1)+1;j2=j1+(r-1)-2; ind=[j1:j2];   
zzitw=lowgam(1:r-2,1:r-2)*zit(ind,:);          
zziw=sum(zzitw)/r;     
ziw_tilde=[ziw_tilde;zziw];
end;


ctg=sum(lowgam')';
c_tilde=sum(ctg)/r;

cc_with=ctg-c_tilde*ones(r-1,1);
c_with=kron(ones(n,1),cc_with);

%%%%%%% log_likelihood %%%%%%%%%%%%%

li_star=[];

%%%%%% li_beta (sum over i) %%%%%%

wkp=[];
for i=1:n
wkp=[wkp;wprime(zbar(i),zbar,bn)./wkernel(zbar(i),zbar,bn,cn)];
end;

lib_A=[];
lig_A=[];
lig_B=[];

for i=1:n
j1=(r-1)*(i-1)+1;j2=j1+(r-1)-1;   ind=[j1:j2];
zitt=zitwith(ind,:);
xitt=xit(ind,:);
yill=yil(ind,:);

lib_A=[lib_A;sum(repmat(zitt,1,k).*xitt./sighat2)];

lig_A=[lig_A;sum(zitt.*yill./sighat2)];
lig_B=[lig_B;(c_tilde/(r-1)*sighat2)*sum(zitt.^2)];

end;

xit_btn=xbar-kron(mean(xbar),ones(n,1));

lib_B=repmat(wkp,1,k).*xit_btn;

li_beta=lib_A-lib_B;

%%%%%% li_gamma (sum over i) %%%%%

xiw_btn=xiw_tilde-kron(mean(xiw_tilde),ones(n,1));
lig_C=wkp.*(xiw_btn*betai+ziw_tilde-c_tilde*zbar);

li_gamma=lig_A+lig_B-lig_C;

%%%%%% li_star %%%%%%
li_star=[li_beta li_gamma];


%%%%%%% Partition of Ihat %%%%%%%%%%%%%
S_wtn=xitwith'*xitwith/n;          
S_btn=xit_btn'*xit_btn/n;
Iw_hat=mean(wkp.^2);

ksi=0;
for ki=1:r-1
for kj=1:r-1
for kk=1:min(ki,kj)-1
ksi=ksi+gammai^(abs(ki-kj)+2*kk);
end;
end;
end;


%%%%% Calculation of I11 matrix
I_11=S_wtn/sighat2+Iw_hat*S_btn;


%%%%% Calculation of I12 matrix
I12_A=sum(repmat(xitw*betai,1,k).*xitwith)/n;    
ctg_long=kron(ones(n,1),ctg);

I12_B=sum(repmat(ctg_long,1,k).*xitwith)/n*mean(zbar);

I12_C=Iw_hat*mean(repmat(xiw_btn*betai,1,k).*xit_btn);           

I_12=(I12_A+I12_B)/sighat2+I12_C;

%%%%% Calculation of I22 matrix
I22_A=((xitw_with*betai)'*(xitw_with*betai))/n/sighat2;  
I22_B=2*(sum(repmat(c_with,1,k).*xitw_with)/n)*betai*mean(zbar)/sighat2;
I22_C=sum(c_with.^2.*kron(zbar.^2,ones(r-1,1)))/n/sighat2;
I22_D=(ksi-r*c_tilde)*sighat2/r^2;
I22_E=(xiw_btn*betai)'*(xiw_btn*betai)/n;

I22_F=0;
for ki=1:r-1
for kj=0:ki-1
I22_F=I22_F+gammai^(2*kj);
end;
end;
I22_F=I22_F*(1-1/r);

I22_G=sum(ctg.^2)/r+2*c_tilde^2/(r-1);

%%%%%% This part added on August 9, 2004 %%%%%%%%%
I22_H =(xiw_tilde*betai+c_tilde*zbar)'*(xiw_tilde*betai+c_tilde*zbar)/n/sighat2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


I_22=I22_A+I22_B+I22_C+Iw_hat*(I22_D+I22_E)+I22_F-I22_G + I22_H;

Ihat = [I_11 I_12';I_12 I_22];


%%%%%%% FINAL STEP %%%%%%%%%%%%%
sli=sum(li_star)';
thetahat = thetatilde + (1/n)*inv(Ihat)*sli;

b_pss3=thetahat(1:k);   
b_pss3l=thetahat(k+1);

%se_pss=sqrt(  diag(inv(Ihat)/n)  );
se_pss=sqrt(  abs(diag(inv(Ihat)/n))  );

se_pss3=se_pss(1:k);  
se_pss3l=se_pss(k+1);      % se_pss(k+1) is the standard error of lagged dependent variable
                           % from PSS3 estimator



ep3=yit-[xit yil]*thetahat; ep3=ep3-mean(ep3);

% Calculation of R2 
ey=y-mean(y);
R2p3=1-(ep3'*ep3)/(ey'*ey);
R2p3=[R2p3;1-(1-R2p3)*(nt-1)/(nt-k-n)];




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INDIVIDUAL EFFECTS ALPHA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

residual=yit-[xit yil]*thetatilde;
PSS3=mean(reshape(residual,r-1,n))';    %  PSS3 = n by 1 matrix
alphapss3=PSS3;
PSS3=PSS3-mean(PSS3);

EPSS3=exp(PSS3-max(PSS3));


% save parameters

save pss3 thetahat yit xit yil alphapss3;


function ww=wkernel(z,zbar,bn,cn)

u=z-zbar;

ww=mean(Kernel1(u,bn))+cn;

function kbn=Kernel1(u,bn)

kbn=(1/bn)*Kernel2(u/bn);

function ku=Kernel2(u)
ku=exp(-u).*(1+exp(-u)).^(-2);


function wp=wprime(z,zbar,bn);

u=z-zbar;

wp=mean(Krnl1(u,bn));

function kbn=Krnl1(u,bn)

kbn=Krnl2(u/bn)/bn^2;

function ku=Krnl2(u)
ku=exp(-u).*(-1+exp(-u))./(1+exp(-u)).^3;


%******** End of Program ***************