sde_demo.m

SIMULATION AND ESTIMATION OF STOCHASTIC DIFFERENTIAL EQUATIONS WITH MATLAB

figure
lbxItoMilstein = perctile(xItoMilstein',2.5); % Lower Bound of the 95% confidence area obtained by taking, at each time, the 2.5% percentile of the simulated trajectories
waitbar(10/15);
ubxItoMilstein = perctile(xItoMilstein',97.5); % Upper Bound of the 95% confidence area obtained by taking, at each time, the 97.5% percentile of the simulated trajectories
waitbar(11/15);
q1xItoMilstein = perctile(xItoMilstein',25); % first quartile of the simulated trajectories
waitbar(12/15);
q3xItoMilstein = perctile(xItoMilstein',75); % third quartile of the simulated trajectories
waitbar(13/15);
meanxItoMilstein = mean(xItoMilstein,2);
waitbar(14/15);
%                 the empirical 95% CI,                                     first and third quartile                         the process empirical mean
plot(OWNTIME,lbxItoMilstein,'k--',OWNTIME,ubxItoMilstein,'k--',OWNTIME, q1xItoMilstein,'k:',OWNTIME,q3xItoMilstein,'k:',   OWNTIME,meanxItoMilstein,'g-') 
waitbar(15/15);
xlabel('t','Fontsize',13,'Rotation',0);
ylabel('X_t','Fontsize',13,'Rotation',0);
titlestring = sprintf('Ito SDE: mean, 95 percent CI, q1-q3 quartiles of the Milstein approximation over %d trajectories',NUMSIM);
title(titlestring);
close(handle); % close the handle defined with the waitbar command

% Plots the histogram of the trajectories at the end-time (Euler-Maruyama approximation)
figure
hist(xItoEuler(end,:),20);
xlabel('X_T','Fontsize',13,'Rotation',0);
titlestring = sprintf('Ito SDE: Histogram of X_t at end-time T=%d with Euler-Maruyama approximation',OWNTIME(end));
title(titlestring,'Fontsize',12);
% Plots the histogram of the trajectories at the end-time (Milstein approximation)
figure
hist(xItoMilstein(end,:),20);
xlabel('X_T','Fontsize',13,'Rotation',0);
titlestring = sprintf('Ito SDE: Histogram of X_t at end-time T=%d with Milstein approximation',OWNTIME(end));
title(titlestring,'Fontsize',12);


%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
%                      STRATONOVICH SDE COMPUTATIONS
%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
SDETYPE = 'Strat';  % specifies the SDE type (Ito or Strat (Stratonovich))

xStratMilstein = zeros(N,NUMSIM*NUMDEPVARS);  
errorT_StratMilstein = zeros(NUMSIM,NUMDEPVARS);
lbxStratMilstein = zeros(NUMSIM,NUMDEPVARS);
ubxStratMilstein = zeros(NUMSIM,NUMDEPVARS);
q1xStratMilstein = zeros(NUMSIM,NUMDEPVARS);
q3xStratMilstein = zeros(NUMSIM,NUMDEPVARS);
meanxStratMilstein = zeros(NUMSIM,NUMDEPVARS);

 
%::: computing trajectories with Milstein scheme::::::::::::
xStratMilstein = SDE_milstein_demo(bigtheta);            % the Milstein solution of the Stratonovich SDE, COMPUTED ON THE SAME BROWNIAN PATHS OF THE TRUE SOLUTION
errorT_StratMilstein = abs(xStratMilstein(end,:)-xTrue(end,:));       % the absolute error with Milstein (Stratonovich) solution at time T 


% plot graphs
handle = waitbar(0,'Stratonovich SDE: computing statistics and plotting simulations...');
figure
plot(OWNTIME,xStratMilstein,'k--'), hold on
waitbar(1/8);
plot(OWNTIME,xTrue,'m-'), hold off
waitbar(2/8);
xlabel('t','Fontsize',13,'Rotation',0);
ylabel('X_t','Fontsize',13,'Rotation',0);
titlestring = sprintf('Stratonovich SDE: Milstein vs analytic solution over %d trajectories',NUMSIM);
title(titlestring);

% Now compute mean, 95% confidence intervals and q1-q3 quartiles for the Milstein approximation
figure
lbxStratMilstein = perctile(xStratMilstein',2.5); % Lower Bound of the 95% confidence area obtained by taking, at each time, the 2.5% percentile of the simulated trajectories
waitbar(3/8);
ubxStratMilstein = perctile(xStratMilstein',97.5); % Upper Bound of the 95% confidence area obtained by taking, at each time, the 97.5% percentile of the simulated trajectories
waitbar(4/8);
q1xStratMilstein = perctile(xStratMilstein',25); % the first quartile of the simulated trajectories
waitbar(5/8);
q3xStratMilstein = perctile(xStratMilstein',75); % the third quartile of the simulated trajectories
waitbar(6/8);
meanxStratMilstein = mean(xStratMilstein,2);
waitbar(7/8);
%                 the empirical 95% CI,                                             first and third quartile                      the process empirical mean
plot(OWNTIME,lbxStratMilstein,'k--',OWNTIME,ubxStratMilstein,'k--',OWNTIME,q1xStratMilstein,'k:',OWNTIME,q3xStratMilstein,'k:',   OWNTIME,meanxStratMilstein,'g-') 
waitbar(8/8);
xlabel('t','Fontsize',13,'Rotation',0);
ylabel('X_t','Fontsize',13,'Rotation',0);
titlestring = sprintf('Stratonovich SDE: mean, 95 percent CI, q1-q3 quartiles of the Milstein approximation over %d trajectories',NUMSIM);
title(titlestring);
close(handle); % close the handle defined with the waitbar command

% Plots the histogram of the trajectories at the end-time (Milstein approximation)
figure
hist(xStratMilstein(end,:),20);
xlabel('X_T','Fontsize',13,'Rotation',0);
titlestring = sprintf('Stratonovich SDE: Histogram of X_t at end-time T=%d with Milstein approximation',OWNTIME(end));
title(titlestring,'Fontsize',12);


%-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
%                             STATISTICS
%-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

fprintf('\n\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n                                                        ITO SDE WITH EULER-MARUYAMA APPROXIMATION: STATISTICS');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');

avgerrorT_ItoEuler = sum(errorT_ItoEuler)/NUMSIM;             % the average absolute error at time T along the NUMSIM simulations with Euler-Maruyama scheme (Ito SDE)
fprintf('\n\nIto SDE: Average (absolute) error at time %d with Euler-Maruyama method: %d',OWNTIME(end),avgerrorT_ItoEuler);
SDE_stats(bigtheta,xItoEuler,PROBLEM,OWNTIME,NUMDEPVARS,NUMSIM,SDETYPE,'EM',0);

fprintf('\n\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n                                                        ITO SDE WITH MILSTEIN APPROXIMATION: STATISTICS');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');

avgerrorT_ItoMilstein = sum(errorT_ItoMilstein)/NUMSIM;       % the average absolute error at time T along the NUMSIM simulations with Milstein scheme (Ito SDE)
fprintf('\n\nIto SDE: Average (absolute) error at time %d with Milstein method: %d',OWNTIME(end),avgerrorT_ItoMilstein);
SDE_stats(bigtheta,xItoMilstein,PROBLEM,OWNTIME,NUMDEPVARS,NUMSIM,SDETYPE,'MIL',0);

fprintf('\n\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n                                                        STRATONOVICH SDE WITH MILSTEIN APPROXIMATION: STATISTICS');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n----------------------------------------------------------------------------------------------------------------------------------------------------------------');

avgerrorT_StratMilstein = sum(errorT_StratMilstein)/NUMSIM;       % the average absolute error at time T along the NUMSIM simulations with Milstein scheme (Stratonovich SDE)
fprintf('\n\nStratonovich SDE: Average (absolute) error at time %d with Milstein method: %d\n\n',OWNTIME(end),avgerrorT_StratMilstein);
SDE_stats(bigtheta,xStratMilstein,PROBLEM,OWNTIME,NUMDEPVARS,NUMSIM,SDETYPE,'MIL',0);