%% basket option
% create two sets of asset prices
S0 = [2 2];
r = 0.01;
sigma = [.5 .2];
strike = 2;
T = 1;
M = 1e4; % # repeats
Z = randn(M,2); % uncorrelated normals
sigma = repmat(sigma,[M,1]);
S = exp( (r-sigma.^2/2)*T + sqrt(T)*sigma.*Z ) * diag(S0);
DiscountedPayoff = exp(-r*T)*max(S(:,1)+S(:,2)-strike,0);
V = mean(DiscountedPayoff)
e = 1.96*std(DiscountedPayoff)/sqrt(M)

%% Control Variates
Y = exp(-r*T)*max(S(:,1)-strike/2,0) + exp(-r*T)*max(S(:,2)-strike/2,0);
eY = bs(S0(1),strike/2,r,sigma(1,1),T,0,'call')...
    + bs(S0(2),strike/2,r,sigma(1,2),T,0,'call');
% The expected value of Y-eY is zero.
c = corrcoef(DiscountedPayoff,Y); c = c(2,1)
cv = cov(DiscountedPayoff,Y); 
theta = cv(2,1)/cv(2,2);
X = DiscountedPayoff+theta*(eY-Y);
VC = mean(X)
eC = 1.96*std(X)/sqrt(M)
%% Other payoffs
Maxoption = exp(-r*T)*max(max(S(:,1),S(:,2))-strike,0);
control = exp(-r*T)*max(S(:,1)-strike,0);
cv = cov(Maxoption,control);
theta = cv(2,1)/cv(2,2);
V = mean(Maxoption), e = std(Maxoption)*1.96/sqrt(M)
X = Maxoption + theta*(bs(S0(1),strike,r,sigma(1),T,0,'call')-control);
VC = mean(X)
eC = 1.96*std(X)/sqrt(M)