//- Premiere partie du TP
//- Illustration du prix de l'actif, du payoff, de la valeur du
//portefeuille et de sa composition

clear;
xdel([1 2 3 4]);
xdel(0);
//------------------------------------
//- Parametres numeriques et financiers
//-------------------------------------
//- Prix a l'instant zero
S0 = 40;
//- Strike
K = 45;
//- Pas de discretisation
N = 10;
//- Maturite
T = 1;
//- Pas de temps
dt = T/N;
//- Taux d interet instantané sans risque
R = 0.03;
//- Rendement sur une periode de durée dt
r = exp(R*dt)-1;     
//- Volatilite de notre modele
vol_mod = 0.2;
//- Hausse et baisse
a_mod = (1+r) * exp(-vol_mod*sqrt(dt))-1;     
b_mod = (1+r) * exp(vol_mod*sqrt(dt))-1;
p = (b_mod-r)/(b_mod-a_mod);

S(1)=S0;

//- Initialisation de l'arbre
Prix=zeros(N+1,N+3);
Prix(N+1,:) = max(S0*(1+a_mod)^N*(((1+b_mod)/(1+a_mod))^(-1:(N+1)))-K,0);
//- Prix(n,i) contient le prix au temps (n-1)dt partant de S0*(1+a)^(n-i+1)*(1+b)^(i-2)
p1 = p/(1+r);
p2 = (1-p)/(1+r);

for j=N:-1:1
Prix(j,1:(j+2))=Prix(j+1,1:(j+2)) * p1+Prix(j+1,2:(j+3)) * p2;
end
//- Fin initialisation arbre

//- Graphe de l'arbre si N<=20
if(N<=20) then
  for j=N+1:-1:2
    xset("window",4)
    plot2d((ones(j-1,1)*[(j-2)*dt,(j-1)*dt])',([Prix(j-1,2:j);Prix(j,2:j)]),style=2*ones(j-1,2))
    plot2d((ones(j-1,1)*[(j-2)*dt,(j-1)*dt])',([Prix(j-1,2:j);Prix(j,3:j+1)]),style=2*ones(j-1,2))
    for i=2:j+1
      xstring((j-1)*dt,Prix(j,i),string(round(Prix(j,i)*10)/10));
    end
  end
xstring(0,Prix(1,2),string(round(Prix(1,2)*10)/10));
end
//- Fin graphe de l'arbre

V(1)=Prix(1,2);
H_1(1)=(Prix(2,3)-Prix(2,2))/(S0*(b_mod-a_mod));
H_0(1)=V(1)-H_1(1)*S(1);
composition(1)=H_1(1)*S(1)/V(1);

nba=0; //- compte le nb de realisation "1+a"

for n=2:N
  Bern=1*(rand(1,1)<p);
  nba=nba+Bern;
  S(n) = S(n-1)*(1+b_mod+(a_mod-b_mod)*Bern);    
  V(n)= H_0(n-1)*(1+r)^(n-1)+H_1(n-1)*S(n);
  H_1(n)= (Prix(n+1,n-nba+2)-Prix(n+1,n-nba+1))/(S(n)*(b_mod-a_mod));
  H_0(n)=(V(n)-H_1(n)*S(n))/(1+r)^(n-1);
  payoff(n)=max(S(n)-K,0);
  composition(n)=H_1(n)*S(n)/V(n);
end
S(N+1)=S(N)*(1+b_mod+(a_mod-b_mod)*(rand(1,1)<p));    
V(N+1)=H_0(N)*(1+r)^(N)+H_1(N)*S(N+1);
payoff(N+1)=max(S(N+1)-K,0);
composition(N+1)=H_1(N)*S(N+1)/V(N+1);

xset("window",1)
xtitle("Actif","Temps","Prix");
plot2d((0:dt:T)',S,style=2)
legends(['Prix CRR'],[2],'ur')
xset("window",2)
plot2d((0:dt:T)',payoff,style=2)
plot2d((0:dt:T)',V,style=3)
legends(['Payoff Call','Valeur Portefeuille'],[2,3],'ur')
xtitle("Payoff du Call et Valeur du portefeuille","Temps");
xset("window",3)
xtitle("Composition du portefeuille","Temps","H1*S/V");
plot2d((0:dt:T)',composition,style=2)
legends(['Composition portefeuille'],[2],'ur')

//- Erreur de couverture a maturite
valeur_option=Prix(1,2)
erreur_couverture = payoff(N+1)-V(N+1)