Stochastic and Random Numerics

Pseudo-random number generators

One of the simplest types of pseudo-random number generators are the linear congruential generators (LCGs). One loop of an LCG consists of a multiplication (by an integer $a$ ), a shift by an integer $b$ , taking the remainder by an integer $m$ (the modulus), and then dividing the result by $m$ :

$x_{i+1} =(a x_i + b) \ \text{mod} \ m$

$u_{i+1} = x_{i+1}/m$

One proposed choice with relatively small numbers is the Park and Miller "minimal standard" LCG which uses $m = 2^{31} - 1$ , $a = 7^5$ , and $b=0$ . The modulus $m$ is a Mersenne prime. This LCG is quite poor by current standards, but it might be sufficient for some purposes (e.g. in a game for which speed is more important than true randomness).

Stochastic ODEs: The Euler-Maruyama Method

In 1905 Albert Einstein famously published four groundbreaking papers, on (1) the photoelectric effect, (2) special relativity, (3) mass-energy equivalence (the " $E = mc^2$ " paper), and (4) Brownian motion. The last of these is perhaps less well-known than the others, but its analysis founded an extremely important branch of statistical mechanics and probabilistic modeling.

Brownian motion can be thought of as a mathematical limit of a random walk. In particular, consider a 1-dimensional process in which the state (a number) changes by a fixed positive or negative increment with equal probability. The figure below shows a 10 realizations of such a process, for 100 steps, in each case starting at 0 and with increments 1 and -1.

Random walk examples — 10 random walks of length 100, with increments 1 and -1.

If we decrease the time taken for each step by a factor of $k$ , the variance of the total of all steps taken in a fixed time period will increase by a factor of $k$ . So if we want the standard deviation of the process to stay constant as we decrease this interval, the standard deviation of each step needs to be decreased by $\sqrt{k}$ . The result is the Brownian motion stochastic process (sometimes called the Weiner process in honor of Norbert Wiener), which starting from $x_0=0$ has the probability distribution:

$p_B(x) = \frac{1}{\sqrt{2 \pi t}} e^{-x^2/(2t)}$

Stochastic differential equations are usually written in terms of increments rather than derivatives, since they are not usually differentiable functions (they are best defined as integral equations). For a SDE of the form

$dX = f(x,t)dt + g(x,t)dB$

we can simulate an initial value problem's solution using the Euler-Marayama method, which generalizes the usual Euler's method. If $x(t_0) = x_0$ , then for a timestep $k = \Delta t$ we approximate using

$x_{i+1} = x_i + f(x_i, t_i) k + g(x_i, t_i) \sqrt{k} w_i$

where $w_i$ is a random sample drawn from a normal distribution with mean $0$ and variance $1$ .

The key point is the scaling of the stochastic component by $\sqrt{k}$ , rather than $k$ .

Example: Stock price model

$dS = r S dt + \sigma S dB_t$

"solution" is $S = S_0 e^{(r-\sigma^2/2)t + \sigma B_t}$ .

The Brownian Bridge

$dy = \frac{y_1 - y}{t_1 - t}dt + dB_t, \ \ \ y(t_0) = y_0$

Exercises

Find the stochastic 'differential' equation for the function $y(t)=2 B_t^2 -2 t + 3$ (including its initial condition), where $B_t$ is the standard Brownian motion process.