%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))