First, this code returns the analytical solution to the Stratonovich stochastic differential equation for geometric Brownian motion:
$$dY = \mu Y dt+\sigma Y \circ dW$$
The full solution to this is:
$$Y(t) = Y_0 \exp(\mu t+\sigma W_t)$$
Your code assumes $Y_0=1$. The solution is slightly different if you use the Itô formulation.
randn returns random variates (values) independently drawn from the normal distribution. So,
randn(1,N) returns a sequence of independent normally-distributed values (500 in this case).
N in this case is just the number of time-steps. These are used to produce the independent Wiener increments, $dW_t$. Thus,
dW = sqrt(dt)*randn(1,N) produces
N independent Wiener increments where the standard deviation is equal to the square root of the time span,
dt, between them $\dagger$. These correspond to the
N points in the vector
t. The code
[0,cumsum(dW)] integrates these increments to produce $W_t$, a Wiener process, which is also called standard Brownian motion
For further details on SDEs, Brownian motion, and simulating them with Matlab I recommend this excellent paper:
Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302
The URL to the Matlab files in the paper won't work, use this one: http://personal.strath.ac.uk/d.j.higham/algfiles.html
$\dagger$ Why is
dW scaled by
sqrt(dt) instead of just
dt like when numerically integrating ODEs? A discrete Wiener process by definition has a variance equal to the time step between the increments, in this case
dt. However, to scale output of
randn properly, we need to express this in terms of a standard deviation.