clear;

rng(500);
mu = -1; sigma = 0.3; Xzero = 1; T = 2;
M = 2e5;
X0 =1;

Ximpl_avg1 = zeros(5,1);
Ximpl_avg2 = zeros(5,1);

for p = 1:6
    Dt = 2^(p-11); N = round(T/Dt); % N steps of size Dt
    rn = randi(2, M, N);
    xi = sqrt(Dt)*(rn==1) - sqrt(Dt)*(rn==2);
    Ximpl1 = X0*(prod((1./(1-Dt*mu -sigma*xi)),2)); % implicit
    Ximpl2 = X0*(prod((1./(1-Dt*(mu-sigma^2)-sigma*xi)),2)); % implicit with adjusted drift
    Ximpl_avg1(p) = mean(Ximpl1);
    Ximpl_avg2(p) = mean(Ximpl2);
end
Xerr1 = abs(Ximpl_avg1 - exp(mu*T));
Xerr2 = abs(Ximpl_avg2 - exp(mu*T));

figure;
Dtvals = 2.^([1:6]-11);
loglog(Dtvals,Xerr1,'b*-', 'Linewidth', 2), 
hold on
loglog(Dtvals,Xerr2,'g*-', 'Linewidth', 2), 
hold on
loglog(Dtvals,0.5*Dtvals,'r--', 'Linewidth', 2), 
hold off
xlabel('\Delta t', 'Fontsize', 20);
ylabel('|E[X_T]-E[X_N]|', 'interpreter','tex','Fontsize', 20);
legend('Sol.','â-Sol','\Delta t')
set(gca, 'FontSize', 20);
saveas(gcf, 'figures/td_10.png');
