function [err,stat] = StrongEulerBS(S0, sigma, T, N, M)
    function y = Vdr(x)
        y = zeros(x);
    endfunction
    function y = Vdi(x)
        y = sigma*x;
    endfunction
    function y = BS(dB)
        B = sum(dB, 2);
        y = S0*exp(sigma*B-0.5*sigma^2*T);
    endfunction
    
    mu = 0;
    mu2 = 0;
    
    for i=1:M
        [Xbar, dB] = EulerTrajectory(S0, Vdr, Vdi, 1, N, T);
        X = BS(dB);
        mu = mu + abs(Xbar-X);
        mu2 = mu2 + (Xbar-X)^2;
    end
    err = mu/M;
    sd = sqrt(mu2/M-err^2);
    stat = 1.65*sd/sqrt(M);
endfunction