%EMWEAK  Test weak convergence of Euler-Maruyama
%
% Solves    dX = lambda*X dt + mu*X dW,   X(0) = Xzero,
%        where lambda = 2, mu = 1 and Xzer0 = 1.
%
% E-M uses different timesteps: 2^(p-10),  p = 1,2,3,4,Nres
% Examine weak convergence at T=1:   | E (X_L) - E (X(T)) |.
%
% Different paths are used for each E-M timestep.
% Code is vectorized over paths.
%
% Uncommenting the line indicated below gives the weak E-M method.

clear *

randn('state',100);                      
lambda = 2; mu = 0.1;                    % problem parameters
Xzero = 1;                               % initial condition
T = 1;                                   % simulation time
N= 2^10;                                 % number of time steps at finest resolution
M = 50000;                               % number of paths sampled
Nres=5;                                  % number of resolution levels

dt=T/N;                                  % time step at finest resolution
Xem = zeros(Nres,1);                        % preallocate arrays
for p = 1:Nres                              % take various Euler timesteps
       Dt = 2^p*dt; L = T/Dt;          % L Euler steps of size Dt 
       Xtemp = Xzero*ones(M,1);            
       for j = 1:L
           Winc = sqrt(Dt)*randn(M,1);   % Brownian increment.  No need to use same realization
                                         % of randomness in different
                                         % simulations because only weak
                                         % error is examined
            %Winc = sqrt(Dt)*sign(rand(M,1)-0.5); %% use for weak E-M %%
           Xtemp = Xtemp + Dt*lambda*Xtemp + mu*Xtemp.*Winc; % Euler-Marayama update
       end
       Xem(p) = mean(Xtemp); % mean of X(T) sampled from the simulations at resolution level p
end
Xerr = abs(Xem - exp(lambda*T)); % error from exact solution:  X(T) should have mean exp(lambda*T) 

Dtvals = 2.^(1:Nres)*dt;                 % Time step resolution levels
subplot(222)                             % top RH picture
loglog(Dtvals,Xerr,'b*-'), hold on
loglog(Dtvals,Dtvals,'r--'), hold off    % reference slope of 1
axis([1e-3 1e-1 1e-4 1])
xlabel('\Delta t'), ylabel('| E(X(T)) - Sample average of X_L |')
title('Euler-Marayama Weak Error','FontSize',10)

%%%% Least squares fit of error = C * dt^q %%%%
A = [ones(Nres,1), log(Dtvals)']; rhs = log(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 weak E-M error in the log-log plot is %f, with residual %f.',q,resid))