%EMSTRONG  Test strong convergence of Euler-Maruyama
%
% Solves    dX = lambda*X dt + mu*X dW,   X(0) = Xzero,
%         where lambda = 2, mu = 1 and Xzer0 = 1.
%
% Discretized Brownian path over [0,1] has dt = 2^(-9).
% E-M uses 5 different timesteps: 16dt, 8dt, 4dt, 2dt, dt.
% Examine strong convergence at T=1:  E | X_L - X(T) |.


clear *
clf

randn('state',100)
lambda = 2; mu = 1; Xzero = 1;    % problem parameters
T = 1; N = 2^9;                   % Simulation time T, number of time steps N
Nres=5;                           % Number of time step resolution levels
M = 1000;                         % number of paths sampled

dt = T/N;                         % Reference time step
Xerr = zeros(M,Nres);                % preallocate array
for s = 1:M,                      % sample over discrete Brownian paths
    dW = sqrt(dt)*randn(1,N);     % Brownian increments
    W = cumsum(dW);               % discrete Brownian path 
    Xtrue = Xzero*exp((lambda-0.5*mu^2)*T+mu*W(end)); % Correct value of solution at end time
    for p = 1:Nres                              % Loop over time step resolution levels
        R = 2^(p-1); Dt = R*dt; L = N/R;     % L Euler steps of size Dt = R*dt
        Xtemp = Xzero;                       % Initial condition
        for j = 1:L
             Winc = sum(dW(R*(j-1)+1:R*j));  % Brownian increment over coarse time step Dt 
             Xtemp = Xtemp + Dt*lambda*Xtemp + mu*Xtemp*Winc; % Euler-Marayama update scheme
        end
        Xerr(s,p) = abs(Xtemp - Xtrue);      % store the error at end time
    end
end

Dtvals = dt*(2.^(0:Nres-1));                       % Time steps in the simulations
subplot(221)                                  % top LH picture
loglog(Dtvals,mean(Xerr),'b*-'), hold on
loglog(Dtvals,(Dtvals.^(.5)),'r--'), hold off % reference slope of 1/2 for theoretical scaling
axis([1e-3 1e-1 1e-4 1])                      % of numerical discretization error
xlabel('\Delta t'), ylabel('Sample average of | X(T) - X_L |')
title('Euler-Marayama Strong Error','FontSize',10)

%%%% Least squares fit of error = C * Dt^q %%%%
A = [ones(Nres,1), log(Dtvals)']; rhs = log(mean(Xerr)');
sol = A\rhs; q = sol(2);
resid = norm(A*sol - rhs);
figure(1)
disp(sprintf('The least squares fit for the slope of the strong E-M error in the log-log plot is %f, with residual %f.',q,resid))