%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;K=50;T=0.5;sigma=0.2;r=0.07;

function [call,exactCall,callError] = MilsteinMC(S0, K, T, sigma, r, N, M);

startTime = cputime;

dt = T/N;
rdt1 = r*dt+1;
sigmadt = sigma*sqrt(dt);
ssdt = sigma^2*dt/2;
randn('state',100); %given an initial state
M0 = M/10;
call = 0;
for i=1:10
    W = randn(M0,N); % an Mx1 matrix containing normal random numbers
    WS = W.^2;
    W1 = repmat(rdt1,M0,N)+sigmadt*W+ssdt*(WS-repmat(1,M0,N));
    W2=cumprod(W1,2);
    ST = S0*W2(1:M0,N);
    callPayoff = max(ST-repmat(K,M0,1),repmat(0,M0,1));
    call = call + exp(-r*T)*sum(callPayoff);
end
call = call/M;
exactCall = blsprice(S0,K,r,T,sigma,0);
callError = abs(call-exactCall);

totalTime = cputime -startTime