clf
rng(500);
mu = 2; sigma = 0.1; Xzero = 1; T = 1;
M = 1e4;

Xem = zeros(5,1);

for p = 1:5   
    Dt = 2^(p-10); N = round(T/Dt); % N steps of size Dt
    Xtemp = Xzero*ones(M,1);
    for j = 1:N
       Winc = sqrt(Dt)*randn(M,1); % Brownian increment for each path
       Xtemp = Xtemp + Dt*mu*Xtemp + sigma*Winc; % EM
    end
    Xem(p) = mean(Xtemp);
   
end
Xerr = abs(Xem - exp(mu*N*Dt));
% true expected value is exp(mu*T)

%% Figures

figure;
Dtvals = 2.^([1:5]-10);
loglog(Dtvals,Xerr,'b*-', 'Linewidth', 2), 
hold on
loglog(Dtvals,2*Dtvals,'r--', 'Linewidth', 2), 
hold off
% reference slope of 1
axis([1e-3 1e-1 1e-3 1])
xlabel('\Delta t'), 
ylabel('| E(X(T)) - E (X_N) |'),
legend('Weak Convergence', 'Order 1', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
set(gca, 'FontSize', 20);
saveas(gcf, 'figures/EM_09.png');
