% Calculate Optimal Bandwidths for PSS estimators using Bootstrapping
% Written by Wonho Song, October 2002
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com

function [s1,s2w,s2g,s3]=bootpss(y,x,n,NB,gr_pss,pss_out,p1,EST)

s1=0; s2w=0; s2g=0; s3=0;

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

sgrid=gr_pss(1):gr_pss(2):gr_pss(3);   % grid search range
[sr,sl]=size(sgrid);                   % sl is the number of s loops


%****** start of grid search ******

crits=[]; table=[]; bws=[]; Label=['bandw   '];

ESTSUM=sum(EST(3:5));

if ESTSUM > 0;
disp(' ')
disp('Calculating optimal bandwidths for PSS estimators')
disp(' ')
end;

for indsg=1:sl; table=[table;sprintf('%5.3f',sgrid(indsg))]; end;


if EST(3)==1;
        
	Label =[Label '   PSS1      '];

        crit1=[]; crs1=[];
	for indsg=1:sl       
	s=sgrid(indsg);
	temp1=bpss1(y,x,p1,n,s,NB);
	crit1=[crit1; temp1];
	crs1=[crs1;sprintf('%13.8f',temp1)];
	end;

	[temp,I]=sort(crit1);
	s1= sgrid(I(1));

	% crits=[crits crit1];

	bws=[bws;sprintf('PSS1 = %5.3f',s1)];
	table=[table crs1];
end;

if EST(4)==1;

        Label=[Label '   PSS2W    ' '    PSS2G     '];

	crit2=[]; crs2=[];
	for indsg=1:sl       
	s=sgrid(indsg);
	[temp21,temp22]=bpss2(y,x,n,s,NB);
	crit2=[crit2; [temp21 temp22]];
	crs2=[crs2;sprintf('%13.8f%13.8f',[temp21 temp22])];
	end;

        [tp1,tp2]=size(crit2);

	if tp1==1; s2w= sgrid; s2g= sgrid;
	else; 
	[temp,I]=sort(crit2); 
	s2w= sgrid(I(1,1));
	s2g= sgrid(I(1,2));
        end;

	% crits=[crits crit2];

	bws=[bws;sprintf('PSS2W= %5.3f',s2w)];
	bws=[bws;sprintf('PSS2G= %5.3f',s2g)];
	table=[table crs2];
end;

if EST(5)==1;

        Label=[Label '   PSS3  '];

        crit3=[]; crs3=[];
	for indsg=1:sl       
	s=sgrid(indsg);
	temp3=bpss3(y,x,n,s,NB);
	crit3=[crit3; temp3];  
	crs3=[crs3;sprintf('%13.8f',temp3)];
	end;

	[temp,I]=sort(crit3);
	s3= sgrid(I(1));

	% crits=[crits crit3];

	bws=[bws;sprintf('PSS3 = %5.3f',s3)];
	table=[table crs3];
end;

%**************end of s loop********************


if ESTSUM > 0;

if pss_out ==1;

disp(' ')
disp('Choose the bandwidths which have the smallest MSEs for PSS estimators.')
disp(' ')
disp(' ')
disp('          Total Bootstrap MSE for PSS1,2,3')
disp(' ')
disp(Label)
disp(table)
disp(' ')
disp('Optimal bandwidths chosen: ')
disp(bws)
disp(' ')

end;

end;




function [crit1]=bpss1(y,x,p1,n,bn,NB)

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

%%%%%%%%%%% ESTIMATION %%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%         BOOTSTRAP DATA GENERATION 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load pss1;

crit=0;  NB1=0;

sighat=std(residual);


for indb=1:NB
epsilon=normrnd(0,sighat,nt,1);

y = x*b_pss1 + kron(alphapss1,ones(t,1)) +epsilon ;


%%%%%%% RERUN OF PROGRAM %%%%%%%%%%%%%%%%%

[b_pss1b,se_pss1,PSS1,EPSS1,R2p1]=pss1(y,x,p1,n,bn);



CHECK=[b_pss1b];
CHECK=[isnan(CHECK);isinf(CHECK)];
CHECK=sum(CHECK);

if CHECK>0; NB1=NB1+1; continue; end;



%%%%%%%%%%% END OF BOOTSTRAP %%%%%%%%%%%%

diff2=(b_pss1b-b_pss1).^2;
crit=crit+sum(diff2);
end;
%************end of bootloop********************

crit1=crit/(NB-NB1);






function [crit2w,crit2g]=bpss2(y,x,n,bn,NB)

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

%%%%%%%%%%% ESTIMATION %%%%%%%%%%%%%%%%%%%%%

[b_pss2w,temp1,b_pss2g,temp3,temp4,temp5,temp6,temp7,temp8,temp9]=pss2(y,x,n,bn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%         BOOTSTRAP DATA GENERATION 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load pss2;

%  start of Bootstrap 

crit2w=0; crit2g=0; NB1=0;

sigeps1b=sqrt(shat2)/(sqrt(1-rhat^2));

for indb=1:NB

eps1i=normrnd(0,sigeps1b,n,1);
temp=eps1i;
epsil=eps1i;
for jt=2:t
   temp=rhat*temp + normrnd(0,sqrt(shat2),n,1);
   epsil=[epsil  temp]; % each column of epsil is a new time period
end
epsiltran=epsil'; % each column of epsiltran is a firm
epsilon=reshape(epsiltran,nt,1); % epsilon is reshaped as column vector by firms   


%%%%%%% RERUN OF PROGRAM %%%%%%%%%%%%%%%%%

%******** With Within as initial *********

y=x*b_pss2w +kron(alphapss2w,ones(t,1)) + epsilon;

[b_pss2wb,temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,temp9]=pss2(y,x,n,bn);

%******** With GLS as initial *********

y=x*b_pss2g +kron(alphapss2g,ones(t,1)) + epsilon;

[temp1,temp2,b_pss2gb,temp3,temp4,temp5,temp6,temp7,temp8,temp9]=pss2(y,x,n,bn);


CHECK=[b_pss2wb;b_pss2gb];
CHECK=[isnan(CHECK);isinf(CHECK)];
CHECK=sum(CHECK);

if CHECK>0; NB1=NB1+1; continue; end;




%%%%%%%%%%% END OF BOOTSTRAP %%%%%%%%%%%%

diff2w=(b_pss2w-b_pss2wb).^2;
crit2w=crit2w+sum(diff2w);

diff2g=(b_pss2g-b_pss2gb).^2;
crit2g=crit2g+sum(diff2g);

end;
%************end of bootloop********************

crit2w=crit2w/(NB-NB1);

crit2g=crit2g/(NB-NB1);



function [crit3]=bpss3(y,x,n,bn,NB)

[nt,k]=size(x);
r=nt/n;

%%%%%%%%%%% ESTIMATION %%%%%%%%%%%%%%%%%%%%%

[b_pss3,b_pss3l,temp1,temp2,temp22,temp3,temp4]=pss3(y,x,n,bn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%         BOOTSTRAP DATA GENERATION 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load pss3;

resid=yit-[xit yil]*thetahat-kron(alphapss3,ones(r-1,1));
sigboot=std(resid);


crit=0;  NB1=0;
for indb=1:NB

y=[];

epsilon=normrnd(0,sigboot,r,n);

beta0=thetahat(1:k);
gam=thetahat(k+1);

for i=1:n
j1=r*(i-1)+1;j2=r*i;
xx=x(j1:j2,:);

   yi=xx(1,:)*beta0+alphapss3(i)+epsilon(1,i);
       for j=2:r
          yi=[yi;(gam*yi(j-1)+xx(j,:)*beta0+alphapss3(i)+epsilon(j,i))];
       end
   y=[y;yi];
end


%%%%%%% RERUN OF PROGRAM %%%%%%%%%%%%%%%%%

[b_pss3,b_pss3l,temp1,temp2,temp22,temp3,temp4]=pss3(y,x,n,bn);

thetahatb = [b_pss3;b_pss3l]; 


CHECK=[thetahatb];
CHECK=[isnan(CHECK);isinf(CHECK)];
CHECK=sum(CHECK);

if CHECK>0; NB1=NB1+1; continue; end;




%%%%%%%%%%% END OF BOOTSTRAP %%%%%%%%%%%%

diff2=(thetahatb-thetahat).^2;
crit=crit+sum(diff2);
end;
%************end of bootloop********************


crit3=crit/(NB-NB1);


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