%EM Euler-Maruyama method on linear SDE
%
% This can be modified to simulate the equation by the Milstein method by
% simply changing one line of the code where the update rule over a time
% step is defined
% SDE is dX = lambda*X dt + mu*X dW, X(0) = Xzero,
% where lambda = 2, mu = 1 and Xzero = 1.
%
clear *
clf
randn('state',100); % seed the random number generator so results are
% reproducible from one run to the next.
% Change this if one wants different
% results (choices of random variables)
% with each run.
lambda = 2; mu = 1; Xzero = 1; % problem parameters
T = 1; N = 2^9; % time interval, number of time steps
M=100; % number of realizations(experiments) to execute
dt = T/N; % time step
XTsum=0; % This is where we will store the running sum of the solution X(T)
% at the ending time, which will be used
% to estimate the mean value of X(T)
% from the simulations
for s=1:M % loop over realizations (experiments)
dW = sqrt(dt)*randn(1,N); % Brownian increments
W = cumsum(dW); % discretized Brownian path
Xtrue = Xzero*exp((lambda-0.5*mu^2)*(dt:dt:T)+mu*W); % exact solution for comparison
% Xem keeps track of the numerically simulated solution
% Xtemp is the numerically simulated value for the solution at the current
% time step
Xem = zeros(1,N); % preallocate for efficiency
Xtemp = Xzero; % start with the initial data
for j = 1:N
Winc = dW(j); % the Brownian increment over a time step
Xtemp = Xtemp + dt*lambda*Xtemp + mu*Xtemp*Winc; % Numerical update: Euler-Marayama
% To do Milstein simulation, simply replace this equation with the
% formula for the Milstein update rule
Xem(j) = Xtemp;
end
XTsum=XTsum+Xem(N); % Add the value X(T) from the current realization (experiment)
% to the running sum over all realizations
end
% These plot the results of the last realization (experiment)
plot(0:dt:T,[Xzero,Xtrue],'m-'), hold on
plot(0:dt:T,[Xzero,Xem],'r--*'), hold off
xlabel('t','FontSize',12)
ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
% This is to calcluate the error between the numerical and exact solution
% in the last realization.
emerr = abs(Xem(end)-Xtrue(end));
disp(sprintf('The numerical error is %f at the end time T=%f.',emerr,T))
figure(1)
XTmean=XTsum/M; % mean value of X(T) as estimated by numerical simulations
XTmeantrue=Xzero*exp(lambda*T); % the theoretically correct average value for the mean value of X(T).
disp(sprintf('The numerical estimate for the mean of X(T) is %f, while the exact answer is %f.',...
XTmean,XTmeantrue))