% Kneip, Sickles, and Song Estimator
% Written by Wonho Song, March 2003
% E-mail: whsong@kiep.go.kr, whsong73@hotmail.com

function [beta1p,se_b1p,KSS1P,L1p,pp,R2k]=kss(y,x,n,LS,gr_kss)
[nt,p] = size(x); t=nt/n;


[L1p,pp]=dim1p(y,x,n,LS,gr_kss); 
[beta1p,se_b1p,KSS1P,wt1p]=method1p(y,x,n,L1p,pp);

ek=y-x*beta1p; ek=ek-mean(ek);

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


%save wt wt1p;

%**********************************************************************
%                   Selection of Dimension for KSS1P
%**********************************************************************

function [L1p,pp]=dim1p(y,x,n,LS,gr_kss);
[nt,p] = size(x); t=nt/n; t1 = (1:t)';
gr_st=gr_kss(1); gr_in=gr_kss(2); gr_en=gr_kss(3);

%*************************** 
%          Step1 
%***************************

ytbar=mean(reshape(y,t,n)')';
y_w=y-kron(ones(n,1),ytbar);

xtbar=[]; for i=1:p
xtbar=[xtbar mean(reshape(x(:,i),t,n)')']; end;
x_w=x-kron(ones(n,1),xtbar);

kk=gr_st:gr_in:gr_en; ccc=[];

for pp=kk
    [Zh,Wh,beta,vi,viz]=vfun(y_w,x_w,n,pp);
    bb=sum((vi-viz).^2)/(n*t)/(1-trace(Zh)/t)^2; 
    ccc=[ccc;bb];
end;


[cc2,I]=sort(ccc);
pp=kk(I(1));


%*************************** 
%         Step2 
%***************************

[Zh,Wh,beta,vi,viz]=vfun(y_w,x_w,n,pp);


temp2=vi-viz;
sig2=sum(temp2.^2)/((n-1)*trace(Wh^2));


%*************************** 
%          Step3 
%***************************

temp=reshape(viz,t,n);
SIGMA=temp*temp'/n;

[vve,vva] = eig(SIGMA);
       va = flipud(diag(vva));   % eigenvalues
       ve = fliplr(vve);         % eigenvectors

Cl=[];
for L=1:LS

       ff = sum(va(L+1:t));
       gg = ve(:,1:L);   % number of principals

       Pl = eye(t)-gg*gg';

%*************************** 
%         Step4 
%***************************

up=n*ff-(n-1)*sig2*trace(Zh*Pl*Zh);
dw=sqrt(2*n*sig2^2*trace((Zh*Pl*Zh)^2));
cc=up/dw;

%disp([up dw cc sig2])

Cl=[Cl;cc];

end;   %*** end of LS loop ***

Cl2=(Cl<2.33);

for i=1:LS
if Cl2(i)==1; L1p=i; break; end;
end;

if sum(Cl2)==0; L1p=LS; end;



%************************************************************************* 
%                           METHOD 1P 
%*************************************************************************

function [beta,se_b,KSS1P,wt]=method1p(y,x,n,L,pp);

[nt,p] = size(x); t=nt/n; t1 = (1:t)';

%*************************** 
%         Step1 
%***************************

ytbar=mean(reshape(y,t,n)')';
y_w=y-kron(ones(n,1),ytbar);

xtbar=[]; for i=1:p
xtbar=[xtbar mean(reshape(x(:,i),t,n)')']; end;
x_w=x-kron(ones(n,1),xtbar);

[Zh,Wh,beta,vi,viz]=vfun(y_w,x_w,n,pp);

%*************************** 
%         Step2
%***************************

temp=reshape(viz,t,n);
SIGMA=temp*temp'/n;

[vve,vva] = eig(SIGMA);
       va = flipud(diag(vva));   % eigenvalues
       ve = fliplr(vve);         % eigenvectors
       gt = sqrt(t)*ve(:,1:L);   % number of principals



%*************************** 
%         Step3 
%***************************


Zh2=gt*inv(gt'*gt)*gt';
Wh2=eye(t)-Zh2;

y_til=[]; x_til=[];
for i=1:n
y_til=[y_til;Wh2*y_w(t*(i-1)+1:t*i)];
x_til=[x_til;Wh2*x_w(t*(i-1)+1:t*i,:)];
end;

beta=inv(x_til'*x_til)*x_til'*y_til;

vi=y_w-x_w*beta;

wt=Zh2*(ytbar-xtbar*beta);

%save wt;

%*************************** 
%         Step4 
%***************************

theta=[]; for i=1:n
theta=[theta inv(gt'*gt)*gt'*vi(t*(i-1)+1:t*i)];
end;

%*************************** 
%         Step5 
%***************************

err=y_til-x_til*beta;
bsig2=err'*err/(nt-n-p);
bcov=bsig2*inv(x_til'*x_til);
se_b=sqrt(diag(bcov));


KSS1P=[];
for i=1:n
KSS1P=[KSS1P gt*theta(:,i)+wt];    
end;

KSS1P=KSS1P-mean(mean(KSS1P));



tv=[];
for i=1:n
v_tmp=vi(t*(i-1)+1:t*i);
tht=theta(:,i);
ee=v_tmp-gt*tht;
s2=ee'*ee/(t-L);
thcov=s2*inv(gt'*gt);
se_th=sqrt(diag(thcov));
tv=[tv tht./se_th];
end;



%*************************************
%    Procedure to calculate viz
%*************************************

function [Zh,Wh,beta,vi,viz]=vfun(y_w,x_w,n,pp);
[nt,p] = size(x_w); t=nt/n; t1 = (1:t)';


Zh=reshape(csaps(t1,eye(t),pp,t1),t,t);
Wh=eye(t)-Zh;

y_til=[]; x_til=[];
for i=1:n
y_til=[y_til;Wh*y_w(t*(i-1)+1:t*i)];
x_til=[x_til;Wh*x_w(t*(i-1)+1:t*i,:)];
end;

beta=inv(x_til'*x_til)*x_til'*y_til;

vi=y_w-x_w*beta;

viz=[];
for i=1:n
viz=[viz;Zh*vi(t*(i-1)+1:t*i)];
end;



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