//-------------------------------------------------------
// Payoff du put
function res=gain_put(x,K); res=max(K-x,0);endfunction

// Payoff du call 
function res=gain_call(x,K); res=max(x-K,0);endfunction 

// Une fonction pour gérer le décalage de 1 dans les vecteurs
function res=inc(n); res=n+1; endfunction;

function res=prix_eu(x_0,N,gain)
// Calcul du prix europeen a l'instant 0
  // ATTENTION AU DECALAGE DE 1 EN TEMPS ET EN ESPACE 
  U=zeros(inc(N),inc(N));
  for k=[0:N] do
    U(inc(N),inc(k)) = gain(x_0 * u^k * d^(N-k),K)/(1+r)^N;
  end;
  for n=[N-1:-1:0] do // le temps decroit de N-1 a 0
    U(inc(n),inc(0):inc(n)) =  ...
        p*U(inc(n+1),inc(1):inc(n+1))+(1-p)*U(inc(n+1),inc(0):inc(n));
  end;
  res=U(inc(0),inc(0));
endfunction;

function res=prix_eu_n(n,x_0,N,gain)
// Calcul du prix a l'instant n
  res=prix_eu(x_0,N-n,gain);
endfunction
//-------------------------------------------------------
function main_1()
  r_0=0.05;sigma=0.3;
  x_0=100;K=100;

  N=50;
  d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  r=r_0/N;

  // Lorsque p=1/2 et K=x_0 le prix du call = le prix du put (exercice!)
  p=1/2;K=x_0;
  p_put  = prix_eu_n(0,x_0,N,gain_put);
  p_call = prix_eu_n(0,x_0,N,gain_call);
  if abs(p_put - p_call) >= 1000 * %eps then
     printf("WARNING: ces deux prix devrait coincider: %f <> %f\n",...
            p_put,p_call);
  else
     printf("Parfait!\n");
  end

  p= (1+r-d)/(u-d);
  printf("Prix du call : %f\n",prix_eu_n(0,x_0,N,gain_call));
  printf("Prix du put : %f\n",prix_eu_n(0,x_0,N,gain_put));
endfunction
//-------------------------------------------------------
function H=histo_discret(samples, maxsize, plot_flag)
  if exists('%nsp') then
    H=histo_discret_nsp(samples, maxsize, plot_flag);
  else
    H=histo_discret_scilab(samples, maxsize, plot_flag);
  end
endfunction

function H=histo_discret_scilab(samples, maxsize, plot_flag)
// histogramme de tirages selon une loi discrète à valeurs entières
// supposé prendre des valeurs entre 0 et max.
// Si  plot_flag=\%t on trace l'histogramme.
  H=0
  for k=0:maxsize
    // Calcul du nbre de tirages valant k / Taille
    H(k+1)=length(find(samples==k)) ./ size(samples,'*'); 
  end; 

  xbasc;plot2d3(0:maxsize,H);
  f=gcf();
  Dessin=f.children(1).children(1).children(1);
  Dessin.thickness=10;Dessin.foreground=5;
endfunction

// version nsp
function H=histo_discret_nsp(samples, maxsize, plot_flag)
  if nargin <= 2 then  plot_flag=%f;end
  // histogramme de tirages selon une loi discr�te � valeurs enti�res
  // suppos� prendre des valeurs entre 0 et max.
  // Si plot_flag=%f pas de dessin
  H=zeros(1,maxsize+1);
  [y,ind,occ]=unique(samples);
  H(y+1) = occ ./ size(samples,'*');
  if plot_flag then
    xbasc;plot2d3(0:maxsize,H,line_color=2,line_thickness=10);
  end;
endfunction
//-------------------------------------------------------
function res=prix_recursif_am(x,k,N,gain)
  if k==N then res=gain(x); return; end;
  //QUESTION: res = <À COMPLÉTER>
  res = p * prix_recursif_am(x*u,k+1,N) + (1-p) * prix_recursif_am(x*d,k+1,N);
  res = max( res / (1+r), gain(x) );
endfunction

function res=prix_slow_am(x,N,gain)
  res=prix_recursif_am(x,0,N,gain);
endfunction
//-------------------------------------------------------
function res=prix_am(x_0,N,gain)
// Calcul du prix americain a l'instant 0
  U=zeros(N+1,N+1);
  x=x_0 * u.^[0:N] .* d.^[N:-1:0]; // les N+1 points de calcul en N
  U(N+1,1:N+1)=gain(x);// Valeur de U en N

  // ATTENTION AU DECALAGE DE 1 en espace et en temps
  for n=[N-1:-1:0] do // le temps decroit de N-1 a 0
    U(n+1,1:n+1)=p*U(n+2,2:n+2) + (1-p)* U(n+2,1:n+1);
    U(n+1,1:n+1)= U(n+1,1:n+1) / (1+r);// actualisation
    x = x_0 * u.^[0:n] .* d.^[n:-1:0]; // les points de calcul en n
    //QUESTION: U(n+1,1:n+1)= POUR UNE OPTION AMERICAINE ?
    U(n+1,1:n+1)=max(U(n+1,1:n+1),gain(x));// calcul du max
  end;
  res=U(1,1);
endfunction;

function res=prix_am_n(n,x_0,N,gain)
  res=prix_am(x_0,N-n,gain);
endfunction
//-------------------------------------------------------
function main_2()
  sigma=0.3; r_0=0.1;
  K=100;x_0=100;

  N=10;
  r=r_0/N;
  d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  p= (1+r-d)/(u-d);//p=1/2;


  prix_am(x_0,N,gain_put)
  prix_slow_am(x_0,N,gain_put)

  // Les deux algos font ils le même chose ?
  // on verifie : prix_slow(x_0,N) \approx prix(x_0,N)
  p1=prix_slow_am(x_0,N,gain_put);
  p2=prix_am(x_0,N,gain_put);
  if (abs(p1 - p2) >= 100 * %eps) then
     printf("WARNING: ces deux prix devrait coincider: %f <> %f\n",...
            p1,p2);
  else
     printf("Parfait!\n");
  end

  N=1000;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  
  prix_am(x_0,N,gain_put)
endfunction;
//-------------------------------------------------------
function main_3()
  // tracé de courbes
  N=50;
  sigma=0.3;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  K=100;x_0=100;
  r_0=0.1;
  r=r_0/N;
  
  vmin=50;
  vmax=150;

  n=0;
  for x=[vmin:vmax] do
    n=n+1; 
    //QUESTION: courbe_am(n) = <A COMPLÉTER>;
    courbe_am(n)=prix_am(x,N,gain_put);
    //QUESTION: courbe_eu(n) = <A COMPLÉTER>;
    courbe_eu(n)=prix_eu(x,N,gain_put);
  end
  // On compare les courbes "Américaines" et "Européennes"
  plot2d([vmin:vmax],courbe_eu,style=2);
  plot2d([vmin:vmax],courbe_am,style=3);
  plot2d([vmin:vmax],gain_put([vmin:vmax]),style=4);
endfunction;
//-------------------------------------------------------
function main_4()
  sigma=0.3;
  K=100;x_0=100;

  N=50;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  r_0=0.1;r=r_0/N;

  vmin=50;vmax=150;
  // évolution de la courbe en fonction de n
  for n=[0,10,20,30,40,45,47,49,50] do
    i=0;
    for x=[vmin:vmax] do
      i=i+1; 
      //QUESTION: courbe_am(i) = <A COMPLÉTER>;
      courbe_am(i)=prix_am_n(n,x,N,gain_put);
    end
    plot2d([vmin:vmax],courbe_am);
  end
endfunction
//-------------------------------------------------------
function main_5()
  // Avec cet algorithme on peut augmenter N
  // mais il faut renormaliser convenablement u et d.
  // Essayer avec N=10,100,200,...,1000
  
  sigma=0.3;
  K=100;x_0=100;

  r_0=0.1;
  p=1/2;
  N=50;

  // La renormalisation convenable
  // pour converger vers un modèle de Black et Scholes
  r=r_0/N;
  d=1-sigma/sqrt(N);
  u=1+sigma/sqrt(N);

  vmin=50;vmax=150;
  for N=[5,10,20,30,50] do
    r=r_0/N;
    d=1-sigma/sqrt(N);
    u=1+sigma/sqrt(N);
    n=0;
    for x=[vmin:vmax] do
      n=n+1; 
      courbe(n)=prix_am(x,N,gain_put);
    end
    plot2d([vmin:vmax],courbe);
    plot2d([vmin:vmax],max(K-[vmin:vmax],0));// Le payoff / l'obstacle
  end
endfunction
//-------------------------------------------------------
function res=simul_bin(N,p)
// sommes partielles nombre de N tirages a pile ou face (p,1-p)
  bin=grand(1,N,"bin",1,p);// tirages a pile ou face (p,1-p)
  //QUESTION: res=<A COMPLÉTER>;// sommes partielles
  res=cumsum(bin);
endfunction;
//-------------------------------------------------------
function V=prix_am_vect(x_0,N,gain)
// On calcule les valeurs de "v(n,x)" (voir question précédente)
// mais, ici, on les conservent dans un vecteur "V"
// ce qui évitera d'avoir à les re-calculer un grand 
// nombre de fois
  
  // ATTENTION AU DECALAGE DE 1 en ESPACE ET EN TEMPS
  V=zeros(N+1,N+1);
  x=x_0 * u.^[0:N] .* d.^[N:-1:0]; // les points de calcul en N
  V(N+1,1:N+1)=gain(x);

  for n=[N-1:-1:0] do // le temps décroit de N-1 a 0
    V(n+1,1:n+1)=p*V(n+2,2:n+2) + (1-p)* V(n+2,1:n+1);
    V(n+1,1:n+1)= V(n+1,1:n+1) / (1+r);// actualisation
    x=x_0 * u.^[0:n] .* d.^[n:-1:0]; // les points de calcul en  n
    V(n+1,1:n+1)=max(V(n+1,1:n+1),gain(x));
  end;
endfunction;
//-------------------------------------------------------
function res=temps_optimal_option(Bin,V,gain)
  // Calcule le temps d'arrêt optimal associé à la trajectoire 
  // |$(X(1),\ldots,X(N))$| où  |$X(n) =$|  x_0*u^Bin(n)*d^(n-Bin(n))
  // Bin est la somme cumulée de tirages de Bernouilli indépendants
  if gain(x_0) == V(1,1)  then  res=0;return;end;
  for n=1:size(Bin,'*')-1 do
    x=x_0*u.^Bin(n)*d.^(n-Bin(n));// Valeur de X(n)
    //QUESTION: if (ÉCRIRE LA CONDITION D'ARRET & (gain(x) > 0) then 
    if (gain(x) == V(n+1,Bin(n)+1)) & (gain(x) > 0) then 
      res=n;return; 
    end;
  end;
  res=N;
endfunction
//-------------------------------------------------------
function tau=un_test(x_0,N)
  sigma=0.3;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  K=100;
  r_0=0.1;
  r=r_0/N;

  V=prix_am_vect(x_0,N,gain_put);
  Nbre=1000;
  for j=1:Nbre do
    Bin=simul_bin(N,p);
    tau(j)=temps_optimal_option(Bin,V,gain_put);
  end
  histo_discret(tau,max(tau),%t);
endfunction;

function main_6()
  // un test simple "at the money" |$(K=x_0)$|
  N=100;
  x_0=100;tau=un_test(x_0,N);
endfunction;

function main_7()
  N=50;
  x_0=60;tau=un_test(x_0,N);// on exerce toujours en 0
  x_0=70;tau=un_test(x_0,N);// on n'exerce plus (jamais) en 0
  x_0=80;tau=un_test(x_0,N);// On exerce de plus en plus tard ...
  x_0=100;tau=un_test(x_0,N);
  x_0=120;tau=un_test(x_0,N);// Le plus souvent en N
endfunction;
//-------------------------------------------------------
// Teste le cas du call |$(x-K)_+$|
// qui ne s'exerce jamais avant l'écheance |$N$|.
// Le fait que le prix eu du call soit toujours plus grand 
// que l'obstacle permet de le prouver rigoureusement.
  
function main_8()
  function res=gain_call(x); 
  // Payoff du call 
    res=max(x-K,0);
  endfunction

  sigma=0.3; r_0=0.1;
  K=100;x_0=100;
  
  N=50;
  r=r_0/N;
  u=(1+r)*exp(sigma/sqrt(N));d=(1+r)*exp(-sigma/sqrt(N));
  p=1/2;  //p= (1+r-d)/(u-d); // en principe pour Cox-Ross

  V_call=prix_am_vect(x_0,N,gain_call);
  Nbre=1000;
  for j=1:Nbre do
    Bin=simul_bin(N,p);
    tau(j)=temps_optimal_option(Bin,V_call,gain_call);
  end
  histo_discret(tau,max(tau),%t);// s'exerce toujours en N
endfunction;
//-------------------------------------------------------
function proba=compute_proba(x_0,N)
  // Proba que tau<N
  sigma=0.3;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  K=100;
  r_0=0.1;
  r=r_0/N;

  V=prix_am_vect(x_0,N,gain_put);
  Nbre=1000;
  for j=1:Nbre do
    Bin=simul_bin(N,p);
    tau(j)=temps_optimal_option(Bin,V,gain_put);
  end

  //QUESTION: proba=QUE VAUT CETTE PROBA;
  proba=1 - (length(find(tau==N)) ./ size(tau,'*')); 
endfunction;

function main_9()
  N = 50;
  for x_0=[60,70,80,100,120] do
     printf("x_0=%f: %f\n",x_0,compute_proba(x_0,N));
  end
endfunction
//-------------------------------------------------------
function s=frontiere(V)
  // Calcule la frontière d'exercice t -> s(t) t\in[0,N]
  N=size(V);
  N=N(1)-1;
  s=0;
  for n=0:N do
    for j=[0:n] do
      x=x_0*u^j*d^(n-j);
      if V(n+1,j+1) > gain_put(x) then
	s(n+1)=x_0*u^(j-1)*d^(n-j+1);break;
      end;
    end
  end;
endfunction

function main_10()  
  N=1000;
  r_0=0.05;r=r_0/N;
  sigma=0.3;
  K=100;x_0=70;

  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);

  V=prix_am_vect(x_0,N,gain_put);
  front=frontiere(V);
  plot2d(1:size(front,'*'),front);
endfunction;
//-------------------------------------------------------
  function u=compute_u(N)
    u=zeros(N,2);
    u(N,:)=[0,1];
    for n=[N-1:-1:1] do
      temp = (n/(n+1))*u(n+1,1) + (1/(n+1))*u(n+1,2);
      // temp >= 0, donc u(n,0)=max(temp,0)=temp
      u(n,:)=[temp,max(temp,n/N)];
    end;
  endfunction;
//-------------------------------------------------------
function main_11()
  N=1000;
  obstacle=[1:N]/N;

  u=compute_u(N);
  plot2d(1:N,obstacle,style=2);
  plot2d(1:N,u(:,2),style=3);// u(n,1)
  plot2d(1:N,u(:,1),style=4);// u(n,0)
endfunction
//-------------------------------------------------------
function main_12()
  // u(1,2) = P(succés pour la stratégie optimale)
  // On vérifie que cette proba tends vers 1/e = 37\%
  i=0;
  courbe=0;
  valeurs=[10,25,50,100,250,500,1000];
  for N=valeurs do
    u=compute_u(N);
    i=i+1;courbe(i)=u(1,2);
  end
  plot2d(valeurs,courbe)
endfunction
//-------------------------------------------------------
  function res=temps_min(N)
   // Calcule le
   // premier temps où |$u(n,1)=f(n,1)=n/N$|
   // à |$\epsilon$| près.
    epsilon= 10*%eps;
    u=compute_u(N);
    for n=[1:N] do
      if (abs(u(n,2) - (n/N)) < epsilon) then 
   	break;
      end
    end
    res=n;
  endfunction
//-------------------------------------------------------
function main_13()
  // Vérification que temps_min vaut environ |$N/e$|.
  Taille=1000;
  x=0;
  for N=[1:Taille] do
    x(N)=temps_min(N)/N;
  end
  x=x-1/exp(1);
  nb=size(x,'*');
  plot2d(1:nb,x);
endfunction
//-------------------------------------------------------
function [R] = Omega2R(omega)
  // Calcule les rangs d'insertion R pour un |$\omega$| donné
  R=zeros(1,length(omega));
  for n=[1:length(omega)] do
    // classe le vecteur omega(1:n) en croissant
    y=gsort(omega(1:n),'g','i');
    // calcule l'indice de omega(n) dans le tableau classe y
    // c'est à dire son classement parmis les n premiers
    R(n)=find(omega(n)==y);
  end
endfunction

function [omega] = R2Omega(R)
// Calcule |$\omega$| connaissant les rangs d'insertion

  N=length(R);
  temp=zeros(1,N);
  for n=[1:N] do 
     // On place le n-ième individu à sa place: la R(n)-ième
     temp=[temp(1:R(n)-1),n,temp(R(n):n-1)]; 
  end;
  // On obtient une application |{\tt temp}| qui au rang associe l'individu.
  // Mais  |$\omega$|, c'est l'application qui à l'individu asscocie son rang:
  // c'est l'inverse de |{\tt temp}| qui se calcule par
  for  n=[1:N] do omega(temp(n))=n;end;
  // ou si l'on est à l'aise en Scilab: 
  // omega(temp)=1:N;//!!
endfunction

function main_14()
  // test sur une permutation
  N=10;
  omega=grand(1,'prm',[1:N]')';// tirage d'une permutation aléatoire.
  R=Omega2R(omega);
  // Retrouve t'on omega ?
  and(R2Omega(R)==omega)
endfunction
//-------------------------------------------------------
function n=temps_optimal(omega,tau_min)
// Calcule le temps d'arret optimum
// pour la permutation |$\omega$| 

  // Calcule la suite S à partir de omega (omega->R->S)
  S=(Omega2R(omega)==1);

  // Le temps optimal se situe après |$\tau_min$| et c'est le premier instant
  // où |$S(n)=1$| apres ce temps, sauf si |$n=N$|, auquel cas on est oblige
  // de prendre  le dernier candidat.
  for n=[tau_min:N] do
    //QUESTION: if CONDITION D'ARRET OPTIMAL then break; end;
     if S(n) then break; end;
  end
  // Si on sort de cette boucle avec n=N et N est bien optimal.
endfunction
//-------------------------------------------------------
function main_15()
  // Verification que la probabilité d'obtenir
  // le meilleur est de l'ordre de 1/e
  N=100;
  Taille=10000;
  tau_min=temps_min(N);
  
  ss=0;
  for i=[1:Taille] do
    omega=grand(1,'prm',[1:N]');
    tirages(i)=temps_optimal(omega,tau_min);
    
    // tirages(i) est il le meilleur ?
    if (omega(tirages(i))==1) then  ss=ss+1;end;
  end
  proba=ss/Taille;
  printf("probabilité d''obtenir le meilleur: %.3f\n",proba);

  // histo_discret(tirages,N,%t);
endfunction
//-------------------------------------------------------