
close all; clear all;
rng('default')

r=0.05; sig=0.2; T=1; S0=110; K=100;

%
% first, estimate variance / covariance
%

N = 1000;
Y = randn(1,N);       % Normal random variable
S = S0*exp((r-sig^2/2)*T + sig*sqrt(T)*Y);
F = exp(-r*T)*max(0,S-K);
C = exp(-r*T)*S;

Fave = sum(F)/N;
Cave = sum(C)/N;

lam = sum((F-Fave).*(C-Cave)) / sum((C-Cave).^2);

%
% now use control variate
%

N = 1e5;
Y = randn(1,N);       % Normal random variable
S = S0*exp((r-sig^2/2)*T + sig*sqrt(T)*Y);
F = exp(-r*T)*max(0,S-K);
C = exp(-r*T)*S;

F2 = F - lam*(C-S0);

Fave  = sum(F)/N;
F2ave = sum(F2)/N;

sd  = sqrt( sum((F -Fave ).^2)/(N*(N-1)) );
sd2 = sqrt( sum((F2-F2ave).^2)/(N*(N-1)) );

disp(sprintf('estimated price (without CV) = %f +/- %f', ...
     Fave, 3*sd));
disp(sprintf('estimated price (with CV)    = %f +/- %f', ...
     F2ave, 3*sd2));
disp(sprintf('exact price                  = %f', ...
     european_call(r,sig,T,S0,K,'value')));

