Vectorizing Brownian Bridge Simulations Across Paths

# 68492


# Efficient method for expanding 1 sim routine to the number of simulations? Brownian Bridge used with multiple underlying assets in a MC simulation,












I believe this is a (fairly) simple question for those familiar with quantitative finance and MC/QMC methods of pricing complex options. Or potentially its just a simple Python loop vectorization question, with no knowledge of quantitative finance needed.

So I'm trying to understand the Brownian-Bridge technique. I work best of code examples, of which a great one is here (credit to author Kenta Oono): https://gist.github.com/delta2323/6bb572d9473f3b523e6e - in the comments is a correction to the routine, which appears correctly implemented from what I've read on the path construction. Here's that rewritten so it makes more sense for what I'm doing. This is for simulating a basket of 12 underlying assets over an averaging period of 21 days. Note it is written for 1 simulation path only:

```
import numpy as np
from matplotlib import pyplot
import timeit

steps = 21
underlyings = 12

#seed = 0 # fix your seed for calcing Greeks
#np.random.seed(seed)

def sample_path_batch(underlyings, steps):
    dt = 1.0 / (steps-1)
    dt_sqrt = np.sqrt(dt)
    B = np.empty((underlyings, steps), dtype=float)
    B[:, 0] = 0
    for n in range(steps - 2):
        t = n * dt
        xi = np.random.randn(underlyings) * dt_sqrt
        B[:, n + 1] = B[:, n] * (1 - dt / (1 - t)) + xi
    B[:, -1] = 0 # set last step to 0
    return B

start_time = timeit.default_timer()
B = sample_path_batch(underlyings, steps)
print('\n' + 'Run time for 1 simulation steps * underlyings: ', int(timeit.default_timer() - start_time), ' seconds')
pyplot.plot(B.T)
pyplot.show()
```

So what's an efficient way (rather than looping over this single routine) to construct this one path to an arbitrary number of simulations? I'm using Sobol so call it 1024 simulations as an example. Although I expect to using far more, hence the desire to remove this loop, as I've put everything into another j loop and made it feed into a 3D NumPy array, but this slows down quickly (with the number of simulations):

```
#Now my inefficient way to generate multiple simulations of the Brownian Bridges:
sims = pow(2,17) # 131,072 simulations
def sample_path_batches(underlyings, steps, sims):
    dt = 1.0 / (steps-1)
    dt_sqrt = np.sqrt(dt)
    B = np.empty((underlyings, steps, sims), dtype=float)
    B[:,:, 0] = 0
    for n in range(steps - 2):
        for j in range(sims):
            t = n * dt
            xi = np.random.randn(underlyings) * dt_sqrt
            B[:, n + 1, j] = B[:, n, j] * (1 - dt / (1 - t)) + xi
            B[:, -1, j] = 0 # set last step to 0
    return B

start_time = timeit.default_timer()
B = sample_path_batches(underlyings, steps, sims)
print('\n' + 'Run time for ', sims, ' simulation steps * underlyings: ', 
int(timeit.default_timer() - start_time), ' seconds')
```

The above takes 13 seconds with pow(2,17) = 131,072 simulations, hence why I'd like to speed this up. Any suggestions are better than none. I think it's more of a general Python question than really a quantitative finance question. Sure I could do the routine in Cython and make the loops parallel but I'm just looking for an efficient way to do this inside of Python and the standard libraries.

## Answer by Matt (score 0, accepted)

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

Potential solution based on the comment by will, which initially (before I changed the np.random.randn to (underlyings, sims)) all the "random" numbers were identical across the simulations. Anyhow this "appears" to work very fast, although the matrices of random numbers (after fixing the seed) are not the same (j loop vs. this vectorized modification). Unfortunately, I'm not competent at plotting 3D arrays with matplotlib. They still appear to be Brownian Bridges when I plot them individually `pyplot.plot(B[:,:,sim#].T); pyplot.show()` and take stats on each simulation, although they differ from the j loop approach (likely, because in the loop, I have iteration-wise random numbers and here, well, they are generated in 2D). They still "appear" to be valid Brownian Bridges, but I did not replicate the same sequence the j loop was producing (i.e. changed the sequence of shock generation). I obviously can't view 131,072 simulations of 12 underlyings across 21 paths 1 by 1... Nonetheless, here's the modification, I suppose I'll know how "good" it is once I try implementing it (not sure if I should be using Mersenne Twister (as shown) or RQMC points for this, but here it is, probably will come in handy for others) - now if I could only vectorize the Brownian loop away, well, maybe on another day:

```
import numpy as np
from matplotlib import pyplot

steps = 21
underlyings = 12
sims = 131072

def sample_path_batches(underlyings, steps, sims):
    dt = 1.0 / (steps-1)
    dt_sqrt = np.sqrt(dt)
    B = np.empty((underlyings, steps, sims), dtype=float)
    B[:, 0, :] = 0 # set first step to 0
    for n in range(steps - 2):
# =============================================================================
# Slow j loop before vectorization:
#
#         for j in range(sims):
#             t = n * dt
#             xi = np.random.randn(underlyings) * dt_sqrt
#             B[:, n + 1, j] = B[:, n, j] * (1 - dt / (1 - t)) + xi
#             B[:, -1, j] = 0 # set last step to 0
# =============================================================================

# Change the generation of random numbers to be done in 1 step of the Brownian Bridge,
# for a huge speedup in computational time; note the random numbers differ from the j loop

        t = n * dt
        xi = np.random.randn(underlyings, sims) * dt_sqrt
        B[:, n + 1, :] = B[:, n, :] * (1 - dt / (1 - t)) + xi
        B[:, -1, :] = 0 # set last step to 0
    return B
```

Run time for 1 simulation steps * underlyings: 0.022 seconds

Run time for 131072 simulation steps * underlyings: 1.131 seconds