command randn(1,N) in matlab

Question

This program is in Matlab to simulate Brownian motion Generating GBM:

T = 1; N = 500; dt = T/N;
t = 0:dt:T;
dW = sqrt(dt)*randn(1,N);
mu = 0.1; sigma = 0.01;
plot(t,exp(mu*t + sigma*[0,cumsum(dW)]))

the command "randn(1,N)" what is it for is this generate array from 1 to N of of independent Gaussian random numbers, all of which are $N(0;1)$. because I know that the command The command will be randn(M,N), which generates an $M \times N $ matrix of independent Gaussian random numbers, all of which are $N(0;1)$. and please if anyone can explain to me what the line "dW = sqrt(dt)*randn(1,N);" means?

Thanks

Answer 1

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.

The function 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: https://www.maths.ed.ac.uk/~dhigham/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.

Answer 2

A normal random variable $\zeta$ distributed like $\mathcal{N}(\mu,\sigma^2)$ can be written in terms of a standard normal random variable $\chi$ very easily as:

$$\zeta = \mu + \sigma \chi.$$

Thus, your line of code is just generating a Gaussian random variable with a different standard deviation.