clear;

%% Data

rng(2022);

alpha = 0.2;

mu = 2;
sigma =1;
X0 = 1;
T = 1;

DeltaTs = [0.1,0.05,0.01,0.005,0.001,0.0005];
nTrajectories = 1;

%% Exact solution

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

dW = sqrt(DeltaT_min)*randn(nTrajectories, nSteps_min);
W = [zeros(nTrajectories, 1), cumsum(dW, 2)];

vecT = 0:DeltaT_min:T;

X_true = zeros(nTrajectories,nSteps_min+1);
X_true(:,1) = X0;

for t = 2:nSteps_min+1
  X_true(:,t) = X0*exp((mu-0.5*sigma^2)*vecT(t) + sigma*W(:,t));
end

%% Rate

error_sup = zeros(length(DeltaTs),1);
error_alpha_sup = zeros(length(DeltaTs),1);


for i = 1:length(DeltaTs)
    
    DeltaT = DeltaTs(i);
    ratio = DeltaT/DeltaT_min;
    nSteps = T/DeltaT;
    
    dW_num = W(:, 1+ratio:ratio:end) - W(:, 1:ratio:end-ratio);
   
    error = zeros(nSteps,1);
    error_alpha = zeros(nSteps,1);
    
    X_n = zeros(1,1+nSteps);
    X_n(:,1) = X0;

    for j = 1:nSteps
        X_n(:,j+1) = X_n(:,j).*(1 + DeltaT*mu + sigma*dW_num(:,j));
        error(j) = abs(X_n(:,j+1) - X_true(:,floor((j*ratio)+1)));
        error_alpha(j) = abs(X_n(:,j+1) - X_true(:,floor((j*ratio)+1)));
    end

    error_sup(i) = max(error);
    error_alpha_sup(i) = max(error_alpha);
end

%% Figures

figure;
loglog(sqrt(DeltaTs), error_sup, 's:', 'Linewidth', 2);
ylabel('Sup Error', 'Fontsize', 20)
xlabel('sqrt({\Delta}t)', 'Fontsize', 20);
legend('sup vs sqrt(Delta)', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
saveas(gcf, 'figures/EM_07_M_1.png');

figure;
loglog((DeltaTs).^(alpha), error_alpha_sup, 'x:', 'Linewidth', 2);
ylabel('Sup Error', 'Fontsize', 20)
xlabel('({\Delta}t)^{\alpha}', 'Fontsize', 20);
legend('sup vs ({\Delta}t)^{alpha}', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
saveas(gcf, 'figures/EM_07_M_1_alpha.png');

%%%%%%%%%

DeltaTs = [0.01,0.001];
nTrajectories = 10^3;

%% Exact solution

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

dW = sqrt(DeltaT_min)*randn(nTrajectories, nSteps_min);
W = [zeros(nTrajectories, 1), cumsum(dW, 2)];
vecT = 0:DeltaT_min:T;

X_true = zeros(nTrajectories, nSteps_min+1);
for t = 1:nSteps_min
  X_true(:,t) = X0.*exp((mu-0.5*sigma^2)*vecT(:,t) + sigma*W(:,t));
end

%% Rate

error_sup_M = zeros(nTrajectories, length(DeltaTs));
error_alpha_sup_M = zeros(nTrajectories, length(DeltaTs));


for i = 1:length(DeltaTs)
        
        DeltaT = DeltaTs(i);
        ratio = DeltaT/DeltaT_min;
        nSteps = T/DeltaT;
        
        dW_num = W(:, 1+ratio:ratio:end) - W(:, 1:ratio:end-ratio);
        X_n = zeros(nTrajectories,1+nSteps);
        X_n(:,1) = X0;
        
        error_M = zeros(nTrajectories, length(nSteps));
        error_alpha_M = zeros(nTrajectories, length(nSteps));
    
        for j = 1:nSteps
            X_n(:,j+1) = X_n(:,j).*(1 + DeltaT*mu + sigma*dW_num(:,j));
            error_M(:,j) = abs(X_n(:,j+1) - X_true(:,floor((j*ratio)+1)));
            error_alpha_M(:,j) = abs(X_n(:,j+1) - X_true(:,floor((j*ratio)+1)));
        end
       
        error_sup_M(:,i) = max(error_M,[],2)/(sqrt(DeltaT));
        error_alpha_sup_M(:,i) = max(error_alpha_M,[],2)/(DeltaT)^alpha;
end


%% Figures

figure;
histogram(error_sup_M(:,1));
saveas(gcf, 'figures/EM_7_M_1000_N_100.png');
legend('sup vs sqrt(Delta) - N=100', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
figure;
histogram(error_sup_M(:,2));
saveas(gcf, 'figures/EM_7_M_1000_N_1000.png');
legend('sup vs sqrt(Delta) - N=1000', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
figure;
histogram(error_alpha_sup_M(:,1));
saveas(gcf, 'figures/EM_7_M_1000_N_100_alpha.png');
legend('sup vs ({\Delta}t)^{alpha} - N=100', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
figure;
histogram(error_alpha_sup_M(:,2));
legend('sup vs ({\Delta}t)^{alpha} - N=1000', 'Orientation', 'vertical', 'Fontsize', 20, 'Location', 'SE');
saveas(gcf, 'figures/EM_7_M_1000_N_1000_alpha.png');
