mat=readmat('data.txt');
mat=mat(180:1000,633:642);
rmat=flipud(mat);
rmat=diff(log(rmat));
rmat=delzeros(rmat);
h=250;


i=1
for Col = rmat
    %%%%%%%%%%%%%%%%% GARCH FIT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [Coeff,Errors,LLF,Innovations,Sigmas,Summary] = garchfit(Col);
    allsigs(1:length(Sigmas),i) = Sigmas;
    s= garchfit(Col);
    future=garchpred(s, Col, h);
    futmat(1:length(future),i) = future;
        
    i = i+1;
end

%%%%%%%%%%%%%%%%%%%%%%%CORRELATION FORECAST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nmat = zscore(rmat); %N(0,1)
[tt,l] = size(Nmat); 
tt = 1/tt;
[Wmat,Pval]=eig(tt*Nmat'*Nmat); % eigen values, eigen vec
Pmat = Nmat*Wmat; % principal component
i=1;
for ColP=Pmat
    %sigma forecast for PCP
    [CoeffPCP,ErrorsPCP,LLFPCP,InnovationsPCP,SigmasPCP,SummaryPCP] = garchfit(ColP);
    
    %%%%%%%%%%%%%%%%%%%%%GARCH PREDICTOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    allsigsP(1:length(SigmasPCP),i) = SigmasPCP;
    sP= garchfit(ColP);
    futureP= garchpred(sP, ColP, h);
    futureP= [SigmasPCP;futureP];
    futmatP(1:length(futureP),i) = futureP;
    i=i+1;
end
%%%%%%%%%%%%%%%%%%%% TRANSFORM BACK TO ORIGINAL RETURNS%%%%%%%%%%%%%%%%%%%
Smat = std(rmat);
Amat = Wmat.*repmat(Smat,length(Smat),1);
% size of A can be smaller than the number of all PCP
dimA = length(Amat);
%loop on the number of forcasts
for k=1:length(futmatP)
    Diagmat1(:,:,k) =diag(futmatP(k,:)');
    Diagmat(:,:,k) = Diagmat1(:,:,k) .*Diagmat1(:,:,k);
    Covforecast(:,:,k)=Amat*Diagmat(:,:,k)*Amat';
    for i=1:dimA
        for j=1:dimA
            Corforecast(i,j,k)=Covforecast(i,j,k)/(sqrt(Covforecast(i,i,k))*sqrt(Covforecast(j,j,k)));
        end
    end
end



% plots the evolution of the correlation between 2 fixed assets
s1=1
s2=3
x1=[1:length(Col)];
x2=[1:length(Corforecast(1,2,:))];
subplot(3,1,1)
plot(x1,rmat(1:end, s1));
axis([0, length(x2), -Inf, Inf])
subplot(3,1,2)
plot(x1,rmat(1:end, s2));
axis([0, length(x2), -Inf, Inf])
xlabel('time');
ylabel('correlation / returns')
subplot
subplot(3,1,3);
plot(x2,squeeze(Corforecast(s1,s2,:)), '-');
axis([0, length(x2), -Inf, Inf])

% plots a tree
%k=plotcortree(Corforecast,h);
% mesh the time evolution of a correlation matrix



