clear all; clc;
tic;
load data.txt;
windowlength=500;
% Below is the GARCH-copula code
%for i=1:size(data,1)-windowlength
%for i=1:1000
%    htx(i,:)=conditionalvar(data(i:i+windowlength-1,1), 1, 1,'NORMAL');
%    hty(i,:)=conditionalvar(data(i:i+windowlength-1,2), 1, 1, 'NORMAL');
%    forecast1(i,:)=GARCHforecast(data(i:i+windowlength-1,1),'NORMAL');
%    forecast2(i,:)=GARCHforecast(data(i:i+windowlength-1,2),'NORMAL');
%end

%Below is the GJR Copula code
for i=1:size(data,1)-windowlength
    htx(i,:)=gjrvar(data(i:i+windowlength-1,1), 1, 1, 1, 'NORMAL');
    hty(i,:)=gjrvar(data(i:i+windowlength-1,2), 1, 1, 1, 'NORMAL');
    forecast1(i,:)=GJRforecast(data(i:i+windowlength-1,1),'NORMAL');
    forecast2(i,:)=GJRforecast(data(i:i+windowlength-1,2),'NORMAL');
end
htx=htx';
hty=hty';
for i=1:size(data,1)-windowlength
    stderrorx(:,i)=data(i:i+windowlength-1,1)./sqrt(htx(:,i));
    stderrory(:,i)=data(i:i+windowlength-1,2)./sqrt(hty(:,i));
    output(i,:)=corr(stderrorx(:,i),stderrory(:,i));
    %[rho(i,:),kappa(:,i)]=copula_tvp(stderrorx(:,i),stderrory(:,i));
    %output(i,:) = P_CopulaForecast(rho(i,:), kappa(:,i),[stderrorx(:,i) stderrory(:,i)],2,13);
end

kappa1=output;
[T,k] = size(stderrorx);
for i=1:k
u(:,i) = empiricalCDF(stderrorx(:,i));
v(:,i) = empiricalCDF(stderrory(:,i));
end
options = optimset('Display','off','TolCon',10^-12,'TolFun',10^-4,'TolX',10^-6);
lower = 0.0001;
theta0 = 1;

for i=1:k
    lower1 = [-0.99 , 2.1 ];
    upper1 = [ 0.999, 100 ];
    theta01(:,i) = [kappa1(i);10];
end
for i=1:k
    kappa8(i,:) = fmincon('tcopulaCL',theta01(:,i),[],[],[],[],lower1,upper1,[],options,[u(:,i),v(:,i)]);
end

kappa8=kappa8(:,1);
%Volatility forecast
htab=kappa8.*sqrt(forecast1.*forecast2);
ratio=htab./forecast2;
toc;