r = 0.1;
sig11 = 0.3;
sig12 = 0.2;     % change later to -0.4
sig21 = 0.1;     % change later to -0.1
sig22 = 0.4;     
S10 = 1;
S20 = 1;
T = 1 ;
B1 = 0.9;        % S1 must stay over this 
B2 = 1.1;        % S2 must stay under this 

N = 100;         % number of subintervals - needs to be justified!
M = 100000;      % number of simulations  -     "       " 
h = T/N ;

S1 = S10*ones(M,N+1);
S2 = S20*ones(M,N+1);
Z1 = randn(M,N);
Z2 = randn(M,N); 

a1 = h*(r - (sig11^2+sig12^2)/2) ; 
a2 = h*(r - (sig21^2+sig22^2)/2) ; 
sqr = sqrt(h);
% generate prices using GBM rules 
for i=1:N
   S1(:,i+1) = S1(:,i).*exp(a1 + sqr*(sig11*Z1(:,i)+sig12*Z2(:,i))) ; 
   S2(:,i+1) = S2(:,i).*exp(a2 + sqr*(sig21*Z1(:,i)+sig22*Z2(:,i))) ; 
end

val = exp(-r*T) * (S1(:,N+1)+S2(:,N+1)) .* (min(S1,[],2)>B1) .* (max(S2,[],2)<B2); 

mean(val)
std(val)/sqrt(M) 

% look also at the number of simulations in which barriers are hit: 
sum( min(S1,[],2)<B1 )                            % number of times B1 hit by S1 
sum( max(S2,[],2)>B2 )                            % number of times B2 hit by S2 
sum( (max(S2,[],2)>B2) .* (min(S1,[],2)<B1) )     % number of times both hit
