Two correlated brownian motions

Question

Is it true (see here, footnote 2, p.22 / p.14, without proof) that we can obtain two discretized brownian motions $W_t^1, W_t^2$ with correlation $\rho$ by doing

$$d W_t^1 \sim \mathcal N(0,\sqrt{dt})$$

$$d W_t^2 = \rho\, d W_t^1 + \sqrt{1-\rho^2} dZ_t$$ with $dZ_t \sim \mathcal N(0,\sqrt{dt})$ ?

If this is true, we could easily simulate them in Python by doing:

import numpy as np
import matplotlib.pyplot as plt

n = 1000000; dt = 0.01; rho = 0.8
dW1 = np.random.normal(0, np.sqrt(dt), n)
dW2 = rho * dW1 + np.sqrt(1 - rho **2) * np.random.normal(0, np.sqrt(dt), n)
W1 = np.cumsum(dW1)
W2 = np.cumsum(dW2)
plt.plot(W1)
plt.plot(W2)
plt.show()

Is this correct?

Answer 1

First you need to correct the formula to: $$ W_t^2 = \rho W_t^1 + \sqrt{1-\rho^2} Z_t, $$ where $Z_t$ is a BM independent of $W_t^1$ If you calculate the variance and the covariance, then you see that it is true: $$ V[W_t^1] = t $$ and $$ V[W_t^2] = \rho^2 V[W_t^1] + (1-\rho^2) V[Z_t] = \rho^2 t + (1-\rho^2) t = t, $$ which is the desired variance.

For the covariance you get $$ Cov[W_t^1,W_t^2] = \rho Cov[W_t^1, W_t^1] + \sqrt{1-\rho^2} Cov[W_t^1, Z_t] $$ which gives (because $Cov[W_t^1, W_t^1] = V[W_t^1]$ and by indpenence of $W_t^1$ and $Z_t$: $$ Cov[W_t^1,W_t^2] = \rho t + 0, $$ and noting that $\sqrt{V[W_t^1]} = \sqrt{V[W_t^2]} = \sqrt{t}$ we get $$ Cor[W_t^1,W_t^2] = \frac{Cov[W_t^1,W_t^2]}{\sqrt{V[W_t^1]} \sqrt{V[W_t^2]} } = \frac{\rho t}{t} = \rho. $$

Answer 2

Here is the general approach you can follow to generate two correlated random variables. Let's suppose, X and Y are two random variable, such that: $$X \sim N(\mu_1, \sigma_1^2)$$ $$Y \sim N(\mu_2, \sigma_2^2)$$ and $$cor(X,Y)=\rho$$ Now consider: $y=bx + e_i$, where $x$ $(=\frac{X-\mu_1}{\sigma_1}$) and $y$ $(=\frac{Y-\mu_2}{\sigma_2}$) both follow standard normal distribution , such that $cor(x,y)=\rho. $ For standard normal variate, $b= \rho$. So we have:

$$y=\rho x + e_i$$

Now, here is the algorithm, you can follow:

1) Generate $n$ standard normal variate for $x$.

2) Since, $e_i \sim N(0, 1-\rho^2)$. So generates $n$ normal variate as $e_i$ from normal distribution with mean 0 and variance $1-\rho^2$.

3) Get $y=\rho x + e_i$.

4) Convert your standard normal numbers back to Normal (remember correlation is independent of change of origin and scale)

R code for simulating correlated GBM:

corGBM <- function(n, r, t=1/365, plot=TRUE) {
#n is number of samples 
#r is correlation
#t is tick step
x <- rnorm(n, mean=0, sd= 1)
se <- sqrt(1 - r^2) #standard deviation of error
e <- rnorm(n, mean=0, sd=se)
y <- r*x + e

X <- cumsum(x* sqrt(t))
Y <- cumsum(y* sqrt(t))
Max <- max(c(X,Y))
Min <- min(c(X,Y))

if(plot) {
plot(X, type="l", ylim=c(Min, Max))
lines(Y, col="blue")
}
return(cor(x,y))
}

#sample result
corGBM(10000,.85)
[1] 0.8523341