% Cornwell-Schmidt Estimator (Within ang GLS)
% Written by Wonho Song, December 2002
% E-mail: whsong@kiep.go.kr, whsong73@hotmail.com

function [betaw,se_bw,betag,se_bg,CSSW,CSSG,R2cw,R2cg]=css(y,x,n)

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


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

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

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

       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/(nt-n-p);
bcov=ssg*inv(xx'*xx);
se_bw=sqrt(diag(bcov));

temp = y-x*betaw;  


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


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(3*(i-1)+1:3*i)];
end;

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


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

ssg=err'*err/(n*(t-3));

%x=[x kron(ones(n,1),tt)];
Q=kron(eye(n),tt);

temp=y-x*betaw;


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(x'*inv(Omega)*x)*x'*inv(Omega)*y;

elseif ch==2;   % please ignore this part 

Pq=Q*inv(Q'*Q)*Q';
Mq=eye(nt)-Pq;
F=Q*(Q'*Q)^(-0.5)*(ssg*eye(n*3)+(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*x; y2=Omega2*y;
CGLS=inv(x2'*x2)*x2'*y2;

end;


betag=CGLS(1:p);


tglscov = ssg*inv(x'*x);
se_bg = sqrt(diag(tglscov));
se_bg = se_bg(1:p);

temp = y-x*CGLS;  

ecg=y-x*betag;
ecg=ecg-mean(ecg);

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



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

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

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