There are two general approaches for managing memory when solving most problems supported by the SDE engine:

Perform a traditional simulation to simulate the underlying variables of interest, specifically requesting and then manipulating the output arrays.

This approach is straightforward and the best choice for small
or medium-sized problems. Since its outputs are arrays, it is convenient
to manipulate simulated results in the MATLAB^{®} matrix-based language.
However, as the scale of the problem increases, the benefit of this
approach decreases, because the output arrays must store large quantities
of possibly extraneous information.

For example, consider pricing a European option in which the
terminal price of the underlying asset is the only value of interest.
To ease the memory burden of the traditional approach, reduce the
number of simulated periods specified by the required input `NPERIODS`

and
specify the optional input `NSTEPS`

. This enables
you to manage memory without sacrificing accuracy (see Optimizing Accuracy: About Solution Precision and Error).

In addition, simulation methods can determine the number of output arguments and allocate memory accordingly. Specifically, all simulation methods support the same output argument list:

[Paths,Times,Z]

where `Paths`

and `Z`

can
be large, three-dimensional time series arrays. However, the underlying
noise array is typically unnecessary, and is only stored if requested
as an output. In other words, `Z`

is stored only
at your request; do not request it if you do not need it.

If you need the output noise array `Z`

, but
do not need the `Paths`

time series array, then you
can avoid storing `Paths`

two ways:

It is best practice to use the

`~`

output argument placeholder. For example, use the following output argument list to store`Z`

and`Times`

, but not`Paths`

:[~,Times,Z]

Use the optional input flag

`StorePaths`

, which all simulation methods support. By default,`Paths`

is stored (`StorePaths`

=`true`

). However, setting`StorePaths`

to`false`

returns`Paths`

as an empty matrix.

Specify one or more end-of-period processing functions to manage and store only the information of interest, avoiding simulation outputs altogether.

This approach requires you to specify one or more end-of-period
processing functions, and is often the preferred approach for large-scale
problems. This approach allows you to avoid simulation outputs altogether.
Since no outputs are requested, the three-dimensional time series
arrays `Paths`

and `Z`

are not stored.

This approach often requires more effort, but is far more elegant and allows you to customize tasks and dramatically reduce memory usage. See Pricing Equity Options.

The following approaches improve performance when solving SDE problems:

**Specifying model parameters as traditional MATLAB arrays and functions, in various combinations.**This provides a flexible interface that can support virtually any general nonlinear relationship. However, while functions offer a convenient and elegant solution for many problems, simulations typically run faster when you specify parameters as double-precision vectors or matrices. Thus, it is a good practice to specify model parameters as arrays when possible.**Use models that have overloaded Euler simulation methods, when possible.**Using Brownian motion (`BM`

) and geometric Brownian motion (`GBM`

) models that provide overloaded Euler simulation methods take advantage of separable, constant-parameter models. These specialized methods are exceptionally fast, but are only available to models with constant parameters that are simulated without specifying end-of-period processing and noise generation functions.**Replace the simulation of a constant-parameter, univariate model derived from the**. Treat the multivariate model as a portfolio of univariate models. This increases the dimensionality of the model and enhances performance by decreasing the effective number of simulation trials.`SDEDDO`

class with that of a diagonal multivariate model### Note

This technique is applicable only to constant-parameter univariate models without specifying end-of-period processing and noise generation functions.

**Take advantage of the fact that simulation methods are designed to detect the presence of**A`NaN`

(not a number) conditions returned from end-of-period processing functions.`NaN`

represents the result of an undefined numerical calculation, and any subsequent calculation based on a`NaN`

produces another`NaN`

. This helps improve performance in certain situations. For example, consider simulating paths of the underlier of a knock-out barrier option (that is, an option that becomes worthless when the price of the underlying asset crosses some prescribed barrier). Your end-of-period function could detect a barrier crossing and return a`NaN`

to signal early termination of the current trial.

The simulation architecture does not, in general, simulate *exact* solutions
to any SDE. Instead, the simulation architecture provides a discrete-time
approximation of the underlying continuous-time process, a simulation
technique often known as a *Euler approximation*.

In the most general case, a given simulation derives directly
from an SDE. Therefore, the simulated discrete-time process approaches
the underlying continuous-time process only in the limit as the time
increment *dt* approaches zero. In other words, the
simulation architecture places more importance on ensuring that the
probability distributions of the discrete-time and continuous-time
processes are close, than on the pathwise proximity of the processes.

Before illustrating techniques to improve the approximation
of solutions, it is helpful to understand the source of error. Throughout
this architecture, all simulation methods assume that model parameters
are piecewise constant over any time interval of length *dt*.
In fact, the methods even evaluate dynamic parameters at the beginning
of each time interval and hold them fixed for the duration of the
interval. This sampling approach introduces *discretization
error*.

However, there are certain models for which the piecewise constant approach provides exact solutions:

Creating Brownian Motion (BM) Models with constant parameters, simulated by Euler approximation (

`simByEuler`

).Creating Geometric Brownian Motion (GBM) Models with constant parameters, simulated by closed-form solution (

`simBySolution`

).Creating Hull-White/Vasicek (HWV) Gaussian Diffusion Models with constant parameters, simulated by closed-form solution (

`simBySolution`

)

More generally, you can simulate the exact solutions for these models even if the parameters vary with time, if they vary in a piecewise constant way such that parameter changes coincide with the specified sampling times. However, such exact coincidence is unlikely; therefore, the previously discussed constant parameter condition is commonly used in practice.

One obvious way to improve accuracy involves sampling the discrete-time
process more frequently. This decreases the time increment (*dt*),
causing the sampled process to more closely approximate the underlying
continuous-time process. Although decreasing the time increment is
universally applicable, however, there is a tradeoff among accuracy,
run-time performance, and memory usage.

To manage this tradeoff, specify an optional input argument, `NSTEPS`

,
for all simulation methods. `NSTEPS`

indicates the
number of intermediate time steps within each time increment *dt*,
at which the process is sampled but not reported.

It is important and convenient at this point to emphasize the
relationship of the inputs `NSTEPS`

, `NPERIODS`

,
and `DeltaTime`

to the output vector `Times`

,
which represents the actual observation times at which the simulated
paths are reported.

`NPERIODS`

, a required input, indicates the number of simulation periods of length`DeltaTime`

, and determines the number of rows in the simulated three-dimensional`Paths`

time series array (if an output is requested).`DeltaTime`

is optional, and indicates the corresponding`NPERIODS`

-length vector of positive time increments between successive samples. It represents the familiar*dt*found in stochastic differential equations. If`DeltaTime`

is unspecified, the default value of 1 is used.`NSTEPS`

is also optional, and is only loosely related to`NPERIODS`

and`DeltaTime`

.`NSTEPS`

specifies the number of intermediate time steps within each time increment`DeltaTime`

.Specifically, each time increment

`DeltaTime`

is partitioned into`NSTEPS`

subintervals of length`DeltaTime`

/`NSTEPS`

each, and refines the simulation by evaluating the simulated state vector at`(NSTEPS - 1)`

intermediate times. Although the output state vector (if requested) is not reported at these intermediate times, this refinement improves accuracy by causing the simulation to more closely approximate the underlying continuous-time process. If`NSTEPS`

is unspecified, the default is`1`

(to indicate no intermediate evaluation).The output

`Times`

is an`NPERIODS + 1`

-length column vector of observation times associated with the simulated paths. Each element of`Times`

is associated with a corresponding row of`Paths`

.

The following example illustrates this intermediate
sampling by comparing the difference between a closed-form solution
and a sequence of Euler approximations derived from various values
of `NSTEPS`

.

Consider a univariate geometric Brownian motion (`GBM`

) model using `gbm`

with constant parameters:

$$d{X}_{t}=0.1{X}_{t}dt+0.4{X}_{t}d{W}_{t}.$$

Assume that the expected rate of return and volatility parameters are annualized, and that a calendar year comprises 250 trading days.

Simulate approximately four years of univariate prices for both the exact solution and the Euler approximation for various values of

`NSTEPS`

:nPeriods = 1000; dt = 1/250; obj = gbm(0.1,0.4,'StartState',100); rng(575,'twister') [X1,T1] = simBySolution(obj,nPeriods,'DeltaTime',dt); rng(575,'twister') [Y1,T1] = simByEuler(obj,nPeriods,'DeltaTime',dt); rng(575,'twister') [X2,T2] = simBySolution(obj,nPeriods,'DeltaTime',... dt,'nSteps',2); rng(575,'twister') [Y2,T2] = simByEuler(obj,nPeriods,'DeltaTime',... dt,'nSteps',2); rng(575,'twister') [X3,T3] = simBySolution(obj,nPeriods, 'DeltaTime',... dt,'nSteps',10); rng(575,'twister') [Y3,T3] = simByEuler(obj,nPeriods,'DeltaTime',... dt,'nSteps',10); rng(575,'twister') [X4,T4] = simBySolution(obj,nPeriods,'DeltaTime',... dt,'nSteps',100); rng(575,'twister') [Y4,T4] = simByEuler(obj,nPeriods,'DeltaTime',... dt,'nSteps',100);

Compare the error (the difference between the exact solution and the Euler approximation) graphically:

clf; plot(T1,X1 - Y1,'red') hold on; plot(T2,X2 - Y2,'blue') plot(T3,X3 - Y3,'green') plot(T4,X4 - Y4,'black') hold off xlabel('Time (Years)') ylabel('Price Difference') title('Exact Solution Minus Euler Approximation') legend({'# of Steps = 1' '# of Steps = 2' ... '# of Steps = 10' '# of Steps = 100'},... 'Location', 'Best') hold off

whos T X Y

As expected, the simulation error decreases as the number of
intermediate time steps increases. Because the intermediate states
are not reported, all simulated time series have the same number of
observations regardless of the actual value of `NSTEPS`

.

Furthermore, since the previously simulated exact solutions are correct for any number of intermediate time steps, additional computations are not needed for this example. In fact, this assessment is correct. The exact solutions are sampled at intermediate times to ensure that the simulation uses the same sequence of Gaussian random variates in the same order. Without this assurance, there is no way to compare simulated prices on a pathwise basis. However, there might be valid reasons for sampling exact solutions at closely spaced intervals, such as pricing path-dependent options.

`bm`

| `cev`

| `cir`

| `diffusion`

| `drift`

| `gbm`

| `heston`

| `hwv`

| `interpolate`

| `sde`

| `sdeddo`

| `sdeld`

| `sdemrd`

| `simByEuler`

| `simBySolution`

| `simBySolution`

| `simulate`

| `ts2func`

- Simulating Equity Prices
- Simulating Interest Rates
- Pricing American Basket Options by Monte Carlo Simulation
- Improving Performance of Monte Carlo Simulation with Parallel Computing
- Base SDE Models
- Drift and Diffusion Models
- Linear Drift Models
- Parametric Models