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

r=0.02;
sigma=0.3;
T=1;
x=100;
K=200;

function [res]=f(g)
   d=prod(size(g));
   S_T=zeros(1,d);
   S_T=x .* exp((r - sigma*sigma/2)*T*ones(1,d)+sigma*sqrt(T)*g);
   res=exp(-r*T)*max(S_T-K,0);
endfunction

function [res]=F(lambda,g)
   res= f(g)^2 * (lambda - g) * exp(-lambda*g+lambda^2/2);
endfunction

function X=estimation(n,C,beta)
  X=zeros(1,n);
  //G=grand(1,n,"nor",0,1);
  lambda0=0;
  //lambda0=(log(K/x)-(r-sigma^2/2)*T)/(sigma*sqrt(T));
  X(1) = lambda0;
  for k=[1:n-1] do
	X(k+1)=X(k)-gamma(k,C,beta)*F(X(k),G(k));
  end;
endfunction

function X=estimation_chen(n,C,beta)
  X=zeros(1,n);
  //G=grand(1,n,"nor",0,1);
  lambda0=0;
  //lambda0=(log(K/x)-(r-sigma^2/2)*T)/(sigma*sqrt(T));
  X(1) = lambda0;
  for k=[1:n-1] do
	X(k+1)=X(k)-gamma(k,C,beta)*F(X(k),G(k));
    if(abs(X(k+1)) >= 5) then X(k+1)=lambda0; end; 
  end;
endfunction

C=0.5;
beta=1.0;
n=10000;

G=grand(1,n,"nor",0,1);
X=estimation(n,C,beta);
Y=estimation_chen(n,C,beta);
plot2d([1:n],X);
plot2d([1:n],Y);