% Antithetic technique for variance reduction for Black-Scholes.
% Uses Euler discretization for the terminal stock price S(T).
clc; clear;

% Inline function for the Black Scholes call
C = inline('s*exp(-q*T)*normcdf((log(s/K) + (r-q+v^2/2)*T)/v/sqrt(T)) - K*exp(-r*T)*normcdf((log(s/K) + (r-q+v^2/2)*T)/v/sqrt(T) - v*sqrt(T))',...
	's','K','r','q','v','T');

% Inline function for the Black Scholes put
P = inline('K*exp(-r*T)*normcdf(-(log(s/K) + (r-q+v^2/2)*T)/v/sqrt(T) + v*sqrt(T)) - s*exp(-q*T)*normcdf(-(log(s/K) + (r-q+v^2/2)*T)/v/sqrt(T))',...
	's','K','r','q','v','T');

S0  = 100;        % Spot price
K   = 100;        % Strike price
r   = 0.05;       % Risk free rate
q   = 0.0;        % Dividend yield
v   = 0.25;       % Volatility
T   = 0.5;        % Time to maturity
PutCall = 'C';    % 'P'ut or 'C'all
N = 20000;        % Number of simulations.

% Initialize the terminal stock price matrices.
ST  = zeros(N,1);
ST1 = zeros(N,1);

% Simulate the stock price and the antithetic stock price.
for i=1:N
	e = randn(1);
	ST(i)  = S0*exp((r-q-v^2/2)*T + v*sqrt(T)*e);
	ST1(i) = S0*exp((r-q-v^2/2)*T + v*sqrt(T)*(-e));
end

% Obtain the option prices.
if strcmp(PutCall,'C')
	BSPrice = C(S0,K,r,q,v,T)
	SimPrice = exp(-r*T)*mean(max(ST-K,0))
	AntiPrice = exp(-r*T)*mean(0.5*(max(ST-K,0) + max(ST1-K,0)))
else
	BSPrice = P(S0,K,r,q,v,T)
	SimPrice = exp(-r*T)*mean(max(K-ST,0))
	AntiPrice = exp(-r*T)*mean(0.5*(max(K-ST,0) + max(K-ST1,0)))
end

% Calculate the errors.
AntiError = abs(AntiPrice - BSPrice)/BSPrice * 100
SimError  = abs(SimPrice - BSPrice)/BSPrice * 100
		