%stochastic intergral simulation
T=20; %final time

dt=0.00001;
t=[0:dt:T];
n=size(t);
n=n(2);

clf reset;
 
Bmt = cumsum(sqrt(dt)*randn(1,n));
plot(t,Bmt(1:n),'Color',[0 0 1]);
hold on;
   
inand=Bmt(1:n-1); %integrand, more choices below
   
%  inand=inand.^2;
%   inand=abs(inand); 
%   inand=sqrt(abs(inand)); 
%   inand=(abs(inand)).^(-0.4); %is in H^2 for power -alpha when alpha<1/2
%  inand=fliplr(inand);  %NOT an Ito integral!   
 
%plot(t(1:n-1), inand, 'Color',[1 0.5 0]);

intor=Bmt(2:n)-Bmt(1:n-1); %integrator
%intor=1/n; %for ordinary integral
pr=inand.*intor;
It=cumsum(pr); %integral

plot(t(1:n-1), It, 'Color',[1 0 0]);
legend('Bm', 'Int','Location', 'northeastoutside', 'AutoUpdate', 'off');
hold on;
% plot(t(1:n-1), It, 'Color',[1 0 0], 'Linewidth', 1);
% plot(t,0.5*(Bmt.^2-t),'Color',[0 1 0]); %check against explicit solition

plot(t(1:n-1), 0*t(1:n-1), 'Color',[0 0 0]);
title('Stochastic Integral');
hold off;

  