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(N,K, fonction_payoff)
  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) * fonction_payoff(a,S_T,K);

  estimation=mean(payoff);          // the Monte-Carlo estimation
  standard_deviation=st_deviation(payoff);  // the estimation of the standard deviation
  method_error=1.96*standard_deviation/sqrt(N);   // half width of the confidence interval

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

stacksize(1000000);

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


K= a'*S_0; // at the money option 
test(10000,K,payoff_call_basket);
test(10000,K,payoff_put_basket);
// the variance is smaller for the put than for the call

K=0.8*a'*S_0; // in the money option 
test(10000,K,payoff_call_basket);
test(10000,K,payoff_put_basket);
// the variance is smaller for the put than for the call

K=1.2*a'*S_0;// out of the money option 
test(10000,K,payoff_call_basket);
test(10000,K,payoff_put_basket);
// the variance is smaller for the call than for the put

K=1.5*a'*S_0;// deep out of the money option
test(10000,K,payoff_call_basket);
test(10000,K,payoff_put_basket);
// the variance is smaller for the call than for the put