clear();

%% Data

rng(2022);

lambda = 1;
mu = 0.1;
X0 = 0.1;
T = 1;

phi = @(x) x;
Z = X0*exp(lambda*T);

L  = 10;

%% Decay of the variance

Nlarray = 2.^(0:L)';
hlarray = T./Nlarray;
Mlarray = (L*hlarray)/(hlarray(end)^2);

phiarray = cell(L + 1, 2);
meanarray = zeros(L + 1, 1);
Meanarray = zeros(L + 1, 1);

disp('Level 0');

Ml = Mlarray(1);  
hl = hlarray(1);  
Nl = Nlarray(1);

dW = sqrt(hl)*randn(Ml, Nl);
W = [zeros(Ml, 1) cumsum(dW, 2)];
DW = W(:, 2:end) - W(:, 1:end-1);

X = prod(1 + hl*lambda + mu*DW, 2).*(X0*ones(Ml, 1));
phiarray{1, 1} = phi(X);
phiarray{1, 2} = zeros(Ml,1);

meanarray(1) = sum(phiarray{1,1} - phiarray{1,2}, 1)/Ml;

X_sol = X0*exp((lambda - (mu^2)/2)*T + mu*W(end)); 
Meanarray(1) = abs(mean(phiarray{1, 1}-phi(X_sol)));

for l = 2:L+1
    
    disp(['Level ', num2str(l-1)]);
    
    Ml = Mlarray(l);
    hl = hlarray(l);
    Nl = Nlarray(l);
    
    dW = sqrt(hl)*randn(Ml, Nl);
    W  = [zeros(Ml, 1) cumsum(dW, 2)];
    
    DW = W(:, 2:end) - W(:, 1:end-1);
    X = X0*prod(1 + hl*lambda + mu*DW, 2);
    phiarray{l, 1} = phi(X);
    
    DW = W(:,3:2:end) - W(:,1:2:end-2);
    X = X0*prod(1 + 2*hl*lambda + mu*DW, 2);
    phiarray{l, 2} = phi(X);   
    
    meanarray(l) = sum(phiarray{l,1} - phiarray{l,2}, 1)/Ml;

    %X_sol = X0*exp((lambda - (mu^2)/2)*T + mu*W(end)); 
    Meanarray(l) = abs(mean(phiarray{l, 1}-phi(X_sol)));

end

barZ = sum(meanarray);
rMSE = realsqrt(mean((Z - barZ)^2, 1));

V = zeros(L + 1, 1);
Varphil = zeros(L + 1, 1);
Var = @(p) mean((p - mean(p,1)).^2, 1);

for l= 1:L+1
    d = phiarray{l,1} - phiarray{l,2};
    V(l) = Var(d);
    Varphil(l) = Var(phiarray{l,1});
end

%% Figure

figure();
ax = axes('nextplot','add','box','on');

semilogy((0:L)', V, ':x', 'Linewidth', 2, 'Markersize', 20);
semilogy((0:L)', Varphil, ':s', 'Linewidth', 2, 'Markersize', 20);
semilogy((0:L)', Meanarray, ':o', 'Linewidth', 2, 'Markersize', 20);

leg = legend({'$\mathrm{Var}(\phi_\ell-\phi_{\ell-1})$', '$\mathrm{Var}(\phi_\ell)$', '$\left|\mathrm{E}[\phi_\ell-\phi]\right|$'});
set(leg, 'box', 'off', 'interpreter', 'latex', 'fontsize', 20);
set(get(gca,'xlabel'), 'string', '$\ell$', 'interpreter', 'latex');
set(ax, 'yscale', 'log', 'Ylim', [10^-9 10^0]);
set(gca, 'FontSize', 20);
saveas(gcf, 'figures/exercise2_2.png');