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