//-------------------------------------------------------
function X=tirage_pf(n,p)
// Effectue n tirage a Pile (P) ou face (F) (p,1-p)
  X=grand(1,n,"bin",1,p);
  X=string(X);// transforme le tableau de nombres en 
              // tableau de caractères
  X=strsubst(X,'0','F');// remplace les 0 par des F
  X=strsubst(X,'1','P');// remplace les 1 par des P  
endfunction

function [MAX] = max_length(U)
// Calcule le nombre maximum de 1 consecutifs
// dans la suite U
  MAX=0;N=0;
  for n=[1:size(U,'*')] do
    //QUESTION: if U(n)=='P' then <À COMPLÉTER> else <À COMPLÉTER> end;
    if U(n)=='P' then N=N+1; else N=0; end;
    MAX=max(MAX,N);
  end
endfunction
//-------------------------------------------------------
function H=histo_discret(samples, maxsize, varargin)
// histogramme de tirages selon une loi discrète à valeurs entières
// supposé prendre des valeurs entre 0 et max.
// Si tracer_histo=%f pas de dessin
  H=0
  for k=0:maxsize
    // Calcul du nbre de tirages valant k / Taille
    H(k+1)=length(find(samples==k)) ./ size(samples,'*'); 
  end; 
  
  [lhs,rhs]=argn(0);
  if (rhs == 3) & (varargin(1)==[%t]) then
    xbasc;plot2d3(0:maxsize,H);
  end;
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 tracer_histo=%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 main_1()
  p=1/2;
  // On teste cette fonction avec N=20
  // rand(1:20) renvoie 100 tirages uniformes sur l'intervalle [0,1]
  U=tirage_pf(20,p)// 20 tirages a pile ou face 1/2,1/2  (1=Pile,0=Face)
  max_length(U)// nombre mximum de P consecutifs

  // On effectue 1000 tirages avec N=100
  N=100;X=0;
  Taille=1000;// nbre de simulation
  for i=[1:Taille] do
    U=tirage_pf(N,p);
    X(i)=max_length(U);
  end;
  
  // On trace l'histogramme si plot_flag==%t
  plot_flag=%t
  if exists('%nsp') then
    histo_discret_nsp(X,20,plot_flag);
  else
    histo_discret(X,20,plot_flag);
  end
endfunction
//-------------------------------------------------------
function [X] = trajectoire(U)
// On calcule la trajectoire de X
  X=0;val=0;
  for n=[1:prod(size(U))] do
    //QUESTION: if U(n)=='P' then <À COMPLÉTER> else <À COMPLÉTER> end;
    if U(n)=='P' then val=val+1; else val=0; end;
    X(n)=val;
  end
endfunction
//-------------------------------------------------------
function main_2()
  N=100;
  p=1/2;
  rectangle=[1,0,N,8];

  // On trace une trajectoire de X
  U=tirage_pf(N,p);
  X=trajectoire(U);
  xbasc;plot2d2(1:N,X,rect=rectangle);
  xclick();// attend un click dans la fenêtre graphique
  
  // 4 autres trajectoires
  for i=1:4 do
    U=tirage_pf(N);
    plot2d2(1:N,trajectoire(U),rect=rectangle,style=i+1);
    xclick();// attend un click dans la fenêtre graphique
  end;
endfunction;
//-------------------------------------------------------
function Pro= proba(N,l,p)
  // Calcule la probabilite de voir au moins l piles consecutifs
  // dans N tirages a pile (p) ou face (1-p)

  // la matrice de transition de la chaîne 
  // de taille (l+1,l+1)
  P=[p*diag(ones(l,1));zeros(1,l-1),1];
  P=[[(1-p)*ones(l,1);0],P];
  PN=P^N;
  //QUESTION: Pro = QUE VAUT LA PROBA DE VOIR AU MOINS l PILES;
  Pro = PN(1,l+1);
endfunction
//-------------------------------------------------------
function [loi]=calculer_loi(N,p)
  MAXMAX=30;// a choisir assez grand (mais pas trop ...) 
  // attention au decalage de 1 pour les vecteurs : 
  //      loi(1)=P(X=0), ..., loi(n+1)=P(X=n)
  
  previous=1;// proba d'avoir au moins 0 pile = 1 !
  // le support de la loi est [0,1,...,N] que l'on tronque en MAXMAX
  for l=0:min(N,MAXMAX) do
    //QUESTION: current=<À COMPLÉTER>;// proba d'avoir au moins l+1 pile
    current=proba(N,l+1,p);// proba d'avoir au moins l+1 pile
    loi(l+1)=previous - current;// proba d'avoir exactement l pile
    previous=current;
  end
  loi=loi(:)';// c'est plus joli comme ca !
endfunction
//-------------------------------------------------------
function main_3()
  // On teste avec N=1 et N=2, p=1/2
  // Pour N=1, 0 pile avec proba 1/2 et 1 pile avec proba 1/2
  calculer_loi(1,1/2)
  // Pour N=2, on doit trouver (1/4,1/2,1/4) pour (0,1,2)
  calculer_loi(2,1/2)
  // en principe ca marche ...
  
  N=100;p=1/2;
  loi=calculer_loi(N,p);
  sum(loi); // on verifie que ca somme a 1
  
  // dessin
  // ATTENTION on est decale de 1, donc 0:20 devient 1:21
  
  if exists('%nsp') then
    xbasc;plot2d3(0:20,loi(1:21),line_thickness=10,line_color=2);
  else
    xbasc;plot2d3(0:20,loi(1:21));
    f=gcf();Dessin=f.children(1).children(1).children(1);Dessin.thickness=10;
  end
  
  // comparaison avec les simulations
  Taille=10000;
  for i=[1:Taille] do
    U=tirage_pf(N);
    X(i)=max_length(U);
  end
  
  histo=0;
  for i=0:20 do
    histo(1,i+1)=sum(X==i)/Taille;
  end
  
  epsilon=norm(histo(1:21)-loi(1:21)) 
          // doit etre "petit", pour bien faire il faudrait faire un
          // test du |$\xi^2$| pour savoir ce que "petit" veut dire ici.
  printf("epsilon=%f, petit en principe.\n",epsilon)
endfunction
//-------------------------------------------------------
function main_4()
// Test du critère lorsque les tirages sont aléatoires
// Ca marche "souvent" mais pas "toujours".
  N=100;
  p=1/2;
  Taille=10;
  for i=[1:Taille] do  
    U=tirage_pf(N,p);
    if max_length(U) >= 4 then 
      printf("OK\n");
    else
      // Ca arrive "rarement" mais ca arrive 3 fois sur 100
      printf("test négatif\n");
    end
  end;
endfunction;
//-------------------------------------------------------
function main_5()
// Calcul du maximum de vraisemblance de la loi 
// imax = indice du maximum,  m = le maximum

  p=1/2;
  printf('N: indice du max -> maximum\n');
  for N=[10,100,1000] do
    loi=calculer_loi(N,p);
    //QUESTION: [m,imax]= ON CHERCHE LE MAX DE LA LOI ET L'INDICE DU MAX
    [m,imax]=max(loi);
    imax=imax-1;// A cause du déclage de 1
    printf('%d: %d -> %f\n',N,imax,m);
  end
endfunction
//-------------------------------------------------------
function s=moyenne(loi)
  N=size(loi,'*')-1;
  //QUESTION: s = CALCULER LA MOYENNE D'UNE V.A. DE LOI 'loi' ?
  s=[0:N]*loi';
endfunction
//-------------------------------------------------------
function main_5_bis()
// La moyenne varie (approximativement) comme C log(N).
// Ca peut se prouver.
  p=1/2;
  i=0;
  for N=[10,50,100,500,1000,5000,10000] do
    i=i+1;
    loi=calculer_loi(N,p);
    x(i)=log(N);
    y(i)=moyenne(loi);
    printf('%d: %f ?=? %f\n',N,y(i),x(i));
  end;
  plot2d(x,y);
endfunction
//-------------------------------------------------------
function X=simul_cox_ross(N,x_0,p,u,d)
  U=grand(1,N,"bin",1,p);// tirages a pile ou face (p,1-p)
  X=x_0 * [1, cumprod(u.^U .* d.^(1-U))];
  // Plus long mais plus comprehensible ...
  // X(1)=x_0;
  // for i=1:prod(length(U)) do
  //  if U(i) then  
  //    X(i+1)=X(i)*u
  //  else
  //    X(i+1)=X(i)*d;
  //  end    
  //end;
endfunction;
//-------------------------------------------------------
function main_6()
  N=50;
  sigma=0.3;
  p=1/2;u=1-sigma/sqrt(N);d=1+sigma/sqrt(N);
  x_0=1;

  X=simul_cox_ross(N,x_0,p,u,d);
  plot2d2(0:N,X);
endfunction
//-------------------------------------------------------
K=100;
function res=f(x)
  res=max(x-K,0);
endfunction
//-------------------------------------------------------
function res=prix_recursif(x,k,N)
  if k==N then res=f(x); return; end;
  //QUESTION: res = RECOPIER L'EQUATION DU COURS !;
  res = p * prix_recursif(x*u,k+1,N) + (1-p) * prix_recursif(x*d,k+1,N);
endfunction

function res=prix_slow(x,N)
  res=prix_recursif(x,0,N);
endfunction
//-------------------------------------------------------
N=10;

// On choisit des paramètres pour converger 
// vers le modèle de Black et Scholes.
sigma=0.3;
p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);

prix_slow(100,N);// Faire N=20 pour savoir ce que slow veut dire !
//-------------------------------------------------------
function res=inc(n);
// Permet de faire comme si les tableaux était indicés 
// à partir de 0 (et non de 1)
 res=n+1;
endfunction;
//-------------------------------------------------------
function res=prix(x_0,N)
  // U=zeros(inc(N),inc(N));
  U=zeros(N+1,N+1);
  for k=[0:N] do
    // U(inc(N),inc(k)) = f(x_0 * u^k * d^(N-k));
    U(N+1,k+1) = f(x_0 * u^k * d^(N-k));
  end;

  for n=[N-1:-1:0] do // le temps decroit de N-1 a 0
    //QUESTION: U(n+1,1:n+1)=RECOPIER L'EQUATION DU COURS;// ATTENTION AU DECALAGE DE 1
    U(n+1,1:n+1) = p*U(n+2,2:n+2) + (1-p)*U(n+2,1:n+1);
    // Une version "vectorisée" qui fait la même chose que
    //     for k=[0:n] do 
    //       U(inc(n),inc(k)) = ...
    //                p*U(inc(n+1),inc(k+1))+(1-p)*U(inc(n+1),inc(k)); 
    //     end;
  end;
  // res=U(inc(0),inc(0));
  res=U(1,1);
endfunction;
//-------------------------------------------------------
function main_7()
  N=10;
  sigma=0.3;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  K=100;x_0=100;

  prix(x_0,N)

  // Les deux algos font ils le même chose ?
  // on verifie : prix_slow(x_0,N) \approx prix(x_0,N)
  printf('différence entre les 2 resultats: %e\n', ...
          abs(prix_slow(x_0,N) - prix(x_0,N)));
endfunction
//-------------------------------------------------------
function main_8()
// Avec cet algorithme on peut augmenter N
// mais il faut renormaliser convenablement u et d pour
// rester borné.
// Essayer avec N=10,100,200,...,1000
  sigma=0.6;
  couleur=1

  plot2d([50:150],max([50:150]-K,0));//xclick;

  for N=[3,5,10,20,50,100,200] do
    d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);

    n=0;
    for x=[50:150] do
      n=n+1;
      courbe(n)=prix(x,N);
    end

    couleur=couleur+1;
    plot2d([50:150],courbe,style=couleur);
    xclick();
  end 
  // Ca converge, mais vers quoi ?
endfunction
//-------------------------------------------------------
function res=f_moy(x,s); res=max((s/N)-K,0);endfunction

function res=prix_moyenne(x,s,k,N)
  if k==N then res=f_moy(x,s); return; end;
  //QUESTION: res = <À COMPLÉTER>;
  res = p * prix_moyenne(x*u,s+x*u,k+1,N)+ (1-p) * prix_moyenne(x*d,s+x*d,k+1,N);
endfunction

function res=prix_slow_moyenne(x,N)
  res=prix_moyenne(x,x,0,N);
endfunction
//-------------------------------------------------------
function main_9()
  N=10;sigma=0.3;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);
  x_0=100;K=100;

  // Ca marche mais ce n'est pas très efficace ...
  printf('Prix option sur moyenne: %f\n',prix_slow_moyenne(x_0,N));
endfunction
//-------------------------------------------------------
function liste=liste_moyenne_rec(x,s,k,N)
// On constitue la liste des points visités par la chaine
// Si un point est visité deux fois, il y figure 2 fois.
  if k==N then 
    liste=[x;s];    //printf("%d: %f %f\n",x,s);
    return; 
  end;
  //QUESTION: liste = <À COMPLÉTER>;
  liste = [liste_moyenne_rec(x*u,s+x*u,k+1,N),liste_moyenne_rec(x*d,s+x*d,k+1,N)];
endfunction

function liste=liste_moyenne(x,N)
  // On part de (x,s=x) a l'instant 0
  liste=liste_moyenne_rec(x,x,0,N)
endfunction
//-------------------------------------------------------
function main_10()
  N=10;
  x_0=100;
  sigma=0.3;
  p=1/2;d=1-sigma/sqrt(N);u=1+sigma/sqrt(N);

  liste=liste_moyenne(x_0,N);

  // Tri des points selon les valeurs de la somme.
  // Les valeurs de x peuvent etre egales, mais pas celle de s.
  // Nous allons le verifier.
  [y,p]=sort(liste(2,:));
  liste1=liste(:,p);

  // On regarde si tous les points sont differents
  // en parcourrant le tableau ainsi classé
  epsilon=0.001;
  Taille=size(liste);
  match=[];
  for i=[1:Taille(2)-1] do
    if (norm(liste1(:,i)-liste1(:,i+1)) < epsilon) then 
      printf("Warning: (%f,%f) ~ (%f,%f)\n",...
                  liste1(1,i),liste1(2,i),liste1(1,i+1),liste1(2,i+1));
      match=[match,[liste1(1,i),liste1(2,i),liste1(1,i+1),liste1(2,i+1)]];    
    end;
  end;
  if size(match,'*') == 0 then 
     printf("Aucun point n''est dupliqué.\n");
  end
endfunction
//-------------------------------------------------------