simulate

Simulate multivariate stochastic differential equations (SDEs) for SDE, BM, GBM, CEV, CIR, HWV, Heston, SDEDDO, SDELD, SDEMRD, Merton, or Bates models

Syntax

[Paths,Times,Z] = simulate(MDL)

[Paths,Times,Z] = simulate(___,Optional,Scheme)

Description

[Paths,Times,Z] = simulate(MDL) simulates NTrials sample paths of NVars correlated state variables, driven by NBrowns Brownian motion sources of risk over NPeriods consecutive observation periods, approximating continuous-time stochastic processes.

[Paths,Times,Z] = simulate(___,Optional,Scheme) adds optional inputs for Optional and Scheme.

The Optional input argument for simulate accepts any variable-length list of input arguments that the simulation method or function referenced by the SDE.Simulation parameter requires or accepts. It passes this input list directly to the appropriate SDE simulation method or user-defined simulation function.

The optional input Scheme lets you specify an approximation scheme used to simulate the sample paths and you can use this optional input with or without an Optional input argument.

Examples

Antithetic Sampling to a Path-Dependent Barrier Option

Consider a European up-and-in call option on a single underlying stock. The evolution of this stock's price is governed by a Geometric Brownian Motion (GBM) model with constant parameters:

Assume the following characteristics:

• The stock currently trades at 105.

• The stock pays no dividends.

• The stock volatility is 30% per annum.

• The option strike price is 100.

• The option expires in three months.

• The option barrier is 120.

• The risk-free rate is constant at 5% per annum.

The goal is to simulate various paths of daily stock prices, and calculate the price of the barrier option as the risk-neutral sample average of the discounted terminal option payoff. Since this is a barrier option, you must also determine if and when the barrier is crossed.

This example performs antithetic sampling by explicitly setting the Antithetic flag to true, and then specifies an end-of-period processing function to record the maximum and terminal stock prices on a path-by-path basis.

Create a GBM model using gbm.

barrier  	= 120;           % barrier
strike   	= 100;           % exercise price
rate     	= 0.05;          % annualized risk-free rate
sigma    	= 0.3;           % annualized volatility
nPeriods 	= 63;            % 63 trading days
dt       	= 1 / 252;       % time increment = 252 days
Time        = nPeriods * dt; % expiration time = 0.25 years
obj 	   = gbm(rate, sigma, 'StartState', 105);

Perform a small-scale simulation that explicitly returns two simulated paths.

rng('default')                % make output reproducible
[X, T] = obj.simBySolution(nPeriods, 'DeltaTime', dt, ...
				'nTrials', 2, 'Antithetic', true);

Perform antithetic sampling such that all primary and antithetic paths are simulated and stored in successive matching pairs. Odd paths (1,3,5,...) correspond to the primary Gaussian paths. Even paths (2,4,6,...) are the matching antithetic paths of each pair, derived by negating the Gaussian draws of the corresponding primary (odd) path. Verify this by examining the matching paths of the primary/antithetic pair.

plot(T, X(:,:,1), 'blue', T, X(:,:,2), 'red')
xlabel('Time (Years)'), ylabel('Stock Price'), ... 
			title('Antithetic Sampling')
legend({'Primary Path' 'Antithetic Path'}, ...
			'Location', 'Best')

To price the European barrier option, specify an end-of-period processing function to record the maximum and terminal stock prices. This processing function is accessible by time and state, and is implemented as a nested function with access to shared information that allows the option price and corresponding standard error to be calculated. For more information on using an end-of-period processing function, see Price European Stock Options Using Monte Carlo Simulation.

Simulate 200 paths using the processing function method.

rng('default')             % make output reproducible
barrier  = 120;            % barrier
strike   = 100;            % exercise price
rate     = 0.05;           % annualized risk-free rate
sigma    = 0.3;            % annualized volatility
nPeriods = 63;             % 63 trading days
dt       = 1 / 252;        % time increment = 252 days
Time     = nPeriods * dt;  % expiration time = 0.25 years
obj    = gbm(rate, sigma, 'StartState', 105);
nPaths = 200;         % # of paths = 100 sets of pairs
f      = Example_BarrierOption(nPeriods, nPaths);
simulate(obj, nPeriods, 'DeltaTime' , dt, ... 
			'nTrials', nPaths, 'Antithetic', true, ...
			'Processes', f.SaveMaxLast);

Approximate the option price with a 95% confidence interval.

optionPrice   = f.OptionPrice(strike, rate, barrier);
standardError = f.StandardError(strike, rate, barrier,...
			 			true);
lowerBound    = optionPrice - 1.96 * standardError;
upperBound    = optionPrice + 1.96 * standardError;

displaySummary(optionPrice, standardError, lowerBound, upperBound);

  Up-and-In Barrier Option Price:   6.6572
         Standard Error of Price:   0.7292
 Confidence Interval Lower Bound:   5.2280
 Confidence Interval Upper Bound:   8.0864

Utility Function

function displaySummary(optionPrice, standardError, lowerBound, upperBound)
fprintf('  Up-and-In Barrier Option Price: %8.4f\n', ...
			optionPrice);
fprintf('         Standard Error of Price: %8.4f\n', ...
			standardError);
fprintf(' Confidence Interval Lower Bound: %8.4f\n', ...
			lowerBound);
fprintf(' Confidence Interval Upper Bound: %8.4f\n', ...
			upperBound);
end

Specify Approximation Scheme with simulate Using a CIR Model

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (SDEMRD): $d{X}_t=S(t)[L(t)-X_t]dt+D(t,X_t^{\frac{1}{2}})V(t)dW$

where $$D$$ is a diagonal matrix whose elements are the square root of the corresponding element of the state vector.

Create a cir object to represent the model: $d{X}_t=0.2(0.1-X_t)dt+0.05X_t^{\frac{1}{2}}dW$.

cir_obj = cir(0.2, 0.1, 0.05)  % (Speed, Level, Sigma)

cir_obj = 
   Class CIR: Cox-Ingersoll-Ross
   ----------------------------------------
     Dimensions: State = 1, Brownian = 1
   ----------------------------------------
      StartTime: 0
     StartState: 1
    Correlation: 1
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
          Sigma: 0.05
          Level: 0.1
          Speed: 0.2

Use the optional name-value argument for 'Scheme' to specify a scheme to simulate the sample paths.

[paths,times] = simulate(cir_obj,10,'ntrials',4096,'scheme','quadratic-exponential');

Use simulate with Quasi-Monte Carlo Simulation with simByEuler Using a CIR Model

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (SDEMRD): $d{X}_t=S(t)[L(t)-X_t]dt+D(t,X_t^{\frac{1}{2}})V(t)dW$

where $$D$$ is a diagonal matrix whose elements are the square root of the corresponding element of the state vector.

Create a cir object to represent the model: $d{X}_t=0.2(0.1-X_t)dt+0.05X_t^{\frac{1}{2}}dW$.

cir_obj = cir(0.2, 0.1, 0.05)  % (Speed, Level, Sigma)

cir_obj = 
   Class CIR: Cox-Ingersoll-Ross
   ----------------------------------------
     Dimensions: State = 1, Brownian = 1
   ----------------------------------------
      StartTime: 0
     StartState: 1
    Correlation: 1
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
          Sigma: 0.05
          Level: 0.1
          Speed: 0.2

Use optional name-value inputs for the simByEuler method that you can call through simulate interface using the simulate 'Scheme' for 'quadratic-exponential'. The optional inputs for simByEuler define a quasi-Monte Carlo simulation using the name-value arguments for 'MonteCarloMethod','QuasiSequence', and 'BrownianMotionMethod'.

[paths,time] = simulate(cir_obj,10,'ntrials',4096,'montecarlomethod','quasi','quasisequence','sobol','BrownianMotionMethod','brownian-bridge','scheme','quadratic-exponential');

Use simulate with Quasi-Monte Carlo Simulation with simByMilstein Using a CIR Model

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (SDEMRD): $d{X}_t=S(t)[L(t)-X_t]dt+D(t,X_t^{\frac{1}{2}})V(t)dW$

where $$D$$ is a diagonal matrix whose elements are the square root of the corresponding element of the state vector.

Create a cir object to represent the model: $d{X}_t=0.2(0.1-X_t)dt+0.05X_t^{\frac{1}{2}}dW$.

cir_obj = cir(0.2, 0.1, 0.05)  % (Speed, Level, Sigma)

cir_obj = 
   Class CIR: Cox-Ingersoll-Ross
   ----------------------------------------
     Dimensions: State = 1, Brownian = 1
   ----------------------------------------
      StartTime: 0
     StartState: 1
    Correlation: 1
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
          Sigma: 0.05
          Level: 0.1
          Speed: 0.2

Use optional name-value inputs for the simByMilstein method that you can call through simulate interface using the simulate 'Scheme' for 'milstein'. The optional inputs for simByMilstein define a quasi-Monte Carlo simulation using the name-value arguments for 'MonteCarloMethod','QuasiSequence', and 'BrownianMotionMethod'.

[paths,time] = simulate(cir_obj,10,'ntrials',4096,'montecarlomethod','quasi','quasisequence','sobol','BrownianMotionMethod','brownian-bridge','scheme','milstein');

Use simulate with Quasi-Monte Carlo Simulation with simByMilstein2 Using a CIR Model

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (SDEMRD): $d{X}_t=S(t)[L(t)-X_t]dt+D(t,X_t^{\frac{1}{2}})V(t)dW$

where $$D$$ is a diagonal matrix whose elements are the square root of the corresponding element of the state vector.

Create a cir object to represent the model: $d{X}_t=0.2(0.1-X_t)dt+0.05X_t^{\frac{1}{2}}dW$.

cir_obj = cir(0.2, 0.1, 0.05)  % (Speed, Level, Sigma)

cir_obj = 
   Class CIR: Cox-Ingersoll-Ross
   ----------------------------------------
     Dimensions: State = 1, Brownian = 1
   ----------------------------------------
      StartTime: 0
     StartState: 1
    Correlation: 1
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
          Sigma: 0.05
          Level: 0.1
          Speed: 0.2

Use optional name-value inputs for the simByMilstein2 method that you can call through simulate interface using the simulate 'Scheme' for 'milstein2'. The optional inputs for simByMilstein2 define a quasi-Monte Carlo simulation using the name-value arguments for 'MonteCarloMethod','QuasiSequence', and 'BrownianMotionMethod'.

[paths,time] = simulate(cir_obj,10,'ntrials',4096,'montecarlomethod','quasi','quasisequence','sobol','BrownianMotionMethod','brownian-bridge','scheme','milstein2');

Input Arguments

MDL — Stochastic differential equation model

object

Stochastic differential equation model, specified as an sde, bates, bm, gbm, cev, cir, hwv, heston,merton sdeddo, sdeld, or sdemrd object.

Data Types: object

Optional — Any variable-length list of input arguments that the simulation method or function referenced by the SDE.Simulation parameter requires or accepts

input arguments

(Optional) Any variable-length list of input arguments that the simulation method or function referenced by the SDE.Simulation parameter requires or accepts, specified as a variable-length list of input arguments. This input list is passed directly to the appropriate SDE simulation method or user-defined simulation function.

Data Types: double

Scheme — Approximation scheme used to simulate the sample paths

character vector with values 'euler', 'solution','quadratic-exponential', 'transition', 'milstein', or 'milstein2' | string with values "euler", "solution","quadratic-exponential", "transition", "milstein", or "milstein2"

(Optional) Approximation scheme used to simulate the sample paths, specified as a character vector or string. When you specify Scheme, the MDL simulation method using the corresponding scheme is used instead of the one defined in the MDL.Simulation property which is the simByEuler method by default.

Note

The supported Scheme depends on which SDE object you use with the simulate method as follows:

simulate Scheme Value	Supported Objects
"euler"	sde, sdeddo, sdeld, sdemrd, bm, cev, gbm, merton, hwv, cir, heston, or bates
"solution"	gbm, merton, or hwv
"transition"	cir, heston, or bates
"quadratic-exponential"	cir, heston, or bates
"milstein"	sdeddo, sdeld, sdemrd, bm, cev, gbm, merton, hwv, cir, heston, or bates
"milstein2"	sdeddo, sdeld, sdemrd, bm, cev, gbm, merton, hwv, cir, heston, or bates

Data Types: char | string

Output Arguments

Paths — Three-dimensional time series array, consisting of simulated paths of correlated state variables

array

Three-dimensional time series array, consisting of simulated paths of correlated state variables, returned as an (NPeriods + 1)-by-NVars-by-NTrials array.

For a given trial, each row of Paths is the transpose of the state vector X_t at time t.

Times — Observation times associated with the simulated paths

column vector

Observation times associated with the simulated paths, returned as an (NPeriods + 1)-by-1 column vector.

Z — Three-dimensional time series array of dependent random variates used to generate the Brownian motion vector (Wiener processes)

array | matrix | table | timetable

Three-dimensional time series array of dependent random variates used to generate the Brownian motion vector (Wiener processes) that drove the simulated results found in Paths, returned as an NTimes-by-NBrowns-by-NTrials array.

NTimes is the number of time steps at which the simulate function samples the state vector. NTimes includes intermediate times designed to improve accuracy, which simulate does not necessarily report in the Paths output time series.

More About

Antithetic Sampling

Simulation methods allow you to specify a popular variance reduction technique called antithetic sampling.

This technique attempts to replace one sequence of random observations with another of the same expected value, but smaller variance. In a typical Monte Carlo simulation, each sample path is independent and represents an independent trial. However, antithetic sampling generates sample paths in pairs. The first path of the pair is referred to as the primary path, and the second as the antithetic path. Any given pair is independent of any other pair, but the two paths within each pair are highly correlated. Antithetic sampling literature often recommends averaging the discounted payoffs of each pair, effectively halving the number of Monte Carlo trials.

This technique attempts to reduce variance by inducing negative dependence between paired input samples, ideally resulting in negative dependence between paired output samples. The greater the extent of negative dependence, the more effective antithetic sampling is.

Algorithms

This function simulates any vector-valued SDE of the form:

$$d{X}_t=F(t,X_t)dt+G(t,X_t)d{W}_t$$

(1)

where:

• X is an NVars-by-1 state vector of process variables (for example, short rates or equity prices) to simulate.

• W is an NBrowns-by-1 Brownian motion vector.

• F is an NVars-by-1 vector-valued drift-rate function.

• G is an NVars-by-NBrowns matrix-valued diffusion-rate function.

References

[1] Ait-Sahalia, Y. “Testing Continuous-Time Models of the Spot Interest Rate.” The Review of Financial Studies, Spring 1996, Vol. 9, No. 2, pp. 385–426.

[2] Ait-Sahalia, Y. “Transition Densities for Interest Rate and Other Nonlinear Diffusions.” The Journal of Finance, Vol. 54, No. 4, August 1999.

[3] Glasserman, P. Monte Carlo Methods in Financial Engineering, New York, Springer-Verlag, 2004.

[4] Hull, J. C. Options, Futures, and Other Derivatives, 5th ed. Englewood Cliffs, NJ: Prentice Hall, 2002.

[5] Johnson, N. L., S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions, Vol. 2, 2nd ed. New York, John Wiley & Sons, 1995.

[6] Shreve, S. E. Stochastic Calculus for Finance II: Continuous-Time Models, New York: Springer-Verlag, 2004.

Version History

Introduced in R2008a

R2023b: Added "milstein2" value to Scheme name-value argument

When the MDL input argument is a sdeddo, sdeld, sdemrd, bm, cev, gbm, merton, hwv, cir, heston, or bates object, you can use a Scheme value of "milstein2".

Use the simByMilstein method to approximate a numerical solution of a stochastic differential equation.

R2023a: Added "milstein" value to Scheme name-value argument

When the MDL input argument is a sdeddo, sdeld, sdemrd, bm, cev, gbm, merton, hwv, cir, heston, or bates object, you can use a Scheme value of "milstein".

Use the simByMilstein method to approximate a numerical solution of a stochastic differential equation.

R2022b: Approximation scheme to simulate the sample paths

When you specify Scheme, the MDL simulation method using the corresponding scheme is used instead of the one defined in the MDL.Simulation property.