This example shows how to improve the performance of a Monte Carlo simulation using Parallel Computing Toolbox™.

Consider a geometric Brownian motion (GBM) process in which you want to incorporate alternative asset price dynamics. Specifically, suppose that you want to include a time-varying short rate and a volatility surface. The process to simulate is written as

$$dS(t)=r(t)S(t)dt+V(t,S(t))S(t)dW(t)$$

for stock price *S(t)*, rate of return
*r(t)*, volatility *V(t,S(t))*, and
Brownian motion *W(t)*. In this example, the rate of return is
a deterministic function of time and the volatility is a function of both time
and current stock price. Both the return and volatility are defined on a
discrete grid such that intermediate values is obtained by linear interpolation.
For example, such a simulation could be used to support the pricing of thinly
traded options.

To include a time series of riskless short rates, suppose that you derive the following deterministic short rate process as a function of time.

```
times = [0 0.25 0.5 1 2 3 4 5 6 7 8 9 10]; % in years
rates = [0.1 0.2 0.3 0.4 0.5 0.8 1.25 1.75 2.0 2.2 2.25 2.50 2.75]/100;
```

Suppose that you then derive the following volatility surface whose columns correspond to simple relative moneyness, or the ratio of the spot price to strike price, and whose rows correspond to time to maturity, or tenor.

```
surface = [28.1 25.3 20.6 16.3 11.2 6.2 4.9 4.9 4.9 4.9 4.9 4.9
22.7 19.8 15.4 12.6 9.6 6.7 5.2 5.2 5.2 5.2 5.2 5.2
21.7 17.6 13.7 11.5 9.4 7.3 5.7 5.4 5.4 5.4 5.4 5.4
19.8 16.4 12.9 11.1 9.3 7.6 6.2 5.6 5.6 5.6 5.6 5.6
18.6 15.6 12.5 10.8 9.3 7.8 6.6 5.9 5.9 5.9 5.9 5.9
17.4 13.8 11.7 10.8 9.9 9.1 8.5 7.9 7.4 7.3 7.3 7.3
17.1 13.7 12.0 11.2 10.6 10.0 9.5 9.1 8.8 8.6 8.4 8.0
17.5 13.9 12.5 11.9 11.4 10.9 10.5 10.2 9.9 9.6 9.4 9.0
18.3 14.9 13.7 13.2 12.8 12.4 12.0 11.7 11.4 11.2 11.0 10.8
19.2 19.6 14.2 13.9 13.4 13.0 13.2 12.5 12.1 11.9 11.8 11.4]/100;
tenor = [0 0.25 0.50 0.75 1 2 3 5 7 10]; % in years
moneyness = [0.25 0.5 0.75 0.8 0.9 1 1.10 1.25 1.50 2 3 5];
```

Set the simulation parameters. The following assumes that the price of the underlying asset is initially equal to the strike price and that the price of the underlying asset is simulated monthly for 10 years, or 120 months. As a simple illustration, 100 sample paths are simulated.

price = 100; strike = 100; dt = 1/12; nPeriods = 120; nTrials = 100;

For reproducibility, set the random number generator to its default, and draw the Gaussian random variates that drive the simulation. Generating the random variates is not necessary to incur the performance improvement of parallel computation, but doing so allows the resulting simulated paths to match those of the conventional (that is, non-parallelized) simulation. Moreover, generating independent Gaussian random variates as inputs also guarantees that all simulated paths are independent.

```
rng default
Z = randn(nPeriods,1,nTrials);
```

Create the return and volatility functions and the GBM model using `gbm`

. Notice that the rate of return
is a deterministic function of time, and therefore accepts simulation time as
its only input argument. In contrast, the volatility must account for the
moneyness and is a function of both time and stock price. Moreover, since the
volatility surface is defined as a function of time to maturity rather than
simulation time, the volatility function subtracts the current simulation time
from the last time at which the price process is simulated (10 years). This
ensures that as the simulation time approaches its terminal value, the time to
maturity of the volatility surface approaches zero. Although far more elaborate
return and volatility functions could be used if desired, the following assumes
simple linear interpolation.

```
mu = @(t) interp1(times,rates,t);
sigma = @(t,S) interp2(moneyness,tenor,surface,S/strike,tenor(end)-t);
mdl = gbm(mu,sigma,'StartState',price);
```

Simulate the paths of the underlying geometric Brownian motion without parallelization.

tStart = tic; paths = simBySolution(mdl,nPeriods,'nTrials',nTrials,'DeltaTime',dt,'Z',Z); time1 = toc(tStart);

Simulate the paths in parallel using a `parfor`

loop.
For users licensed to access the Parallel
Computing Toolbox, the
following code segment automatically creates a parallel pool using
the default local profile. If desired, this behavior can be changed
by first calling the `parpool`

function. If a parallel
pool is not already created, the following simulation will likely
be slower than the previous simulation without parallelization. In
this case, rerun the following simulation to assess the true benefits
of parallelization.

tStart = tic; parPaths = zeros(nPeriods+1,1,nTrials); parfor i = 1:nTrials parPaths(:,:,i) = simBySolution(mdl,nPeriods,'DeltaTime',dt,'Z',Z(:,:,i)); end time2 = toc(tStart);

If you examine any given path obtained without parallelization to the corresponding path with parallelization, you see that they are identical. Moreover, although relative performance varies, the results obtained with parallelization will generally incur a significant improvement. To assess the performance improvement, examine the runtime of each approach in seconds and speedup gained from simulating the paths in parallel.

time1 % elapsed time of conventional simulation, in seconds time2 % elapsed time of parallel simulation, in seconds speedup = time1/time2 % speedup factor

time1 = 6.1329 time2 = 2.5918 speedup = 2.3663

`bm`

| `cev`

| `cir`

| `diffusion`

| `drift`

| `gbm`

| `heston`

| `hwv`

| `interpolate`

| `parpool`

| `sde`

| `sdeddo`

| `sdeld`

| `sdemrd`

| `simByEuler`

| `simBySolution`

| `simBySolution`

| `simulate`

| `ts2func`

- Simulating Equity Prices
- Simulating Interest Rates
- Pricing American Basket Options by Monte Carlo Simulation
- Base SDE Models
- Drift and Diffusion Models
- Linear Drift Models
- Parametric Models