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

function [res] = BS_Call(S_0,K,sigma,r,T)
  d1=(log(S_0/K)+(r+sigma^2/2)*T)/(sigma*sqrt(T));
  d2=d1-sigma*sqrt(T);
  res=S_0*erf(d1)-K*exp(-r*T)*erf(d2);
endfunction

function [y] = black_scholes(lambda, Mean,var)
  m = Mean + log(lambda);
  d = (m-log(K))/sqrt(var);
  y= exp(m + (var/2)) * erf(d+sqrt(var)) - K*erf(d)
endfunction

function [I_T,Z_T] = ...
	  control_variate_simulation(r,sigma,Sigma,S_0,a)
  // Simulation of I_T and Z_T the control variate
  n=size(Sigma);n=n(2);
  d=prod(size(S_0));

  // Simulation of the basket price I_T
  W_T = T * diag(r - sigma^2/2) * ones(d,N) ...
			 + sqrt(T) * Sigma * rand(n,N,"gauss");
  S_T = diag(S_0) * exp(W_T);
  I_T= a' * S_T;

  // Simulation of the control variate
  H_0 = a .* S_0 / I_0; 
              // = a_i S_0^i / I_0
  Z_T = H_0' * W_T;
endfunction

function [] = test_call_basket(d,N)
  // definition of the model
  r=0.05;
  T = 1;
  rho=0.5;
  sigma=0.3*ones(d,1);
  Rho = (1 - rho) *eye(d,d) + rho * ones(d,d);
  Gamma = diag(sigma) * Rho * diag(sigma);
  Sigma=sqroot(Gamma);

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

  I_0 = a' * S_0;  // at the money
  K = 1.0 * I_0;
  // K = 1.2 * I_0; // out of the money
  // K = 0.8 * I_0; // in the money

  [I_T,Z_T] = control_variate_simulation(r,sigma,Sigma,S_0,a);
  
  // lambda exp(Z_T) approximate S_T
  // we compute lambda, the mean and the variance of 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; 
  Mean= T * (H_0' * (r-sigma^2/2)); // Mean of Z_T
  var = sum(diag(J_0) * Rho * diag(J_0)); // variance of Z_T
  
  // The Monte-Carlo method without variance reduction
  Y=exp(-r*T) * max(I_T-K,0);
  estimation= mean(Y);
  ecart_type=st_deviation(Y); 
  method_error=1.96*ecart_type/sqrt(N); 
  printf("Without variance reduction N=%d, %f +- %f\n",N, estimation, method_error);

  // The Monte-Carlo method with variance reduction
  Y=exp(-r*T) * (max(I_T-K,0) ...
                    - max(lambda * exp(Z_T) - K,0));
                      // here is the control variate

  estimation= mean(Y) + exp(-r*T) * black_scholes(lambda,Mean,var);
  ecart_type=st_deviation(Y); 
  method_error=1.96*ecart_type/sqrt(N); 
  printf("With variance reduction N=%d, %f +- %f\n",N, estimation, method_error);
endfunction 

stacksize(10^8);
d=10;
N=1000;
test_call_basket(d,N)