% generates  Wiener Process paths:
% place all of the paths with a starting point of 0
% for no particularly good reason:
randn('state',0)
dt = 0.01;  % Delta t: You should try varying this
N = 10000; % Number of time steps (number of increments)
% t_i = t_{i-1} + dt;
Tf = N*dt; % final time
tt = linspace(0,Tf,1000); % time array
% W(t_i) = \sum_j=1^N (dW(j))
dW = sqrt(dt)*randn(N,1); % a vector of normal variantes, with var = dt;
W = zeros(N,1); % collects W(t)
for ii = 1:N-1
    W(ii+1,1) = sum(dW(1:ii,1)');
end
plot(tt,W','linewidth',2)
xlabel('time')
ylabel('W')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Comparison of analytical and computational
% estimates of the Ito calculation of 
% \int_0^t W(s) dW_s
intwdw = zeros(N,1);
for ii = 1:N-1
    intwdw(ii+1,1) = intwdw(ii,1) + dW(ii,1)^2;
end
figure(2)
plot(tt,intwdw,'o',tt,tt,'r')
xlabel('time')
ylabel('\int WdW')
legend('sum','estimate')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finally, we approximate 
% \int_0^t exp(-0.1 s) dW(s)

exps = exp(-0.1*tt');
intexpdw = zeros(N,1);
for ii = 1:N-1
    intexpdw(ii+1,1)  = intexpdw(ii,1) + exps(ii,1)*dW(ii,1);
end
figure(3)
plot(tt,intexpdw)
xlabel('time')
ylabel('\int exp(-t/10) dW')