clear;

%% Data

rng(2022);

lambda = 2;
mu = 1;
X0 = 1;

T = 1;
nSteps = 2^8;
dt = 1/nSteps;

%% Exact solution

dW = sqrt(dt)*randn(1, nSteps);
W = [0, cumsum(dW, 2)];
tVec = 0:dt:T;
X = X0*exp((lambda - (mu^2)/2)*tVec + mu*W); 

%% Euler Maruyama solution

tRatio_vec = 2.^[2, 4];

figure;
plot(tVec, X, 'DisplayName', 'Exact solution', 'Linewidth', 2);
hold on;
for i = 1:length(tRatio_vec)
    tRatio = tRatio_vec(i);
    DeltaT = tRatio*dt;
    X_EM = nan(1, 1 + nSteps/tRatio);
    X_EM(1) = 1;
    dW_EM = W(:, 1+tRatio:tRatio:end) - W(:, 1:tRatio:end-tRatio);
    X_EM(:, 2:end) = X0*cumprod(1 + lambda*DeltaT + mu*dW_EM);
    plot(tVec(1:tRatio:end), X_EM, 'DisplayName', ['{\Delta}t = ', num2str(DeltaT)], 'Linewidth', 2);
end
xlabel('t', 'Fontsize', 20);
legend('Location', 'best', 'Fontsize', 20);
set(gca, 'FontSize', 20);
%saveas(gcf, 'figures/exercise2.png');