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


function [b_pss1,se_pss1,PSS1,EPSS1,R2p1]=pss1(y,x,p1,n,h);

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


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

xbar   = reshape(mean(reshape(x,t,n*p)),n,p);
ybar   = mean(reshape(y,t,n))';

xtilde = x - kron(xbar,ones(t,1));
ytilde = y - kron(ybar,ones(t,1));

ttilde = inv(xtilde'*xtilde)*xtilde'*ytilde;

utilde = ytilde - xtilde*ttilde;

stilde2= utilde'*utilde/(n*(t-1));
ttildecv = stilde2*inv(xtilde'*xtilde);
sbar=ybar-xbar*ttilde;

sbarm=  mean(sbar);
ssbar = sbar - sbarm;
sb2 = ssbar'*ssbar/n;
theta = (stilde2/(t*sb2))^(.5);

z = x(:,p1+1:p);
zbar = xbar(:,p1+1:p);

nx = x(:,1:p1);
nxbar= xbar(:,1:p1);

nxtilde=xtilde(:,1:p1);

ztilde = z -kron(zbar,ones(t,1));

% Hausman - Taylor

xqv = [kron(nxbar,ones(t,1)) xtilde];


% x2hat = xqv*inv(xqv'*xqv)*xqv'*z;
                    %  we do not use this line because it does not work 
                    % if the matrix is too large. Instead, we use the following lines.

tmp1=inv(xqv'*xqv);
tmp2=xqv'*z;
x2hat = xqv*tmp1*tmp2;



xx = [nx x2hat ones(n*t,1)];
[temp,xp]=size(xx);

xxbar   = reshape(mean(reshape(xx,t,n*xp)),n,xp);

xstar = xx - kron(xxbar,ones(t,1))+kron((theta*xxbar),ones(t,1));
ystar = y - kron(ybar,ones(t,1)) + kron((theta*ybar),ones(t,1));

bht = inv(xstar'*xstar)*xstar'*ystar;
tglscov = stilde2*inv(xstar'*xstar);


% Get Hausman Parameter estimates 
bht = bht(1:p);

nxbbar=mean(nxbar);
zbbar=mean(zbar);

zstar=[ones(n,1) zbar];

phatz=zstar*inv(zstar'*zstar)*zstar'*nxbar;


nxbtil = nx - kron(nxbar,ones(t,1)) - (kron(nxbar,ones(t,1)) - kron(nxbbar,ones(n*t,1)));
nzbtil = z - kron(zbar,ones(t,1)) - (kron(zbar,ones(t,1)) -  kron(zbbar,ones(n*t,1)));

SBXZ=(nxbar-phatz)'*(nxbar-phatz)/(n);

SBX=SBXZ;
SWX     = nxbtil'*nxbtil/(n);
SWZ     = nzbtil'*nzbtil/(n);
SWXZ    = nzbtil'*nxbtil/(n);
SW      = [[SWX;SWXZ] [SWXZ';SWZ]];
SWINV   =inv(SW);
h1 = h;

kernel = estpdfk(utilde,utilde,h1,1);
ker = estpdfk(utilde,utilde,h1,2)/h1;
fpdf   = ker./ kernel;

fI0     = mean(fpdf.^2);
cat = [sbar zbar];

% Binwidth of the distribution of the effects and regressors 

kernel = estpdfk(cat,cat,h1,1);
ker = estpdfk(sbar,sbar,h1,2)/h1;
kerp = estpdfk(zbar,zbar,h1,1);
kerprime = ker.*kerp;


wpdw   = kerprime ./ kernel;
I0     = mean(wpdw.^2);

Ihat = [[(SWX*fI0 + I0*SBX) fI0*SWXZ'];[fI0*SWXZ zeros(v,v)]];
Ihat(p1+1:p,p1+1:p)=SWZ*fI0;


Ihatinv = inv(Ihat);

correc1 = (nxbtil)'*fpdf;
correc2 = correc1 + nxbbar'*sum(wpdw);
correc3 = [correc2 - nxbar'*wpdw ; ((nzbtil)'*fpdf)];

b_pss1 = bht + inv(Ihat)*correc3/n;
se_pss1=sqrt(diag(Ihatinv/n));


residual=y-x*b_pss1;
residual=residual-mean(residual);

ep1=residual;

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


PSS1=mean(reshape(residual,t,n))'; alphapss1=PSS1;
PSS1=PSS1-mean(PSS1);
EPSS1=exp(PSS1-max(PSS1));

residual=y-x*b_pss1-kron(PSS1,ones(t,1));

save pss1 alphapss1 residual;





% Kernels and kernel estimator procedures 
function [kk]=kerneln(z);

[temp,k]=size(z);
if k==1;
kk=(1/sqrt((2*pi)^k))*exp(-.5*((z.*z)));  
else;
kk=(1/sqrt((2*pi)^k))*exp(-.5*sum((z.*z)')');  
end;



function [kk]= kernpri(z);
kk=-z.*kerneln(z);




function [kk]=kerneln4(z);

[temp,k]=size(z);

kk=prod(((1.5 - .5*(z.*z)).*(1/sqrt((2*pi)^k)).*exp(-.5*(z.*z)))');




function [kk]=kern4pri(z);

kk=(1.5 - .5*(z.*z)) .* (kernpri(z)) - (z.*kerneln(z));

function [kk]=kernelu(z);
kk=prod((abs(z)<(.5))');




function [pdf]=estpdfk(x,x0,h,c);

[n0,temp]=size(x0);
pdf = zeros(n0,1);

for i = 1:n0
psi = ( x-kron(x0(i,:),ones(n0,1)) )/h;


if c==1;
w = kerneln(psi);
elseif c==2;
w = kernpri(psi);
end;


pdf(i,1) = mean(w)/h;
end;