clear;

%% Data

rng(2022);

X0 = 1;

T = 1;
DeltaTs = 2.^(-8:-2);
nTrajectories = 10^4;

%% Exact solution

DeltaT_min = DeltaTs(1);
nSteps_min = T/DeltaT_min;

dW1 = sqrt(DeltaT_min)*randn(nTrajectories, nSteps_min);
W1 = [zeros(nTrajectories, 1), cumsum(dW1, 2)];
dW2 = sqrt(DeltaT_min)*randn(nTrajectories, nSteps_min);
W2 = [zeros(nTrajectories, 1), cumsum(dW2, 2)];

X = exp(-3*T/2 + W1(:,end) + W2(:,end));

%% Errors

strongError = nan(1, length(DeltaTs));

for i = 1:length(DeltaTs)
    
    DeltaT = DeltaTs(i);
    ratio = DeltaT/DeltaT_min;
    nSteps = T/DeltaT;
    
    dW1_num = W1(:, 1+ratio:ratio:end) - W1(:, 1:ratio:end-ratio);
    dW2_num = W2(:, 1+ratio:ratio:end) - W2(:, 1:ratio:end-ratio);
    X_MP = X0*prod(1 - DeltaT/2 + dW1_num + dW2_num + (dW1_num.^2 - DeltaT)/2 + (dW2_num.^2 - DeltaT)/2 + dW1_num.*dW2_num, 2);
    
    strongError(i) = mean(abs(X_MP - X));

end

%% Figures

figure;
loglog(DeltaTs, strongError, 's:', 'Linewidth', 2);
hold on;
loglog(DeltaTs, DeltaTs/5, '--', 'Linewidth', 2);
xlabel('{\Delta}t', 'Fontsize', 20);
xticks([1e-2, 1e-1]);
legend('Milstein-Platen', 'Order 1', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
set(gca, 'FontSize', 20);
saveas(gcf, 'figures/exercise7.png');