% Pricing a max option on an asset that obeys GBM 
% using both EM and the exact GBM formula to get prices
% (the error from taking the "max" over a discrete set
%  of points is more significant than the EM error) 

r = 0.02;
mu = 0.02;
sigma = 0.2;
T = 2; 
S0 = 1;
K = 1.2;

M = 1000000;     % keep large as want to see the deterministic error 
N = 50;          % vary to see the deterministic error 
h = T/N;
t = h*(0:N);

S = S0*ones(M,N+1);
Z = randn(M,N); 

% Euler Maruyama 
%for i=1:N
%    S(:,i+1) = S(:,i).*(1 + h * r  + sigma * sqrt(h) * Z(:,i));
%end

% GBM 
for i=1:N
    S(:,i+1) = S(:,i).*exp( (r-sigma^2/2)*h + sigma * sqrt(h)*Z(:,i));
end

MM = max(S,[],2)-K;
NN = (MM > 0).*MM*exp(-mu*T) ; 
[mean(NN), std(NN)/sqrt(M)]