%A Monte Carlo simulation for a European call
%S0: the stock price at time zero
%K: the strike price, T: the time to maturity in years
%sigma: the volatility of the annualized continuously 
%       compounded rate of return of the stock
%r: the risk-free interest rate
%M: the number of realized paths
%N: the number of time steps
%dt: the time step length
%the equation is: dS = r*S*dt+sigma*S*dW

%S0 = 50;T=0.5;sigma=0.2;r=0.07;

function SE = geopath(S0, T, sigma, r, N,M);

startTime = cputime;

dt = T/N;
rdt1 = r*dt+1;
sigmadt = sigma*sqrt(dt);
ssdt = sigma^2*dt/2;
time = [0:dt:T];

randn('state',100); %given an initial state
W = randn(M,N); % an Mx1 matrix containing normal random numbers
WS = W.^2;
WE = repmat(rdt1,M,N)+sigmadt*W;
WM = WE+ssdt*(WS-repmat(1,M,N));
WA = cumsum(W,2);
tmp = exp((r-sigma^2/2)*time(2:N+1));
SA = S0*repmat(tmp,M,1).*exp(sigmadt*WA);
WE=S0*cumprod(WE,2);
WM=S0*cumprod(WM,2);

S0 = repmat(S0,M,1);
SA = [S0,SA];
SE = [S0,WE];
SM = [S0,WM];

%mean error
ErrorEuler = mean(abs(SE-SA))
ErrorMilstein = mean(abs(SM-SA))

plot(time,SE,'g',time,SM,'r',time,SA,'b');
xlabel('Time');
ylabel('S');

figure(2);
plot(time, ErrorEuler,'g',time, ErrorMilstein, 'r');
xlabel('Time');
ylabel('Expected absolute error');

totalTime = cputime -startTime