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


function [beta,se_b,KSS,L,pp,R2k,zz]=kss(y,x,n,LI,LS,gr_kss)
[nt,p] = size(x); t=nt/n;

%% Note: 1-pp is equivalent to kappa in the paper. 
%% MATLAB gives (1-pp) coefficient to the penalty part of the spline smoothing.


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


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

xtbar=zeros(t,p); 
for i=1:p
xtbar(:,i)=mean(reshape(x(:,i),t,n)')'; 
end;

x_w=x-kron(ones(n,1),xtbar);


[L,pp]=dim(y_w,x_w,n,LI,LS,gr_kss,ytbar,xtbar); 
[beta,se_b,KSS,wt,zz]=method(y_w,x_w,n,L,pp,ytbar,xtbar);

save wt wt;


ek=y-x*beta; 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)];



%**********************************************************************
%                   Selection of Dimension for KSS
%**********************************************************************

function [L,pp]=dim(y_w,x_w,n,LI,LS,gr_kss,ytbar,xtbar);
[nt,p] = size(x_w); t=nt/n; t1 = (1:t)';
gr_st=gr_kss(1); gr_in=gr_kss(2); gr_en=gr_kss(3);

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

kk=gr_st:gr_in:gr_en; 


Cl=10*ones(LS,2);   % we select L such that Cl is less than 2.33
%disp('Dimensionality Test of KSS Estimator')

for L=LI:LS

%disp([L])

%% 'Leave-one-out' Cross-validation %%
ccc=[];  

for pp=kk                        %%%%%%%  pp = smoothing parameter

sp=zeros(nt,1);
	for i= 1:n
	   y_w2=y_w; x_w2=x_w;
	   curr=1+t*(i-1):t*i;
	   y_w2(curr)=[];
	   x_w2(curr,:)=[];  
	   [beta,gt]=bg2(y_w2,x_w2,n-1,L,pp);
	   vic=y_w(curr)-x_w(curr,:)*beta;
	   sp(curr)=(vic-gt*inv(gt'*gt)*gt'*vic); 
	end;
    bb=sp'*sp;
 
    ccc=[ccc;bb];
end;
[cc2,I]=sort(ccc);
pp=kk(I(1));

%ccc

clear sp;

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

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


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


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

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

	ff = sum(va(L+1:t));
	gg = ve(:,1:L);   % number of principals
	Pl = eye(t)-gg*gg';

       gt = sqrt(t)*ve(:,1:L);   % number of principals

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


y_til=zeros(nt,1); 
tempy=reshape(y_w,t,n);
tempy=Wh2*tempy;
y_til=reshape(tempy,nt,1);

x_til=zeros(nt,p);
for i=1:p
tempx=reshape(x_w(:,i),t,n);
tempx=Wh2*tempx;
x_til(:,i)=reshape(tempx,nt,1);
end;

beta=inv(x_til'*x_til)*x_til'*y_til;
vi=y_w-x_w*beta;

viz=zeros(nt,1);
tempv=reshape(vi,t,n);
tempv=Zh*tempv;
viz=reshape(tempv,nt,1);


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



	%*************************** 
	%         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(L,:)=[cc pp];

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


Cl2=(Cl(:,1)<2.33);

for i=LI:LS
if Cl2(i)==1; L=i; pp=Cl(i,2); return; end;
end;

if sum(Cl2)==0; L=LS; 
disp('please increase LS ') 
end;



%************************************************************************* 
%                           ESTIMATION
%*************************************************************************

function [beta,se_b,KSS,wt,zz]=method(y_w,x_w,n,L,pp,ytbar,xtbar);

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

[beta,gt,vi,theta,wt,y_til,x_til,va]=bg(y_w,x_w,n,L,pp,ytbar,xtbar);

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


KSS=zeros(t,n);
for i=1:n
KSS(:,i)=gt*theta(:,i)+wt;    
end;

KSS=KSS-mean(mean(KSS));


%%  S.e. of theta, the coefficients of factors
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;


%% dimensionality test %%

bigI=ones(t,1); 
[b_with,se_with,sigf,ef,FIX,EFIX,R2w]=panelf(y_w,x_w,n);
Zh=reshape(csaps(t1,eye(t),pp,t1),t,t);

if L==1;

nume=((bigI-gt)'*(bigI-gt))/t-sigf*trace(Zh*(eye(t)-bigI*bigI'/t)*Zh)/va(1)/n;
deno=sigf*sqrt( 2*trace( ( Zh*( eye(t) - bigI*bigI'/t)*Zh )^2 ) )/va(1)/n;
zz = nume/deno;

else; zz=10;
end;

%*******************************************************************
%                       Beta and gt estimation
%*******************************************************************

function [beta,gt,vi,theta,wt,y_til,x_til,va]=bg(y_w,x_w,n,L,pp,ytbar,xtbar);

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

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


[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=zeros(nt,1); 
tempy=reshape(y_w,t,n);
tempy=Wh2*tempy;
y_til=reshape(tempy,nt,1);

x_til=zeros(nt,p);
for i=1:p
tempx=reshape(x_w(:,i),t,n);
tempx=Wh2*tempx;
x_til(:,i)=reshape(tempx,nt,1);
end;

beta=inv(x_til'*x_til)*x_til'*y_til;
vi=y_w-x_w*beta;
wt=Zh*(ytbar-xtbar*beta);


viz=zeros(nt,1);
tempv=reshape(vi,t,n);
tempv=Zh*tempv;
viz=reshape(tempv,nt,1);


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

theta=zeros(L,n); 
theta=inv(gt'*gt)*gt'*reshape(vi,t,n);



%*******************************************************************
%                       Beta and gt estimation II
%*******************************************************************

function [beta,gt]=bg2(y_w,x_w,n,L,pp);

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

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


[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=zeros(nt,1); 
tempy=reshape(y_w,t,n);
tempy=Wh2*tempy;
y_til=reshape(tempy,nt,1);

x_til=zeros(nt,p);
for i=1:p
tempx=reshape(x_w(:,i),t,n);
tempx=Wh2*tempx;
x_til(:,i)=reshape(tempx,nt,1);
end;


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


refine=0;

if refine==1;

vi=y_w-x_w*beta;
viz=zeros(nt,1);
tempv=reshape(vi,t,n);
tempv=Zh*tempv;
viz=reshape(tempv,nt,1);


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

[vve,vva] = eig(SIGMA);
       va = flipud(diag(vva));   % eigenvalues
       ve = fliplr(vve);         % eigenvectors
       ff = sum(va(L+1:t));
       gg = ve(:,1:L);   % number of principals
       Pl = eye(t)-gg*gg';
       gt = sqrt(t)*ve(:,1:L);   % number of principals

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


y_til=zeros(nt,1); 
tempy=reshape(y_w,t,n);
tempy=Wh2*tempy;
y_til=reshape(tempy,nt,1);

x_til=zeros(nt,p);
for i=1:p
tempx=reshape(x_w(:,i),t,n);
tempx=Wh2*tempx;
x_til(:,i)=reshape(tempx,nt,1);
end;

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


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=zeros(nt,1); 
tempy=reshape(y_w,t,n);
tempy=Wh*tempy;
y_til=reshape(tempy,nt,1);

x_til=zeros(nt,p);
for i=1:p
tempx=reshape(x_w(:,i),t,n);
tempx=Wh*tempx;
x_til(:,i)=reshape(tempx,nt,1);
end;

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

vi=y_w-x_w*beta;

viz=zeros(nt,1);
tempv=reshape(vi,t,n);
tempv=Zh*tempv;
viz=reshape(tempv,nt,1);


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