clear;
c=10;
lambda=10;
tau=[10^3,10^3,10^3];
teta1=[0.5,0.5];
teta2=[1,1];

R = zeros(tau(3),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:3 %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:3
        for t=1:tau(2)
            x(t,j,bd) = X(j,1,bd,t)-X(j,2,bd,t);
%            x(j,bd,t)= X(j,1,bd,t)-X(j,2,bd,t);
        end

        for t=1:tau(2)
            M(j,bd)=0;
            G(j,bd)=1;
            for t=1:tau(2)
                M(j,bd)=x(t,j,bd)+M(j,bd);
                G(j,bd)=(1+x(t,j,bd))*G(j,bd);
            end
        end
        M(j,bd)=M(j,bd)/tau(2);
        G(j,bd)=real(G(j,bd)^(1/tau(2))-1);
            
    end
    
    R(j,1)=corr(x(:,j,1),x(:,j,2));
    R(j,2)=corr(x(:,j,1),x(:,j,3));
    
end

RS = sum(R)/tau(3);

r=0.08;
S(1)=40;
S(2)=60;
alpha=[0:0.2:1];
K=[20:10:80];
base=zeros(length(K)*length(alpha),8);

for l=1:length(alpha)

Pm0=[S(1)*exp(M(:,1)),S(2)*exp(M(:,2)),S(2)*exp(M(:,3))];
Pg0=[S(1)*exp(G(:,1)),S(2)*exp(G(:,2)),S(2)*exp(G(:,3))];
Pm1=[Pm0(:,1)*alpha(l)+Pm0(:,2)*(1-alpha(l)),Pm0(:,1)*alpha(l)+Pm0(:,3)*(1-alpha(l))];
Pg1=[Pg0(:,1)*alpha(l)+Pg0(:,2)*(1-alpha(l)),Pg0(:,1)*alpha(l)+Pg0(:,3)*(1-alpha(l))];
  
    for s=1:length(K)
        
        Cm=[max(Pm1(:,1)-K(s),0),max(Pm1(:,2)-K(s),0)]*exp(-0.08/252*tau(2));
        Cg=[max(Pg1(:,1)-K(s),0),max(Pg1(:,2)-K(s),0)]*exp(-0.08/252*tau(2));
        Pm=[max(K(s)-Pm1(:,1),0),max(K(s)-Pm1(:,2),0)]*exp(-0.08/252*tau(2));
        Pg=[max(K(s)-Pg1(:,1),0),max(K(s)-Pg1(:,2),0)]*exp(-0.08/252*tau(2));
    
        AmC=[mean(Cm(:,1)),mean(Cm(:,2))];
        AgC=[mean(Cg(:,1)),mean(Cg(:,2))];
        AmP=[mean(Pm(:,1)),mean(Pm(:,2))];
        AgP=[mean(Pg(:,1)),mean(Pg(:,2))];

        base(s+length(K)*(l-1),1:8)=[AmC,AgC,AmP,AgP];   
    end
end