% Check the strong convergence rate of Euler-Maruyama method

% Wiener part
lambda=2; mu=1; Xzero=1;
T=1; N=2^9; dt=T/N;
M=1000; % number of Wiener paths

err=zeros(M,5); % preallocation of the matrix containing the values
                % of the errors in the final integration point
                % for each computed path and in correspondence
                % of 5 different values for the stepsize of EM method
for s=1:M
    dW=sqrt(dt)*randn(1,N);
    W=cumsum(dW);
    Xtrue=Xzero*exp((lambda-mu^2/2)*T+mu*W(end));
    for p=1:5
        R=2^(p-1); h=R*dt; L=N/R;
        Xtemp=Xzero;
        for j=1:L
            Winc=sum(dW((j-1)*R+1:j*R));
            Xtemp=Xtemp+h*lambda*Xtemp+mu*Xtemp*Winc;
        end
        err(s,p)=abs(Xtemp-Xtrue); % storing the error in T=1 and 
                                   % and in correspondence to the 
                                   % current path
    end
end

% Pattern of the average error in the final integration point
hvec=dt*(2.^(0:4));
loglog(hvec,mean(err),'b*-')
hold on
loglog(hvec,hvec.^0.5,'r--') % reference slope for order 1/2