Euler-Maruyama simulation of an SDE

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?

• Instead of working with a 1 by L vector Xem, try using a 100 by L matrix Xem, where each row corresponds to one path. Your code wont have to change very much. – David Nov 23 '15 at 21:51
• @David, that you for the answer. You are correct, this is the way I finally managed to do it. – peter5 Nov 26 '15 at 13:32
• Great to hear! :) – David Nov 26 '15 at 21:32

randn('state',100)
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);
for i=1:M
Xtemp(i,1) = Xzero;
for j = 1:L
Winc(i,j) = sum(i,dW(R*(j-1)+1:R*j));
Xtemp(i,j+1) = Xtemp(i,j) - Dt * lambda * Xtemp(i,j) + mu*Winc(i,j);
Xem(i,j) = Xtemp(i,j);
end
end