// TP6 Calcul de la DIC, DOC, DIP, DOP
clear;
  
n=12;T=1;sigma=0.2;deltat=T/n;S0=150;
K=S0;
//
L=130;
//L=K-1;
up=exp(sigma*(sqrt(deltat)));down=(1/up);
//SS=zeros(n+1,n+1);
//SS(1,1)=S0;
//for i=1:n
//  SS(i+1,1)=SS(i,1)*down; // ici j=0
//  for j=1:i
//    SS(i+1,j+1)=SS(i,j)*up;
//  end;
//end;
//
function s=S(i,j)
    s=S0*up^j.*down^(i-j);
endfunction;


function phi=phi(S);
phi=max(S-K,0);
endfunction;

r=0.05;R=exp(r*deltat);
p=((exp(r*deltat)-down)/(up-down));

CC=zeros(n+1,n+1);
for j=0:n
CC(1+n,1+j)=phi(S(n,j));
end;
for i=n-1 :-1 :0
for j=0 :i
CC(1+i,1+j)=(p*CC(1+i+1,1+j+1)+(1-p)*CC(1+i+1,1+j))/R;
end;
end;

CCC=zeros(n+1,n+1);
CCC(n+1,:)=phi(S(n,0:n));
for i=0:n-1
  for j=0:i
    CCC(i+1,j+1)=(phi(S(n,j:j+n-i))*binomial(p,n-i)')/R^(n-i);
  end
end



/// Indicatrice de la région sous la barriere/////

function indi=sousL(i,j);
    if S(i,j)<=L then indi=1; else indi=0;end;
endfunction;

////Calcul des valeurs de la DIC////

DICC=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==1 then DICC(n+1,j+1,k+1)=phi(S(n,j));
        else DICC(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          DICC(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*DICC(1+i+1, 1+j+1,1+max(k, sousL(i+1,j+1)))..
           + (1-p)*DICC(1+i+1,1+j,1+max(k,sousL(i+1,j)))   )/R;
        end;
    end;
end;
DIC=DICC(1,1,1)
DICC(1,1,2)
disp("DIC="); disp(DIC); //On suppose que S0>L
////Calcul des valeurs de la DOC////
D0CC=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==0 then DOCC(n+1,j+1,k+1)=phi(S(n,j));
        else DOCC(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          DOCC(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*DOCC(1+i+1, 1+j+1,1+max(k, sousL(i+1,j+1)))..
           + (1-p)*DOCC(1+i+1,1+j,1+max(k,sousL(i+1,j)))   )/R;
        end;
    end;
end;
DOC=DOCC(1,1,1); //On suppose que S0>L
disp("DOC="); disp(DOC);
disp("Call="); disp(CC(1,1));
DIC
DOC
DIC+DOC
DIC+DOC-DICC(1,1,2)
CC(1,1)
DICC(1,1,2)
////////////////////// DIC+DOC=C //////////////////
max(abs(DICC(:,:,1)+DOCC(:,:,1)-CC))
//////////////////// PP /////////////
function psi=psi(S);
psi=max(K-S,0);
endfunction;

PP=zeros(n+1,n+1);
for j=0:n
PP(1+n,1+j)=psi(S(n,j));
end;
for i=n-1 :-1 :0
for j=0 :i
PP(1+i,1+j)=(p*PP(1+i+1,1+j+1)+(1-p)*PP(1+i+1,1+j))/R;
end;
end;

 ////////////////DIP,DOP

////Calcul des valeurs de la DIP////

DIPP=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==1 then DIPP(n+1,j+1,k+1)=psi(S(n,j));
        else DIPP(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          DIPP(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*DIPP(i+1+1, j+1+1, max(k, sousL(i+1,j+1))+1)..
           + (1-p)*DIPP(i+1+1,j+1,max(k,sousL(i+1,j))+1)     )/R;
        end;
    end;
end;
DIP=DIPP(1,1,1);
disp("DIP="); disp(DIP); //On suppose que S0>L
////Calcul des valeurs de la DOC////

DOPP=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==0 then DOPP(n+1,j+1,k+1)=psi(S(n,j));
        else DOPP(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          DOPP(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*DOPP(i+1+1, j+1+1, max(k, sousL(i+1,j+1))+1)..
           + (1-p)*DOPP(i+1+1,j+1,max(k,sousL(i+1,j))+1)     )/R;
        end;
    end;
end;
DOP=DOPP(1,1,1); //On suppose que S0>L
disp("DOP="); disp(DOP);
//
disp("Put="); disp(CC(1,1)-S0+exp(-r*T)*K);//relation Call-Put
disp("DIP+DOP=");disp(DIP+DOP);
max(abs(DIPP(:,:,1)+DOPP(:,:,1)-PP))




////////////////// UIC,UOC,UIP,UOP //////////////////
////Calcul des valeurs de la DIC////
L=170
function indi=susL(i,j);
    if S(i,j)>=L then indi=1; else indi=0;end;
endfunction;
UICC=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==1 then UICC(n+1,j+1,k+1)=phi(S(n,j));
        else UICC(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          UICC(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*UICC(1+i+1, 1+j+1,1+max(k, susL(i+1,j+1)))..
           + (1-p)*UICC(1+i+1,1+j,1+max(k,susL(i+1,j)))   )/R;
        end;
    end;
end;
UIC=UICC(1,1,1)
UICC(1,1,2)
////Calcul des valeurs de la UOC////
U0CC=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==0 then UOCC(n+1,j+1,k+1)=phi(S(n,j));
        else UOCC(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          UOCC(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*UOCC(1+i+1, 1+j+1,1+max(k, susL(i+1,j+1)))..
           + (1-p)*UOCC(1+i+1,1+j,1+max(k,susL(i+1,j)))   )/R;
        end;
    end;
end;
UOC=UOCC(1,1,1) //On suppose que S0<L

////////////////////// UIC+UOC=C //////////////////
max(abs(UICC(:,:,1)+UOCC(:,:,1)-CC))
//////////////////// PP /////////////
function psi=psi(S);
psi=max(K-S,0);
endfunction;

PP=zeros(n+1,n+1);
for j=0:n
PP(1+n,1+j)=psi(S(n,j));
end;
for i=n-1 :-1 :0
for j=0 :i
PP(1+i,1+j)=(p*PP(1+i+1,1+j+1)+(1-p)*PP(1+i+1,1+j))/R;
end;
end;

 ////////////////DIP,DOP

////Calcul des valeurs de la UIP////

UIPP=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==1 then UIPP(n+1,j+1,k+1)=psi(S(n,j));
        else UIPP(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          UIPP(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*UIPP(i+1+1, j+1+1, max(k, susL(i+1,j+1))+1)..
           +(1-p)*UIPP(i+1+1,j+1,max(k,susL(i+1,j))+1)     )/R;
        end;
    end;
end;
UIP=UIPP(1,1,1)
////Calcul des valeurs de la UOC////

UOPP=zeros(n+1,n+1,2);
for j=0 :n //On commence avec la valeur à t=T
    for k=0 :1 //k=1 code "on est passé en dessous de L"
        if k==0 then UOPP(n+1,j+1,k+1)=psi(S(n,j));
        else UOPP(n+1,j+1,k+1)=0;
        end;
    end;
end;
for i=n-1 :-1 :0 //récurrence "backward"
    for j=0 :i
        for k=0 :1 //k'=1 code "on est passé en dessous de L"
          UOPP(i+1,j+1,k+1)=.. // ne s'enclenche que pour t>0
             (  p*UOPP(i+1+1, j+1+1, max(k, susL(i+1,j+1))+1)..
           + (1-p)*UOPP(i+1+1,j+1,max(k,susL(i+1,j))+1)     )/R;
        end;
    end;
end;
UOP=UOPP(1,1,1); //On suppose que S0<L

max(abs(UIPP(:,:,1)+UOPP(:,:,1)-PP))
