function [ms,mil] = milstein(mu,sigma,dt)
%compares Euler and Milstein scheme for geometric BM , dt=size of step.
n=dt*10000;
m=1/dt;
w=cumsum(normrnd(0,1/100,1,10000));
y=exp(((mu-.5*sigma^2)*.0001*(1:10000)) + sigma*w);
w1 = [0,w(1,n*(1:m))];
dw=diff(w1);
eul = 1*ones(1,m+1);
mil = 1*ones(1,m+1);
for i=1:m
eul(1,i+1) = eul(1,i) + mu*dt*eul(1,i) + sigma*eul(1,i)*dw(1,i);
mil(1,i+1) = mil(1,i) +mil(1,i)*(mu*dt + sigma*dw(1,i) + .5*sigma^2*((dw(1,i)^2 - dt)));
end
x=dt*[0,1:m];
z=.0001*(1:10000);
plot(z,y,x,eul,'r',x,mil,'--k')
grid on
title(['comparison of Euler (red) and Milstein(dashed) schemes: dt= ' num2str(dt)])
%mean squared differences
y=[0 y(n:n:10000)];
ms=[mean((y-eul).^2),mean((y-mil).^2)];
%xlabel(['mean sqared differences: Euler= ' num2str(ms(1)) ' Milstein= ' num2str(ms(2))])