% Computing the singular factors of models A, B, C and using it for
%prediction at various horizons.

%----------------------
%Parameters
%----------------------
N=2 %number of factors
model = 1 %1 is for model C, 2 is for model B, 3 is for model A

%--------------------------------------------------------------------------------
% Load and transform the data
%--------------------------------------------------------------------------------
load ED_rates_CRB.txt
data=100-ED_rates_CRB;
dataini=data;
%data=data(1:5:2486,:); %selecting weekly data
%data=data(1:21:2486,:); %selecting monthly data
%data=data(1:136:2486,:); %selecting 6-month data
%data=data(1:252:2486,:); %selecting yearly data
[m0,n]=size(data); %m0 is number of dates, %n is number of points in the curve
mdata=mean(data);
data=data-kron(ones(m0,1),mdata);%demean the data
data=data';
change=(data(:,2:m0)-data(:,1:m0-1));

%--------------------------------------------
%Restricting sample for out of sample analysis
%----------------------------------------------

%m=round(m0/2) %using half of sample
m=m0-1

%----------------------------------------------------------------
% Estimate covariance and cross-covariance operators
%----------------------------------------------------------------




if model == 1 
   autocov=data(:,1:m)*data(:,1:m)'/m;%autocovariance
   crosscov=change(:,1:m-1)*data(:,1:m-1)'/(m-1);%crosscovariance G21 (model C: the future change on the current curve)
else if model == 2
      autocov=change(:,1:m-1)*change(:,1:m-1)'/(m-1);%autocovariance
      crosscov=change(:,2:m-1)*change(:,1:m-2)'/(m-2); %(model B: the future change on the past change)
   else if model == 3
         autocov=data(:,1:m)*data(:,1:m)'/m;%variance of the predictor
         crosscov=data(:,2:m)*data(:,1:m-1)'/(m-1); %(model A: the future value on the past value)
      end
   end
end


Vp=crosscov'*crosscov;
Vf=crosscov*crosscov';
[R,D]=eig(Vp);  %Ep is singular factors, Dp is singular values
D=diag(D);
D=flipdim(D,1);
R=flipdim(R,2);
R=R(:,1:N);%first N eigenvectors
D=D(1:N,1)%first N eigenvalues
F=(crosscov*R)*diag(D.^(-1/2));	%Ef is singular factor loadings, normalized to unit norm.
G=autocov*R;

subplot(1,3,1)
plot([1:n]',R(:,1),'b-',[1:n]',R(:,2),'m:')
subplot(1,3,2)
plot([1:n]',F(:,1),'b-',[1:n]',F(:,2),'m:')
subplot(1,3,3)
plot([1:n]',G(:,1),'b-',[1:n]',G(:,2),'m:')




if model == 1 
   FactorV=G'*data; %singular factor variates (models C)
else if model == 2
      FactorV=G'*change; %(model B)
   else if model == 3
         FactorV=G'*data; %(model A)
      end
   end
end

%---------------
%predicting
%-------------
predict=F*diag(D.^(1/2))*FactorV;
%--------------
%evaluating predictions
%(this criterion aims to compare performance of a given predition model to that of random walk model)
%--------------
if model == 1 
   criterion=sum(sum(change(:,m:m0-1).*predict(:,m:m0-1))')/(m-1)*n %evaluation (model C)
else if model == 2
      criterion=sum(sum(change(:,m+1:m0-1).*predict(:,m:m0-2))')/(m-1)*n %evaluation (model B)
   else if model == 3
       criterion=sum(sum(data(:,m+1:m0).*(predict(:,m:m0-1)-data(:,m:m0-1)))')/(m-1)*n  %evaluation (model A)
      end
   end
end