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
2 Answers
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.
The URL to the Matlab files in the paper won't work, use this one:
$\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.
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.
