Euler-Maruyama simulation of an SDE

Question

The Euler-Maruyama method for the following SDE

\begin{align} dX_{t} &= -\lambda X_{t}dt + \mu dW_{t}\\ X_{0}&=x>0 \end{align}

where $\lambda,\mu$ are given constants, is (according to Higham):

randn('state',100)
lambda = 1; mu = 0.7; Xzero = 0.5; 
T = 1; N = 2^8; dt = T/N;
dW = sqrt(dt)*randn(1,N); 
W = cumsum(dW); 
R = 4; Dt = R*dt; L = N/R; 
Xem = zeros(1,L);
Xtemp = Xzero;

for j = 1:L
    Winc = sum(dW(R*(j-1)+1:R*j));
    Xtemp = Xtemp - Dt * lambda * Xtemp + mu*Winc;
    Xem(j) = Xtemp;
end

plot([0:Dt:T],[Xzero,Xem],'r--*')

which works perfectly. My question is, how to simulate 100 such paths? To be more precise, I guess I have to simulate 100 Brownian motions and for each one to simulate the solution of the SDE.

dW = sqrt(dt)*randn(M,N);
W = cumsum(dW,2); 

I think that the above lines generate M different BM paths. Any specific ideas on the implementation?

Answer 1

This is a programming question rather than a math question, but regardless the required code is below. Read the Higham paper again as he also addresses this point. Also see code for the paper updated for MATLAB 2022 on GitHub.

rng(100,'v5normal');
lambda = 1; mu = 0.7; Xzero = 0.5; 
T = 1; N = 2^8; dt = T/N; M = 100;

dW = sqrt(dt)*randn(M,N); 
W = cumsum(dW,2); 
R = 1; Dt = R*dt; L = N/R; 
Xem = zeros(M,L);
Winc = zeros(M,L);
Xtemp = zeros(M,L);
for i = 1:M
    Xtemp(i,1) = Xzero;
    for j = 1:L
        Winc(i,j) = sum(dW(R*(j-1)+1:R*j),i);
        Xtemp(i,j+1) = Xtemp(i,j)-Dt*lambda*Xtemp(i,j)+mu*Winc(i,j);
        Xem(i,j) = Xtemp(i,j);
    end
end