% Cornwell-Schmidt Estimator (Within ang GLS)
% Written by Wonho Song, December 2002
% Updated by Wonho Song, May 2014
% E-mail: whsong@cau.ac.kr, whsong73@hotmail.com

function [betaw,se_bw,betag,se_bg,CSSW,CSSG,R2cw,R2cg]=css(y,xo,n,zp)    % updated May 2014

   [nt,p] = size(xo);

if zp==0; x=xo; Z=[];                                % updated May 2014
elseif zp>0; x=xo(:,1:(p-zp)); Z=xo(:,(p-zp+1):p);    % updated May 2014
end;                                                % updated May 2014

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


%%%%%%% CSS WITHIN %%%%%%%

t1 = 1:t; t1 = t1';  t2 = t1.^2;    
tt = [ones(t,1) t1 t2];

%tt = [ones(t,1)];            % with this line, CSS and Fixed and Random effects estimators are equivalent. 

wd = size(tt,2); 


Pt=tt*inv(tt'*tt)*tt';
Qt=eye(t)-Pt;

%LQ = kron(eye(n),tt);             
%PQ = LQ*inv(LQ'*LQ)*LQ'; 
%MQ = eye(nt) - PQ;           % eye(nt) is a huge matrix. We may have 'out of memory' error. 
%yy = y-PQ*y;
%xx = x-PQ*x; 

       yy = []; 
       xx = [];

    for i = 1:n
       yy = [yy;Qt*y(t*(i-1)+1:t*i,1)];    
       xx = [xx;Qt*x(t*(i-1)+1:t*i,:)];
       end;

betaw = inv(xx'*xx)*xx'*yy;

err=yy-xx*betaw;
ssg=err'*err/(n*(t-wd)-p);                      % updated May 2014
bcov=ssg*inv(xx'*xx);
se_bw=sqrt(diag(bcov));

temp = y-x*betaw;  


% Calculation of R2 
ey=y-mean(y);
Rer=y-x*betaw; 
R2cw=1-(Rer'*Rer)/(ey'*ey);
%R2cw=[R2cw;1-(1-R2cw)*(nt-1)/(nt-p-n)];
R2cw=[R2cw;1-(1-R2cw)*(nt-1)/(n*(t-wd)-p)];      % updated May 2014

thetaw=[];
for i=1:n
thetaw=[thetaw;inv(tt'*tt)*tt'*temp(t*(i-1)+1:t*i,1)];
end;

CSSW=[];
for i=1:n
CSSW=[CSSW tt*thetaw(wd*(i-1)+1:wd*i)];
end;

CSSW=CSSW-mean(mean(CSSW));


%%%%%%% CSS GLS %%%%%%%


W = kron(ones(n,1),tt);         % updated May 2014
LX=[x Z W];                     % updated May 2014
pp = size(LX,2);                % updated May 2014

Q=kron(eye(n),tt);


PW = W*inv(W'*W)*W'; 
temp=y-x*betaw;             
temp = temp-PW*temp; 
%mean(temp)


Delta=0;
for i=1:n
ind=t*(i-1)+1:t*i;
ee=temp(ind);               
Delta=Delta+ (inv(tt'*tt)*tt'*ee*ee'*tt*inv(tt'*tt) -ssg*inv(tt'*tt))/n ;
end;

ch=1;

if ch==1;

Omega=ssg*eye(nt)+Q*kron(eye(n),Delta)*Q';
CGLS=inv(LX'*inv(Omega)*LX)*LX'*inv(Omega)*y;    % updated May 2014

elseif ch==2;                   % not recommended becasue of 'out of memory' error

Pq=Q*inv(Q'*Q)*Q';
Mq=eye(nt)-Pq;
F=Q*(Q'*Q)^(-0.5)*(ssg*eye(n*wd)+(Q'*Q)^(1/2)*kron(eye(n),Delta)*(Q'*Q)^(1/2))^(-0.5)*(Q'*Q)^(-0.5)*Q';
Omega2=Mq/sqrt(ssg)+F;
x2=Omega2*LX; y2=Omega2*y;
CGLS=inv(x2'*x2)*x2'*y2;

end;


betag=CGLS(1:p);

if zp>0;                                 % updated May 2014
   gammag = CGLS((p+1):(p+zp));     
   delta0g= CGLS((p+zp+1):pp);     
   err=y-x*betag -Z*gammag - W*delta0g; 
elseif zp==0;
   delta0g= CGLS((p+zp+1):pp);          
   err=y-x*betag -W*delta0g;            
end;


tglscov = inv(LX'*inv(Omega)*LX);
se_bg = sqrt(diag(tglscov));
se_bg = se_bg(1:p);


% Calculation of R2 

if zp>0;                                 % updated May 2014
   err=y-x*betag -Z*gammag ; 
elseif zp==0;
   err=y-x*betag;            
end;

ey=y-mean(y);
R2cg=1-(err'*err)/(ey'*ey);                   % updated May 2014
%R2cg=[R2cg;1-(1-R2cg)*(nt-1)/(nt-p-n)];
R2cg=[R2cg;1-(1-R2cg)*(nt-1)/(n*(t-wd)-pp)];   % updated May 2014


thetag=[];
for i=1:n
thetag=[thetag;inv(tt'*tt)*tt'*err(t*(i-1)+1:t*i,1)];
end;

CSSG=[];
for i=1:n
CSSG=[CSSG tt*thetag(wd*(i-1)+1:wd*i)];
end;

CSSG=CSSG-mean(mean(CSSG));