Brownian Path Simulation using Matlab

%BPATH1  Brownian path simulation

randn('state',100)           % set the state of randn to be 100
T = 1; N = 500; dt = T/N;    % Length = 1, do 500 steps, each step measures dt=1/500
dW = zeros(1,N);             % preallocate arrays ...
W = zeros(1,N);              % for efficiency

dW(1) = sqrt(dt)*randn;      % first approximation outside the loop ...
W(1) = dW(1);                % since W(0) = 0 is not allowed
for j = 2:N
 dW(j) = sqrt(dt)*randn;   % general increment, dW ~ sqrt(dt)
 W(j) = W(j-1) + dW(j);    % then go a step of dW(j) from W(j-1) to W(j)
end

plot([0:dt:T],[0,W],'r-')    % plot W against t: plot from time 0 to time T, increment dt
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

Note: Matlab starts array from 1, not 0. So in order to include the starting point W(0) = 0 into the graph, we need to include it in the plotting command by doing plot([0:dt:T],[0,W])

Following is a more elegant and more efficient version of bpath1.m :
The "for" loop is replaced by a "vectorized" command
(KEY to writing efficient MATLAB)code (See DJHigham & NJHigha, MATLAB guide, SIAM, Philly 2000)


%BPATH2  Brownian path simulation: vectorized

randn('state',100)          % set the state of randn
T = 1; N = 500; dt = T/N;

dW = sqrt(dt)*randn(1,N);   % generate the array of increments
W = cumsum(dW);             % cumulative sum, W(j)= dW(1) + ... + dW(j)

plot([0:dt:T],[0,W],'r-')   % plot W against t, red line
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

A higher simulation of Brownian path:
Average from about 1000 different paths. Sketch the average path and 5 of the brownian paths.
Here we evaluate the function u(W(t)) = exp(t + W(t)/2) along 1000 discretized Brownian paths.
The average of u(W(t)) over these paths is plotted with a solid blue line.
Five individual paths are also plotted using a dashed red line.

%BPATH3  Function along a Brownian path

randn('state',100)                                % set the state of randn
T = 1; N = 500; dt = T/N; t = [dt:dt:1];          % t = row vector, N components, from dt to 1, increments dt.

M = 1000;                                         % M paths simultaneously
dW = sqrt(dt)*randn(M,N);                         % increments for M paths. 1 path 1 row. Each row N increments ...
                                                % i.e., dW is an M-by-N array st dW(i,j) gives the increment dWj

W = cumsum(dW,2);                                 % cumulative sum accross the second dimension (column)
                                                % i.e., each row of W is the cumsum of the corresp. row of dW

U = exp(repmat(t,[M 1]) + 0.5*W);                 % repmat(t,[M 1]) form a "matrix" M row 1 column, each "entry" is t.
                                                % With t a N-dim row vector, this put M t's stack together,
                                                % forming a M by N matrix.
                                                % then add it to the matrix .5*W
                                                  % then take exp of this matrix
Umean = mean(U);                                  % compute columwise averages
plot([0,t],[1,Umean],'b-'), hold on               % plot mean over M paths
                                                 % x-axis is 0 and the vector t
                                                 % y-axis is 1. At time 0 the path begin at y=1. WHY???
plot([0,t],[ones(5,1),U(1:5,:)],'r--'), hold off  % plot 5 individual paths
                                                 % discrete "time" from 0 to 1
                                                 % over 5 paths, all begin at 1
                                                 % paths are first 5 rows of U
xlabel('t','FontSize',16)
ylabel('U(t)','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
legend('mean of 1000 paths','5 individual paths',2)

averr = norm((Umean - exp(9*t/8)),'inf')          % sample error. Infinity norm of matrix
                                                 % largest row sum of A.
                                                 % w.o. 'inf' - 1-norm, largest column sum

To be explored....