Monte-carlo Simulation Notes-1

\(d\)-dimensional Normal random variables

• Is characterized by \(d\)-vector \(\mu\), \(d\times d\) covariance matrix \(\Sigma\). \(\Sigma\) should be positive definite and symmetric. The density of the \(d\)-dimensional normal random variable is given by:

\[\phi_{\mu, \Sigma} = \frac{1}{(2\pi)^{d/2}\vert \Sigma \vert^{1/2}} \exp{\left(-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)\right)}.\]

• If \(Z \sim \mathcal(0, I)\), then \(X = \mu + A Z \sim \mathcal(\mu, AA^{T})\). It remains to find \(A\) such that \(AA^{T} = \Sigma\). We make use of Cholesky factorization for this.

Cholesky factorization

• For a \(2\times 2\), \(\Sigma\), we find a lower triangular matrix with the first row begin \(\sigma_1, 0\) and the second row \(\rho \sigma_2\), \(\sqrt{1-\rho^2}\sigma_2\).
• In case of bivariate normal, we use the following equations:

\[\begin{align*} X_1 &= \mu_1 + \sigma_1 Z_1, \\ X_2 &= \mu_2 + \rho \sigma_2 Z_1 + \sqrt{1-\rho^2}\sigma_2 Z_2 \\ \end{align*}\]

• For a general case, with covariance matrix \(\Sigma\), we get the following relations

\[\begin{align*} A_{ij} & = \frac{1}{A_jj}\left(\Sigma_{ij} - \sum_{k=1}^{j-1} A_{ik}A_{jk}\right), j<i \\ A_{ii} &= \sqrt{\Sigma_{ii} - \sum_{k=1}^{j-1} A^2_{ik}} \end{align*}\]

Brownian motion

• Wiener process A Wiener process (or Brownian motion; \(W_t\) or \(W\)) is a time continuous process with the following properties:

  1. \(W_0 = 0\).
  2. \(W_t \sim \mathcal N(0, t)\) for all \(t \ge 0\). That is, for each \(t\), the random variable \(W_t\) is distributed normally with mean \(E(W_t) = 0\) and variance \(\text{Var}(W_t) = E(W_t^2) = t\).
  3. All increments \(\Delta W_t := W_{t+\Delta t} + W_t\) on non-overlapping time intervals are independent.
  4. \(W_t\) depends continuously on \(t\).

Algorithm for simulation of standard Wiener process

• Start: \(t_0 = 0, W_0 = 0; \Delta t\)
• loop $j = 1,2, $:
• \(t_j = t_{j-1} + \Delta t\)
• draw \(Z \sim N(0,1)\).

• \(W_j = W_{j-1} + Z \sqrt{\Delta t}\)

• Almost all realization of Wiener processes are nowhere differentiable.

• For non-standard Brownian motion (\(\text{BM} (\mu,\sigma^2)\), the values are generated by the following equation:

\[ X(t_{i+1}) = X(t_i) + \mu(t_{i+1}-t_{i}) + \sigma \sqrt{t_{i+1} - t_{i}} Z_{i+1}, \quad i = 0, 1, \cdots, n-1.\]

• With time-dependent coefficients, the recursion becomes

\[X(t_{i+1}) = X(t_i) + \int_{t_i}^{t_{i+1}} \mu(s) \, \text ds + \sqrt{\int_{t_i}^{t_{i+1}} \sigma^2(u) \text du}\cdot Z_{i+1}, \quad i = 0, \cdots, n-1.\]

• The Euler approximation of the above formula is given by

\[ X(t_{i+1}) = X(t_i) + \mu(t_i)(t_{i+1}-t_i) + \sigma(t_i) \sqrt{t_{i+1} - t_{i}} Z_{i+1}, i = 0, \cdots, n-1\]

Geometric Brownian motion

• Geometric Brownian motion is valuable in modelling of stock prices. An ordinary Brownian motion can take negative values, whereas the Geometric Brownian motion can take only positive values which is why it is useful in modelling stock markets.
• If \(S \sim \text{GB}(\mu, \sigma^2)\) and if \(S\) has an initial value of \(S(0)\), then \(S(t) = S(0) \exp((\mu - \frac{1}{2} \sigma^2) t + \sigma\cdot W(t))\). or more generally, if \(u < t\), then

\[ S(t) = S(u) \exp((\mu - \frac{1}{2} \sigma^2)\cdot (t-u) + \sigma(W(t) - W(u)).\]

• A recursion for simulating values is,

\[ S(t_{i+1}) = S(t_i) \exp((\mu - \frac{1}{2} \sigma^2) \cdot (t_{i+1} - t_{i}) + \sigma (\sqrt{t_{i+1} - t_{i}}) \cdot Z_i)\]

• We also try to simulate jump at various time period. The algorithm is listed below
• Generate \(Z \sim \mathcal N(0, 1)\)
• Generate \(N \sim \text{Poisson}(\lambda(t_{i+1} - t_{i})\)
• If \(N \neq 0\), then
  • Generate \(log(Y_1), \cdots, log(T_N)\) from their common distribution
  • Set $ M = (Y_i)$
• If \(N \neq 0\), then
  • Set \(M = 0\).
• Set \(X(t_{i+1}) = X(t_i) + \left(\mu - \frac{1}{2} \sigma^2\right)(t_{i+1} - t_{i}) + \sigma \sqrt{t_{i+1} - t_{i}}\cdot Z + M.\)

Code in R

{% highlight R %}

sigma = 0.3 mu = 0.06 N = 1000 flag = 0 lambda = c(0.01, 0.05, 0.1, 0.5, 1) xs = c(5) Xs = matrix(5, nrow = 5, ncol = N) zs = rnorm(N) for (i in 2:N) { xs[i] = xs[i-1] + (mu - 0.5sigmasigma)/N + sigmasqrt(1/N)zs[i] } for (j in 1:5) { for (i in 2:N) { n = rpois(1, lambda[j]1/N) M = 0 if(n != 0) { flag <<- flag+1 ys = 0.1rnorm(n, mean = 0, sd = 1)+1 M = 0 for (k in 1:n) M = M + log(ys[k]) } Xs[j, i] = Xs[j,i-1] + (mu-0.5sigmasigma)/N + sigmasqrt(1/N)zs[i] + M } } yM = max(max(xs), max(Xs)) ym = min(min(xs), min(Xs)) plot(xs, cex = 0.05, type=“l”, ylim=c(ym, yM)) lines(1:N, Xs[1,1:N], col=“red”) lines(1:N, Xs[2,1:N], col=“blue”) lines(1:N, Xs[3,1:N], col=“green”) lines(1:N, Xs[4,1:N], col=“violet”) lines(1:N, Xs[5,1:N], col=“orange”)

{% endhighlight %}