function [res]=gamma(n,C,beta)
  res= C / (n^beta);
endfunction

function [res]=F(x,g,proba)
   if (abs(g) < x) then res=1-proba; 
      else res= - proba; end;
endfunction

function X=estimation(n,C,beta)
  proba=0.95;
  X=zeros(1,n);
  G=grand(1,n,"nor",0,1);
  X(1) = 1.0;
  for k=[1:n-1] do
	X(k+1)=X(k)-gamma(k,C,beta)*F(X(k),G(k),proba);
  end;
endfunction

n=10000;

// Beaucoup de variance, la condition initiale est vite oubliee 
C=1;
beta=0.5;
X=estimation(n,C,beta);
plot2d([1:n],X);

// Moins de variance, la condition initiale est moins vite oubliee 
C=1;
beta=0.8;
X=estimation(n,C,beta);
plot2d([1:n],X);


//   Pas de variance, mais la condition initiale 
//   a beaucoup de mal a etre oubliee !
C=0.05;
beta=1;
X=estimation(n,C,beta);
plot2d([1:n],X);

//   Ca va un peu mieux
C=0.5;
beta=1;
X=estimation(n,C,beta);
plot2d([1:n],X);

//   encore mieux
C=5.0;
beta=1;
X=estimation(n,C,beta);
plot2d([1:n],X);