function [Y] = payoff_call_basket(a,S_T,K)
  // S_T is a d*N matrix 
  //   d=number of assets
  //   N=number of drawings
  // I_T and Y are 1*N vector
  I_T = a' * S_T;
  Y   = max(I_T-K,0);
endfunction


function [Y] = payoff_put_basket(a,S_T,K)
  // S_T is a d*N matrix
  // I_T and Y are 1*N vector
  I_T= a' * S_T;
  Y=max(K-I_T,0);
endfunction

function []=test_call_wcv(N,K)
  // without control variate
  r=0.05;
  T = 1;
  d=10;
  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);
  Gamma = diag(sigma) * Rho * diag(sigma);
  Sigma=sqroot(Gamma);
  a=(1/d)*ones(d,1);
  I_0= a'*S_0;
  
  S_T=black_scholes(N,T,S_0, r,sigma, Sigma);
  payoff=exp(-r*T) * payoff_call_basket(a,S_T,K);

  estimation=mean(payoff);
  ecart_type=st_deviation(payoff); 
  erreur=1.96*ecart_type/sqrt(N); 

  printf("Direct N=%d, %f +- %f\n",N, estimation, erreur);
endfunction

function []=test_call_cv(N,K)
  // using control variate
  r=0.05;
  T = 1;
  d=10;
  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);
  Gamma = diag(sigma) * Rho * diag(sigma);
  Sigma=sqroot(Gamma);
  a=(1/d)*ones(d,1);
  I_0= a'*S_0;
  
  S_T=black_scholes(N,T,S_0, r,sigma, Sigma);
  I_T=a' * S_T;


  // Here we use the control variate I_T
  payoff=exp(-r*T) * (payoff_call_basket(a,S_T,K) - a' * S_T);

  estimation=mean(payoff) + I_0;
  ecart_type=st_deviation(payoff); 
  erreur=1.96*ecart_type/sqrt(N); 

  printf("Direct N=%d, %f +- %f\n",N, estimation, erreur);
endfunction

stacksize(1000000);

d=10;
a=(1/d)*ones(d,1);
S_0=100*ones(d,1);

K= 1.0 * a'*S_0; // at the money option
test_call_wcv(10000,K);
test_call_cv(10000,K); // in this case the control variate works but is not very efficient

K= 1.2 * a'*S_0; // out of the money option
test_call_wcv(10000,K);
test_call_cv(10000,K); // the control variate increase the variance ! 

K= 0.8 * a'*S_0; // in the money option
test_call_wcv(10000,K);
test_call_cv(10000,K); // in this case the control variate is efficient

K= 0.5 * a'*S_0;
test_call_wcv(10000,K);
test_call_cv(10000,K); // the smaller K is and the better the control variate is