function temple5044_euler_maruyama
%TEMPLE5044_EULER_MARUYAMA
%   Conducts a numerical error convergence study of the
%   Euler-Maruyama method applied to the stochastic
%   differential equation [geometric Brownian motion]
%     dX = mu*X dt + sigma*X dW.
%
%   (C) 2019/04/15 by Benjamin Seibold
%       http://www.math.temple.edu/~seibold/

mu = 1; % drift
sigma = 1; % diffusion
X0 = 1; % initial value
T = 1; % final time
n_mc = 1e4; % number of Monte-Carlo samples
nt_vec = 25*2.^(0:7); % numbers of time steps

clf
subfigs = length(nt_vec); % number of plots
n_cols = ceil(sqrt(subfigs+1)); % number of subplot columns
n_rows = ceil((subfigs+1)/n_cols); % number of subplot rows
strong_error = nt_vec*0; weak_error = nt_vec*0;
for i = 1:subfigs % loop over choices of time step
    nt = nt_vec(i); % number of time steps
    k = T/nt; % time step size
    t = linspace(0,T,nt+1); % time vector
    R = randn(n_mc,nt);
    W = zeros(n_mc,nt+1);
    X_EM = X0*ones(n_mc,nt+1);
    for n = 1:nt % numerical integration
        W(:,n+1) = W(:,n) + sqrt(k)*R(:,n); % Brownian motions
        dW = W(:,n+1)-W(:,n); % random increment
        X_EM(:,n+1) = X_EM(:,n) + ... % Euler-Maruyama update
            mu*X_EM(:,n)*k + sigma*X_EM(:,n).*dW;
    end
    X = X0*exp(sigma*W+(mu-sigma^2/2)*(ones(n_mc,1)*t)); % true sol.
    strong_error(i) = mean(abs(X_EM(:,end)-X(:,end)));
    weak_error(i) = max(abs(mean(X_EM)-mean(X)));

    % Plot computational results
    subplot(n_rows,n_cols,i)
    ind_plot = round(linspace(1,n_mc,200)); % plot a few paths
    hl = plot(t,X_EM(ind_plot,:)-X(ind_plot,:),'c.-',...
        [0,T],[0,0],'k-'); hl = hl(1); hold on
    hl = [hl;plot(t,mean(abs(X_EM-X)),'b*-',...
        t,mean(X_EM)-mean(X),'r.-')];
    legend(hl,'Y-X (approx.-truth)','E[|X(T)-Y(T)|]',...
        'E[X(T)]-E[Y(T)]','Location','SouthWest')
    xlabel('t'), ylabel('Y-X')
    axis([0,T,[-1,1]*strong_error(i)*5])
    title(sprintf('Euler-Maruyama with k=1/%0.0f',1/k))
    text(T*.02,strong_error(i)*4.3,sprintf(...
        'Error: strong=%0.2e; weak=%0.2e',...
        strong_error(i),weak_error(i)))
end
% Plot convergence curves
k_range = [T/max(nt_vec),T/min(nt_vec)];
subplot(n_rows,n_cols,subfigs+1)
loglog(k_range,k_range.^0.5,'m--',k_range,k_range,'c--',...
    T./nt_vec,strong_error,'r.-',T./nt_vec,weak_error,'b.-')
xlim(k_range), grid on
legend('slope 1/2','slope 1','strong error','weak error',...
    'Location','SouthEast')
xlabel('k'), ylabel('error')
title('Error convergence')
