Simulate 1,000 geometric brownian motions in MATLAB

Question

I currently have code to simulate a geometric Brown motion, courtesy of http://www-math.bgsu.edu/~zirbel/sde/matlab/index.html.

However, I would like to generate 1,000 simulations and to be to display them in a graph.

The codes I have at the moment to generate a single simulation are as follows:

% geometric_brownian(N,r,alpha,T) simulates a geometric Brownian motion 
% on [0,T] using N normally distributed steps and parameters r and alpha

function [X] = geometric_brownian(N,r,alpha,T)

t = (0:1:N)'/N;                   % t is the column vector [0 1/N 2/N ... 1]

W = [0; cumsum(randn(N,1))]/sqrt(N); % S is running sum of N(0,1/N) variables

t = t*T;
W = W*sqrt(T);

Y = (r-(alpha^2)/2)*t + alpha * W;

X = exp(Y);

plot(t,X);          % plot the path
hold on
plot(t,exp(r*t),':');
axis([0 T 0 max(1,exp((r-(alpha^2)/2)*T+2*alpha))])
title([int2str(N) '-step geometric Brownian motion and its mean'])
xlabel(['r = ' num2str(r) ' and alpha = ' num2str(alpha)])
hold off

Solution

That code cannot be used directly to simulate 1,000 paths/simulations. Unfortunately, it has not been vectorized. The easiest way to do what you want is to use a for loop:

N = 1e3;
r = 1;
alpha = 0.1;
T = 1;
npaths = 1e3;          % Number of simulations

rng(0);                % Always set a seed
X = zeros(N+1,npaths); % Preallocate memory
for i = 1:n
    X(:,i) = geometric_brownian(N,r,alpha,T);
    hold on
end
t = T*(0:1:N).'/N;
plot(t,exp(r*t),'r--');

This is rather slow and inefficient. You will need to modify the function a lot to vectorize it. One thing that would improve performance is if you at least removed the plotting code from inside the function and ran that separately after the loop.

Another alternative might be to use the sde_gbm function in my SDETools toolbox, which is fully-vectorized and much faster:

N = 1e3;
r = 1;
alpha = 0.1;
T = 1;
npaths = 1e3;        % Number of simulations

t = T*(0:1:N)/N;     % Time vector
y0 = ones(npaths,1); % Vector of initial conditions, must match number of paths
opts = sdeset('RandSeed',0,'SDEType','Ito'); % Set seed
y = sde_gbm(r,alpha,t,y0,opts);

figure;
plot(t,y,'b',t,y0*exp(r*t),'r--');
xlabel('t');
ylabel('y(t)');
title(['Geometric Brownian motion and it's mean: ' int2str(npaths) ...
       ' paths, r = ' num2str(r) ', \alpha = ' num2str(alpha)]);

In either case, one obtains a plot that looks something like this

               

Comments

[horchler]

I think the OP is asking how to generate 1,000 independent simulations (or paths in Brownian motion parlance) for 0 to T, not 1,000 time-steps from a single simulation.