// We want to sample S_T
// For efficience reason, we assume that Sigma 
// has already been computed 

// Sigma is the matrix of the BS model

// sigma is the vector of volatility of the assets

function [S_T] = black_scholes(N,T,S_0, r,sigma, Sigma)
  // Sigma is the matrix of the BS model
  // sigma is the vector of volatility of the assets
  // N is the number of drawings

  // d is the size of the basket
  d=prod(size(S_0));
  n=size(Sigma);n=n(2);
  // A little bit hard to understand ! Sorry.
  // I want that my procedure works even when Rho is non invertible
  // (this is the case when, for instance, rho = 1).
  // The matrix is then degenerate and scilab give you a matrix dxn 
  // instead of dxd (n is the rank of the matrix)

  // After that, it is easy
  S_T =  exp(T * diag(r - sigma^2/2) * ones(d,N) ...
			 + sqrt(T) * Sigma * rand(n,N,"gauss"));
  S_T = diag(S_0) * S_T;
endfunction

N=10;
T=1;
r=0.05;

d=3; // or d = 10;
rho=0.5; 

// The initial values of each assets is equal to 100
S_0=100*ones(d,1);

//  Rho,  1 on the diagonal,  rho outside
Rho = (1 - rho) *eye(d,d) + rho * ones(d,d);

// The volatility of each assets is equal to 0.3
sigma=0.3*ones(d,1);

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

black_scholes(N,T,S_0, r,sigma, Sigma)'