% Euler and Milstein discretization for Black-Scholes.

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');

% Option features.
S0  = 60;         % Spot price
K   = 50;         % Strike price
r   = 0.00272;    % Risk free rate
q   = 0.0;        % Dividend yield
v   = 0.347;      % Volatility
mat = 0.25;       % Time to maturity
PutCall = 'P';    % 'P'ut or 'C'all

% Simulation settings.
N = 5000;        % Number of simulations.
T = 100;         % Nuber of time steps.
dt = mat/T;      % Time increment

% Initialize the terminal stock price matrices 
% for the Euler and Milstein discretization schemes.
Sm = zeros(N,T);          % Milstein
SE = zeros(N,T);          % Euler stock price
Se = zeros(N,T);          % Euler log stock price
Sm(:,1) = S0;
SE(:,1) = S0;
Se(:,1) = S0;

% Simulate the stock price and the log-stock price under the Euler scheme
% Simulate the stock price under the Milstein schemes
for n=1:N
	for t=2:T
		Z = randn(1);
		W = sqrt(dt)*Z;
		% Euler discretization of stock price
		SE(n,t) = SE(n,t-1) + (r-q)*SE(n,t-1)*dt + v*SE(n,t-1)*W;
		% Euler discretization of log-stock price
		Se(n,t) = Se(n,t-1)*exp((r-q-v^2/2)*dt + v*W);
		% Milstein discretization
		Sm(n,t) = Sm(n,t-1) + (r-q)*Sm(n,t-1)*dt + v*Sm(n,t-1)*W...
  			     + 0.5*v^2*Sm(n,t-1)*(Z^2-1)*dt;
	end
end

% Terminal stock prices.
STE = SE(:,end);        % Euler stock price
STe = Se(:,end);        % Euler log stock price
STm = Sm(:,end);        % Milstein stock price


% Obtain the simulated option prices.
if strcmp(PutCall,'C')
	BSPrice  = C(S0,K,r,q,v,mat)
	EulerSt  = exp(-r*mat)*mean(max(STE-K,0))
	EulerLn  = exp(-r*mat)*mean(max(STe-K,0))
	Milstein = exp(-r*mat)*mean(max(STm-K,0))
else
	BSPrice = P(S0,K,r,q,v,mat)
	EulerSt  = exp(-r*mat)*mean(max(K-STE,0))
	EulerLn  = exp(-r*mat)*mean(max(K-STe,0))
	Milstein = exp(-r*mat)*mean(max(K-STm,0))
end

% Calculate the errors.
EulerErrorSt   = abs(EulerSt  - BSPrice)/BSPrice * 100
EulerErrorLn   = abs(EulerLn  - BSPrice)/BSPrice * 100
MilsteinError  = abs(Milstein - BSPrice)/BSPrice * 100

