clear;
c=10;
lambda=10;
tau=[10^3,10^3,10^3];

tet=zeros(8,1);
tet(1,1:2)=[10^-3,10^-3];
tet(2,1:2)=[10^-2,10^-2];
tet(3,1:2)=[10^-1,10^-1];
tet(4,1:2)=[10^-0,10^-0];
tet(5,1:2)=[10^1,10^1];
tet(6,1:2)=[10^2,10^2];
tet(7,1:2)=[10^3,10^3];
tet(8,1:2)=[10^4,10^4];

for i=1:8

    teta1=tet(i,1:2);

for j=1:tau(3) %quantidade de processos a serem simulados
    
    for ud=1:2 %processo de cima de de baixo
            for bd=1:2 %processo base ou derivado
                gama(j,ud,bd,1)=0; k=1;
                if bd==1 % processo base
                    while gama(j,ud,bd,k)<tau(1) %precisão da simulação
                        k=k+1; gama(j,ud,bd,k)=gama(j,ud,bd,k-1)+exprnd(1);
                    end
                    n(j,ud)=k-1;
                    F=rand([n(j,ud),1]);
                    U=rand([n(j,ud),2]);
                    V=exprnd(1,[n(j,ud),2]);
                elseif bd==2 % cópula de cleyton
                    for z=2:n(j,ud) %primeiro termo igual a zero
                        gama(j,ud,bd,z)=gama(j,ud,1,z)*(F(z)^(-teta1(ud)/(1+teta1(ud)))-1)^(-1/teta1(ud));
                    end
                else % cópula de frank
                    for z=2:n(j,ud) %primeiro termo igual a zero
                        gama(j,ud,bd,z)=1/teta2(ud)*log(1+exp(teta2(ud)*gama(j,ud,1,z)))*(F(z)/(1-F(z)));
                    end
                end

                for t=1:tau(2) %precisão da discretização do processo
                    if bd==1
                            Z(1)=0;
                        for z=2:n(j,ud)
                            Z(z)=Z(z-1)+(U(z,1)<=(t/tau(2)))*exp(-gama(j,ud,bd,z)/c)/lambda*V(z,1);
                        end
                        X(j,ud,bd,t)=Z(n(j,ud));
                    else %gerador aleatório U e V são iguais para ambas as cópulas
                            Z(1)=0;
                        for z=2:n(j,ud)
                            Z(z)=Z(z-1)+(U(z,2)<=(t/tau(2)))*exp(-gama(j,ud,bd,z)/c)/lambda*V(z,2);
                        end
                        X(j,ud,bd,t)=Z(n(j,ud));
                    end
                end

            end
    end

    for bd=1:2
        for t=1:tau(2)
            x(t,j,bd)= X(j,1,bd,t)-X(j,2,bd,t);
        end
    end

    correl(i,j)=corr(x(:,j,1),x(:,j,2));
        
end

r(i) = sum(correl(i,:))/tau(3);

subplot(4,2,i)
plot(x(tau(2),:,1),x(tau(2),:,2),'.')
title(strcat('\theta_{u}=', num2str(teta1(1)),',','\theta_{d}=', num2str(teta1(2)),',','\rho=',num2str(r(i))))
xlabel('Processo original')
ylabel('Processo dependente')

end
