% Monte Carlo simulations to price European call options and to study two
% variance reduction techniques.
clear;
% The following specifies the parameters in the option:
r = 0.01;         % risk-free interest rate
sigma = 0.25;     % volatility
T = 0.5;          % expiration
S0 = 100;         % current stock price
K = 100;          % strike of the option

p = 100;           % number of points in the convergence study

% The following is the exact Black-Scholes price:
Vbs = bsoption(S0,K,r,0,sigma,T,'call');  

Serror = zeros(p,1);
Serror_at = zeros(p,1);
Serror_mm = zeros(p,1);

for k=1:p
    M = 100*k^2;       % number of paths
    xi = randn(M,1);   % simulated normal distribution draws
    xstd = sqrt(var(xi)); % standard deviation of the draws

    ST1 = S0*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*xi);
    % simulated prices without any variance reduction
    ST2 = S0*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*(-xi));  
    % image prices of ST1
    ST3 = S0*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*xi/xstd);
    % prices using moment matching draws
    ST4 = S0*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*(-xi)/xstd);
    % image prices of ST3
    
    % the following are the corresponding payoffs
    VT1 = 0.5*((ST1 - K)+abs(ST1-K));
    VT2 = 0.5*((ST2 - K)+abs(ST2-K));
    VT3 = 0.5*((ST3 - K)+abs(ST3-K));
    VT4 = 0.5*((ST4 - K)+abs(ST4-K));

    V0 = exp(-r*T)*mean(VT1);
    % simulated result without variance reduction
    V0at = exp(-r*T)*0.5*(mean(VT1)+mean(VT2));
    % result with anti-thetic technique
    V0mm = exp(-r*T)*0.5*(mean(VT3)+mean(VT4));
    % result with second moment matching
    
    % the following are errors for all three cases:
    Serror(k) = abs(V0-Vbs);
    Serror_at(k) = abs(V0at-Vbs);
    Serror_mm(k) = abs(V0mm-Vbs);
end

plot(10:p,Serror(10:p),'+',10:p,Serror_at(10:p),'x',10:p,Serror_mm(10:p),'o');
xlabel('sqrt(M/1000)'); ylabel('Error'); title('MC convergence');
legend('w/o variance reduction','anti-thetic','2nd moment matching');

