% Cornwell-Schmidt Estimator (Within only)
% 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,CSSW,R2cw]=cssw(y,x,n)

   [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];

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;                   % This is a huge matrix. We may have 'out of memory' error. 
%yy = MQ*y;
%xx = MQ*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-3)-p);                % updated on 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-3)-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(3*(i-1)+1:3*i)];
end;

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