Applying Quasi-Monte Carlo and Variance Reduction to Option Pricing

# Speeding up computations: when to use Quasi and standard Monte-Carlo in pricing


# Speeding up computations: when to use Quasi and standard Monte-Carlo in pricing












I am familiar with the theory of Monte-Carlo techniques in the numerical integration, and recently I have started my experiments with these methods applied to derivatives pricing. I am using Glasserman's book as my reference.

I've started by computing simply ATM vanilla call price. My setting is very basic: I have 100 time steps, I use Black-Scholes model (so that I can check the validity of the results) and Euler-Maryuama, Runge Kutta and Milstein numerical integration schemes. It appears that I need about $10^6$ samples to get stable results, and in Python that takes about 6 minutes. I have a rather fast machine, and I find this speed to be really low. Initially I thought the reason is in random sampler. However from what I've checked, sampling normal random variables does take some time, but $80\%$ of the computational time comes from the arithmetic operations (addition and multiplication) - even in case of basic Euler-Maryuama scheme. I do understand that in case of Geometric Brownian Motion there are much smarter simulation techniques available, but I'd like to keep it generic since I ultimately would like to work with more general diffusions.

Perhaps that's just an issue with Python, but that issue made me thought of speeding up my Monte-Carlo simulations. In Glasserman's book there are two types of speeding procedures:

- Variance reduction, including control variates, antithetic sampling and importance sampling. I've used the latter before, and my best bet is on that.

- Quasi Monte Carlo (low discrepancy sequences). I've never used it before, and it seems now to me that it can only be applied to compute integrals in the usual form, say when the densities are given explicitly just the analytic formula for the integral is unknown.

So far my goal is to be able to price exotic options with Monte-Carlo. Say, the model is 1-dimensional, and exotic payoff is extremely path-dependent. Perhaps later I'd need to price also an American version of such option. My questions are as follows:

- Which of the variance reduction methods shall I pay special attention to?

- Is QMC at all possible to apply in such problem, or do I really need to have an integral expression (over a finite-dimensional domain) for the option price?

To clarify my question regarding QMC. I see three ways (Q)MC can be applied in pricing.

- to sample increments of Brownian motion/jump components over the real line

- to sample whole paths from the path space

- if the value of the option have been found to have shape $$ \tag{1} \int_D p(S_1,\dots,S_n)f(S_1,\dots,S_n)\;\mathrm dS_1\dots \mathrm dS_n $$ where $p$ is a payoff function, $S_1,\dots, S_n$ are payoff variable (multiple stocks or multiple time measurements of stock or whatever else), then (Q)MC can also be used to sample from $D$ to compute the integral in $(1)$ numerically.

As far as I understood from Glasserman's book, he only considers QMC as a good solution for the very last method, where he provides some evidence that it beats usual random MC. In contrast, for the first two methods he talks about usual MC and variance reduction techniques. So my question regarding QMC is: can it be successfully applied in the first two methods, or it is not well-designed for them and I shall not expect much benefits from using QMC in the first two methods?

## Answer by Olaf (score 8, accepted)

https://quant.stackexchange.com/a/18069

By definition the fair value of an option is given by an expectation value of the payoff, $\mathbf{E}\left[\textrm{payoff}(\textit{paths})\right]$. The probability distribution of the `paths` is the risk neutral measure. This is just an integral expression of the form you wrote. This applies to all option prices. Many options are, of course, special in the sense that they only depend on the terminal value of each path.

Quasi Monte Carlo can therefore always be applied in these situations, and the sampling of paths is essentially the same thing as computing an integral (i.e. the expectation value of the option payoff). Glasserman's comment most likely refers to the fact that if you simply use QMC to generate paths, then it is possible that you are not receiving the full benefit of the low discrepancy of the random numbers generated by QMC. Put differently: the samples have a low discrepancy, but the generated paths might not. So the variance reduction might not be that big. This is all relative, in the sense that other variance reduction methods might produce better results (i.e. lower variance in the price estimate).

Now, your actual problem doesn't have anything to do with QMC. It's Python which is holding you back. Schemes such as Euler-Maryuama, Runge Kutta and Milstein are all schemes which cannot be vectorized. Therefore I suspect that you might have some code snippets of the form:

```
import numpy as np
... some definitions ...
for i in range(1, N_timesteps):
    for j in range(N_paths):
        S[i, j] = S[i - 1, j] * np.exp((r-0.5*sigma**2)*dt
                 + sigma*np.sqrt(dt) * np.random.randn())
```

The problem is that this is usually very slow in Python. Things you can do to speed this up:

- Sample all your random numbers outside the loop.

- Pre-compute as many factors outside the loops as possible.

- Vectorize inner loops, if possible.

- Avoid functions calls inside the loop as much as possible.

So in general, we want to do as little as possible inside the for loops, since these calls are very expensive in Python.

For example, we could write the above as:

```
import numpy as np
... some definitions ...
Z = sigma*np.sqrt(dt) *np.random.randn(N_timesteps -1, N_paths)
drift = (r-0.5*sigma**2)*dt
for i in range(1, N_timesteps):
        S[i] = S[i - 1] * np.exp(drift + Z[i-1])
```

So the inner loop has been "vectorized" away, the sampling is all done before the loops, etc.

If that is still not enough (which is very likely if you have many for loops), then you will need to perform some more "advanced" optimization. This usually comes down to taking the slowest piece of your code (almost always the for loop), put that in a function and write this function in e.g. Cython or Numba.

For example, with Numba we might define the for loop in a decorated function as:

```
from numba import autojit

@autojit
def my_loop(S, Z, drift):
    for i in range(1, S.shape[0]):
        for j in range(S.shape[1]):
            S[i, j] = S[i - 1, j] * np.exp(drift + Z[i-1, j])

... some definitions ...
Z = sigma*np.sqrt(dt) *np.random.randn(N_timesteps -1, N_paths)
drift = (r-0.5*sigma**2)*dt
my_loop(S, Z, drift)
```

The first time my_loop is called Numba will compile the function (takes half a second). This results in a function that runs much, much faster. It's a great way to optimize certain routines in your code.

But to really help you we'll need to see some code!