function [y]=NN(x)
  [y,Q]=cdfnor("PQ",x,0,1);
endfunction

function [y] = formule_black_scholes(m,var)
  d = (m-log(K))/sqrt(var);
  y= exp(m + (var/2)) * NN(d+sqrt(var)) - K*NN(d)
endfunction

function [I_T,Z_T] = simulation_variable_de_controle(r,sigma,Sigma,S_0,a)
  // Simulation de 
  //   - I_T la valeur du panier
  //   - Z_T la variable de controle
  n=size(Sigma);n=n(2);
  d=prod(size(S_0));
  LOG_S_T_sur_S_0 = T * diag(r - sigma.^2/2) * ones(d,N) ...
			              + sqrt(T) * Sigma * rand(n,N,"gauss");
  S_T = diag(S_0) * exp(LOG_S_T_sur_S_0);

  I_T= a' * S_T; // Calcul de la valeur de l'indice
  H_0 = a .* S_0 / I_0; // = a_i S_0^i / I_0
  Z_T = ****** ... *******;// A COMPLETER  Valeur de la variable de controle
endfunction

function [] = test_call_basket(d,N,K_sur_I0)

  // Les données du problème
  r=0.05;
  T = 1;
  rho=0.5;
  sigma=0.3*ones(d,1);
  S_0=100*ones(d,1);
  Rho = (1 - rho) *eye(d,d) + rho * ones(d,d);
  a=(1/d)*ones(d,1);

  Gamma = diag(sigma) * Rho * diag(sigma);
  Sigma=sqroot(Gamma);


  I_0 = a' * S_0;
  K = K_sur_I0 * I_0;

  [I_T,Z_T] = simulation_variable_de_controle(r,sigma,Sigma,S_0,a);
  
  // lambda exp(Z_T) est une approximation de S_T
  // on calcule lambda, la moyenne et la variance de Z_T
  H_0 = a .* S_0 / I_0; // = a_i S_0^i / I_0
  J_0 = H_0 .* sigma;   // = a_i S_0^i sigma_i / I_0 

  lambda = I_0; 
  moyenne= T * (H_0' * (r-sigma.^2/2)); // moyenne de Z_T
  var = sum(diag(J_0) * Rho * diag(J_0)); // variance de Z_T
  
  // Methode de Monte Carlo sans reduction de variance
  Y=exp(-r*T) * max(I_T-K,0);
  estimation= mean(Y);
  ecart_type=stdev(Y); // estimation de l'ecart type
  erreur=1.96*ecart_type/sqrt(N); // demi-largeur de l'intervalle de confiance
  printf("Sans réduction de variance : N=%d, %f +- %f\n",N, estimation, erreur);

  
  // Methode de Monte Carlo avec reduction de variance
  Y=exp(-r*T) * (max(I_T-K,0) - max(lambda * exp(Z_T) - K,0));
  estimation= mean(Y) + exp(-r*T) * formule_black_scholes(moyenne+log(lambda),var);
  ecart_type=stdev(Y); // estimation de l'ecart type
  erreur=1.96*ecart_type/sqrt(N); // demi-largeur de l'intervalle de confiance
  printf("Avec réduction de variance : N=%d, %f +- %f\n",N, estimation, erreur);

endfunction 

stacksize(10^8);

d=10;
N=10000;

K_sur_I0=1.2;test_call_basket(d,N,K_sur_I0);
K_sur_I0=1.0;test_call_basket(d,N,K_sur_I0);
K_sur_I0=0.8;test_call_basket(d,N,K_sur_I0);