Simulating Equity Prices

Simulating Multidimensional Market Models

This example compares alternative implementations of a separable multivariate geometric Brownian motion process that is often referred to as a multidimensional market model. It simulates sample paths of an equity index portfolio using sde, sdeddo, sdeld, cev, and gbm objects.

The market model to simulate is:

$$d{X}_t=\mu X_{t}dt+D(X_t)\sigma d{W}_t$$

(1)

where:

• μ is a diagonal matrix of expected index returns.

• D is a diagonal matrix with X_t along the diagonal.

• σ is a diagonal matrix of standard deviations of index returns.

Represent Market Models Using SDE Models

Create an sde object using sde to represent the equity market model.

Load the Data_GlobalIdx2 data set from the matlabroot/examples/finance/data directory.

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
Dataset.NIK Dataset.FTSE Dataset.SP];

Convert daily prices to returns.

returns =  tick2ret(prices);

Compute data statistics to input to simulation methods.

nVariables  = size(returns, 2);
expReturn   = mean(returns);
sigma       = std(returns);
correlation = corrcoef(returns);
t           = 0;
X           = 100;
X           = X(ones(nVariables,1));

Create simple anonymous drift and diffusion functions accessible by$\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right).$

F = @(t,X) diag(expReturn) * X;
G = @(t,X) diag(X) * diag(sigma);

Use these functions with sde to create an sde object to represent the market model in

SDE = sde(F, G, 'Correlation', correlation, 'StartState', X)

SDE = 
   Class SDE: Stochastic Differential Equation
   -------------------------------------------
     Dimensions: State = 6, Brownian = 6
   -------------------------------------------
      StartTime: 0
     StartState: 100 (6x1 double array) 
    Correlation: 6x6 double array 
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler

The sde object requires additional information to determine the dimensionality of the model, because the functions passed to the sde object are known only by their $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$ interface. In other words, the sde object requires only two inputs: a drift-rate function and a diffusion-rate function, both accessible by passing the sample time and the corresponding state vector $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$.

In this case, this information is insufficient to determine unambiguously the dimensionality of the state vector and Brownian motion. You resolve the dimensionality by specifying an initial state vector, StartState. The SDE engine has assigned the default simulation method, simByEuler, to the Simulation parameter.

Represent Market Models Using SDEDDO Models

Create an sdeddo object using sdeddo to represent the market model in

Load the Data_GlobalIdx2 data set from the matlabroot/examples/finance/data directory.

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
Dataset.NIK Dataset.FTSE Dataset.SP];

Convert daily prices to returns.

returns =  tick2ret(prices);

Compute data statistics to input to simulation methods.

nVariables  = size(returns, 2);
expReturn   = mean(returns);
sigma       = std(returns);
correlation = corrcoef(returns);

Create drift and diffusion objects using drift and diffusion.

F = drift(zeros(nVariables,1), diag(expReturn))

F = 
   Class DRIFT: Drift Rate Specification  
   -------------------------------------  
      Rate: drift rate function F(t,X(t)) 
         A: 6x1 double array
         B: 6x6 diagonal double array

G = diffusion(ones(nVariables,1), diag(sigma))

G = 
   Class DIFFUSION: Diffusion Rate Specification 
   --------------------------------------------- 
       Rate: diffusion rate function G(t,X(t))  
      Alpha: 6x1 double array
      Sigma: 6x6 diagonal double array

Pass the drift and diffusion objects to sdeddo.

SDEDDO = sdeddo(F, G, 'Correlation', correlation, ...
'StartState', 100)

SDEDDO = 
   Class SDEDDO: SDE from Drift and Diffusion Objects
   --------------------------------------------------
     Dimensions: State = 6, Brownian = 6
   --------------------------------------------------
      StartTime: 0
     StartState: 100 (6x1 double array) 
    Correlation: 6x6 double array 
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
              A: 6x1 double array
              B: 6x6 diagonal double array
          Alpha: 6x1 double array
          Sigma: 6x6 diagonal double array

The sdeddo object requires two input objects that provide more information than the two functions from step 4 of Represent Market Models Using SDE Models. Thus, the dimensionality is more easily resolved. In fact, the initial price of each index is as a scalar (StartState = 100). This is in contrast to the sde object, which required an explicit state vector to uniquely determine the dimensionality of the problem.

Once again, the class of each object is clearly identified, and parameters display like fields of a structure. In particular, the Rate parameter of drift and diffusion objects is identified as a callable function of time and state, $\mathit{F}\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$and $\mathit{G}\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$, respectively. The additional parameters, A, B, Alpha, and Sigma, are arrays of appropriate dimension, indicating static (non-time-varying) parameters. In other words, $\mathit{A}\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$, $\mathit{B}\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$, $\mathrm{Alpha}\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$, and $\mathrm{Sigma}\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$are constant functions of time and state.

Represent Market Models Using SDELD, CEV, and GBM Objects

Create sdeld, cev, and gbm objects to represent the market model in

Load the Data_GlobalIdx2 data set from the matlabroot/examples/finance/data directory.

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
    Dataset.NIK Dataset.FTSE Dataset.SP];

Convert daily prices to returns.

returns =  tick2ret(prices);

Compute data statistics to input to simulation methods.

nVariables  = size(returns, 2);
expReturn   = mean(returns);
sigma       = std(returns);
correlation = corrcoef(returns);
t           = 0;
X           = 100;
X           = X(ones(nVariables,1));

Create an sdeld object using sdeld.

SDELD = sdeld(zeros(nVariables,1), diag(expReturn), ...
    ones(nVariables,1), diag(sigma),'Correlation', ... 
    correlation, 'StartState', X)

SDELD = 
   Class SDELD: SDE with Linear Drift
   ----------------------------------------
     Dimensions: State = 6, Brownian = 6
   ----------------------------------------
      StartTime: 0
     StartState: 100 (6x1 double array) 
    Correlation: 6x6 double array 
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
              A: 6x1 double array
              B: 6x6 diagonal double array
          Alpha: 6x1 double array
          Sigma: 6x6 diagonal double array

Create a cev object using cev.

CEV = cev(diag(expReturn), ones(nVariables,1), ... 
    diag(sigma), 'Correlation', correlation, ...
    'StartState', X)

CEV = 
   Class CEV: Constant Elasticity of Variance
   ------------------------------------------
     Dimensions: State = 6, Brownian = 6
   ------------------------------------------
      StartTime: 0
     StartState: 100 (6x1 double array) 
    Correlation: 6x6 double array 
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
         Return: 6x6 diagonal double array
          Alpha: 6x1 double array
          Sigma: 6x6 diagonal double array

Create a gbm object using gbm.

GBM = gbm(diag(expReturn), diag(sigma), 'Correlation', ...
    correlation, 'StartState', X)

GBM = 
   Class GBM: Generalized Geometric Brownian Motion
   ------------------------------------------------
     Dimensions: State = 6, Brownian = 6
   ------------------------------------------------
      StartTime: 0
     StartState: 100 (6x1 double array) 
    Correlation: 6x6 double array 
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
         Return: 6x6 diagonal double array
          Sigma: 6x6 diagonal double array

Note the succession of interface restrictions:

• sdeld objects require you to specify A, B, Alpha, and Sigma.

• cev objects require you to specify Return, Alpha, and Sigma.

• gbm objects require you to specify only Return and Sigma.

However, all three objects represent the same multidimensional market model.

Also, cev and gbm objects display the underlying parameter B derived from the sdeld object as Return. This is an intuitive name commonly associated with equity models.

Simulate Equity Markets Using the Default Simulate Method

Load the Data_GlobalIdx2 data set from the matlabroot/examples/finance/data directory and use sde to specify the SDE model as in

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
    Dataset.NIK Dataset.FTSE Dataset.SP];

returns =  tick2ret(prices);

nVariables  = size(returns,2);
expReturn   = mean(returns);
sigma       = std(returns);
correlation = corrcoef(returns);
t           = 0;
X           = 100;
X           = X(ones(nVariables,1));

F = @(t,X) diag(expReturn)* X;
G = @(t,X) diag(X) * diag(sigma);

SDE = sde(F, G, 'Correlation', ...
    correlation, 'StartState', X);

Simulate a single path of correlated equity index prices over one calendar year (defined as approximately 250 trading days) using the default simulate method.

nPeriods = 249;      % # of simulated observations
dt       =   1;      % time increment = 1 day
rng(142857,'twister')
[S,T] = simulate(SDE, nPeriods, 'DeltaTime', dt);

whos S

  Name        Size            Bytes  Class     Attributes

  S         250x6             12000  double              

The output array S is a 250-by-6 = (NPeriods + 1-by-nVariables-by-1) array with the same initial value, 100, for all indices. Each row of S is an observation of the state vector ${\mathit{X}}_{\mathit{t}}^{}$at time t.

Plot the simulated paths.

plot(T, S), xlabel('Trading Day'), ylabel('Price')
title('Single Path of Multi-Dimensional Market Model')
legend({'Canada' 'France' 'Germany' 'Japan' 'UK' 'US'}, ... 
    'Location', 'Best')

Simulate Equity Markets Using simByEuler Method

Because simByEuler is a valid simulation method, you can call it directly, overriding the Simulation parameter's current method or function (which in this case is simByEuler). The following statements produce the same price paths as generated in Simulate Equity Markets Using the Default Simulate Method.

Load the Data_GlobalIdx2 data set from the matlabroot/examples/finance/data directory and use sde to specify the SDE model as in

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
    Dataset.NIK Dataset.FTSE Dataset.SP];

returns =  tick2ret(prices);

nVariables  = size(returns,2);
expReturn   = mean(returns);
sigma       = std(returns);
correlation = corrcoef(returns);
t           = 0;
X           = 100;
X           = X(ones(nVariables,1));

F = @(t,X) diag(expReturn)* X;
G = @(t,X) diag(X) * diag(sigma);

SDE = sde(F, G, 'Correlation', ...
    correlation, 'StartState', X);

Simulate a single path using simByEuler.

nPeriods = 249;      % # of simulated observations
dt       =   1;      % time increment = 1 day
rng(142857,'twister')
[S,T] = simByEuler(SDE, nPeriods, 'DeltaTime', dt);

Simulate 10 trials with the same initial conditions, and examine S.

rng(142857,'twister')
[S,T] = simulate(SDE, nPeriods, 'DeltaTime', dt, 'nTrials', 10);

whos S

  Name        Size               Bytes  Class     Attributes

  S         250x6x10            120000  double              

Now the output array S is an NPeriods + 1-by-nVariables-by-NTrials time series array.

Plot the first paths.

plot(T, S(:,:,1)), xlabel('Trading Day'), ylabel('Price')
title('First Path of Multi-Dimensional Market Model')
legend({'Canada' 'France' 'Germany' 'Japan' 'UK' 'US'},...
    'Location', 'Best')

The first realization of

is identical to the paths in the plot.

Simulate Equity Markets Using GBM Simulation Methods

Consider simulation using GBM simulation methods. Separable GBM models have two specific simulation methods:

• An overloaded Euler simulation method, simByEuler, designed for optimal performance

• A function, simBySolution, provides an approximate solution of the underlying stochastic differential equation, designed for accuracy

Load the Data_GlobalIdx2 data set from the matlabroot/examples/finance/data directory and use sde to specify the SDE model as in Represent Market Models Using SDE Models, and the GBM model as in Represent Market Models Using SDELD, CEV, and GBM Objects.

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
    Dataset.NIK Dataset.FTSE Dataset.SP];

returns =  tick2ret(prices);

nVariables  = size(returns,2);
expReturn   = mean(returns);
sigma       = std(returns);
correlation = corrcoef(returns);
t           = 0;
X           = 100;
X           = X(ones(nVariables,1));

F = @(t,X) diag(expReturn)* X;
G = @(t,X) diag(X) * diag(sigma);

SDE = sde(F, G, 'Correlation', ...
    correlation, 'StartState', X);

GBM = gbm(diag(expReturn),diag(sigma), 'Correlation', ...
    correlation, 'StartState', X);

To illustrate the performance benefit of the overloaded Euler approximation method, increase the number of trials to 10000.

nPeriods = 249;      % # of simulated observations
dt       =   1;      % time increment = 1 day
rng(142857,'twister')
[X,T] = simulate(GBM, nPeriods, 'DeltaTime', dt, ...
    'nTrials', 10000);

whos X

  Name        Size                     Bytes  Class     Attributes

  X         250x6x10000            120000000  double              

The output X is a much larger time series array.

Using this sample size, plot a histogram to examine the terminal distribution of Canada's TSX Composite to verify qualitatively the lognormal character of the data.

histogram(squeeze(X(end,1,:)), 30), xlabel('Price'), ylabel('Frequency')
title('Histogram of Prices after One Year: Canada (TSX Composite)')

Simulate 10 trials of the solution and plot the first trial.

rng(142857,'twister')
[S,T] = simulate(SDE, nPeriods, 'DeltaTime', dt, 'nTrials', 10);
rng(142857,'twister')
[X,T] = simBySolution(GBM, nPeriods,...
    'DeltaTime', dt, 'nTrials', 10);
subplot(2,1,1)
plot(T, S(:,:,1)), xlabel('Trading Day'),ylabel('Price')
title('1st Path of Multi-Dim Market Model:Euler Approximation')
subplot(2,1,2)
plot(T, X(:,:,1)), xlabel('Trading Day'),ylabel('Price')
title('1st Path of Multi-Dim Market Model:Analytic Solution')

In this example, all parameters are constants, and simBySolution does indeed sample the exact solution. The details of a single index for any given trial show that the price paths of the Euler approximation and the exact solution are close, but not identical.

Plot the difference between the two methods.

subplot(1,1,1)
plot(T, S(:,1,1) - X(:,1,1), 'blue'), grid('on')
xlabel('Trading Day'), ylabel('Price Difference')
title('Euler Approx Minus Exact Solution:Canada(TSX Composite)')

The simByEuler Euler approximation literally evaluates the stochastic differential equation directly from the equation of motion, for some suitable value of the dt time increment. This simple approximation suffers from discretization error. This error is attributed to the discrepancy between the choice of the dt time increment and what in theory is a continuous-time parameter.

The discrete-time approximation improves as DeltaTime approaches zero. The Euler method is often the least accurate and most general method available. All models shipped in the simulation suite have this method.

In contrast, the simBySolution method provides a more accurate description of the underlying model. This method simulates the price paths by an approximation of the closed-form solution of separable models. Specifically, it applies a Euler approach to a transformed process, which in general is not the exact solution to this GBM model. This is because the probability distributions of the simulated and true state vectors are identical only for piecewise constant parameters.

When all model parameters are piecewise constant over each observation period, the simulated process is exact for the observation times at which the state vector is sampled. Since all parameters are constants in this example, simBySolution does indeed sample the exact solution.

For an example of how to use simBySolution to optimize the accuracy of solutions, see Optimizing Accuracy: About Solution Precision and Error.

Induce Dependence and Correlation Between States

This example illustrates two techniques that induce dependence between individual elements of a state vector. It also illustrates the interaction between Sigma and Correlation.

The first technique generates correlated Gaussian variates to form a Brownian motion process with dependent components. These components are then weighted by a diagonal volatility or exposure matrix Sigma.

The second technique generates independent Gaussian variates to form a standard Brownian motion process, which is then weighted by the lower Cholesky factor of the desired covariance matrix. Although these techniques can be used on many models, the relationship between them is most easily illustrated by working with a separable GBM model (see Simulate Equity Markets Using GBM Simulation Methods).

The market model to simulate is:

where μ is a diagonal matrix.

Load the data set.

load Data_GlobalIdx2
prices  = [Dataset.TSX Dataset.CAC Dataset.DAX ...
    Dataset.NIK Dataset.FTSE Dataset.SP];

Convert the daily prices to returns using tick2ret.

returns =  tick2ret(prices);

Specify Sigma and Correlation using the first technique where you specify Sigma as a diagonal matrix of asset return standard deviations.

expReturn   = diag(mean(returns));  % expected return vector
sigma       = diag(std(returns));   % volatility of returns

Specify Correlation as the sample correlation matrix of those returns. In this case, the components of the Brownian motion are dependent.

correlation = corrcoef(returns);    
GBM1        = gbm(expReturn,sigma,'Correlation',...
                  correlation);

Specify Sigma and Correlation using the second technique where you specify Sigma as the lower Cholesky factor of the asset return covariance matrix.

covariance = cov(returns);
sigma      = cholcov(covariance)';

Set Correlation to an identity matrix.

GBM2       = gbm(expReturn,sigma);

Here, sigma captures both the correlation and magnitude of the asset return uncertainty. In contrast to the first technique, the components of the Brownian motion are independent. Also, this technique accepts the default assignment of an identity matrix to Correlation, and is more straightforward.

Simulate a single trial of 1000 observations (roughly four years of daily data) using both techniques. By default, all state variables start at 1.

rng(22814,'twister')
[X1,T] = simByEuler(GBM1,1000);  % correlated Brownian motion
rng(22814,'twister')
[X2,T] = simByEuler(GBM2,1000);  % standard Brownian motion

When based on the same initial random number state, each technique generates identical asset price paths.

subplot(2,1,1)
plot(T, X1)
title('Sample Paths from Correlated Brownian Motion')
ylabel('Asset Price')
subplot(2,1,2)
plot(T, X2)
title('Sample Paths from Standard Brownian Motion')
xlabel('Trading Day')
ylabel('Asset Price')

Dynamic Behavior of Market Parameters

As discussed in Creating SDE Objects, object parameters may be evaluated as if they are MATLAB® functions accessible by a common interface. This accessibility provides the impression of dynamic behavior regardless of whether the underlying parameters are truly time-varying. Furthermore, because parameters are accessible by a common interface, seemingly simple linear constructs may in fact represent complex, nonlinear designs.

For example, consider a univariate geometric Brownian motion (GBM) model of the form:

In this model, the return, μ(t), and volatility,σ(t), are dynamic parameters of time alone. However, when creating a gbm object to represent the underlying model, such dynamic behavior must be accessible by the common $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$interface. This reflects the fact that GBM models are restricted parameterizations that derive from the general SDE class.

As a convenience, you can specify parameters of restricted models, such as GBM models, as traditional MATLAB® arrays of appropriate dimension. In this case, such arrays represent a static special case of the more general dynamic situation accessible by the $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$ interface.

Moreover, when you enter parameters as functions, object constructors can verify that they return arrays of correct size by evaluating them at the initial time and state. Otherwise, object constructors have no knowledge of any particular functional form.

The following example illustrates a technique that includes dynamic behavior by mapping a traditional MATLAB® time series array to a callable function with a $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$ interface. It also compares the function with an otherwise identical model with constant parameters.

Because time series arrays represent dynamic behavior that must be captured by functions accessible by the $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$ interface, you need utilities to convert traditional time series arrays into callable functions of time and state. The following example shows how to do this using the conversion function ts2func (time series to function).

Load a daily historical data set containing three-month Euribor rates and closing index levels of France's CAC 40 spanning the time interval February 7, 2001 to April 24, 2006.

load Data_GlobalIdx2

Simulate risk-neutral sample paths of the CAC 40 index using a geometric Brownian motion (GBM) model:

where r(t) represents evolution of the risk-free rate of return.

Furthermore, assume that you need to annualize the relevant information derived from the daily data (annualizing the data is optional, but is useful for comparison to other examples), and that each calendar year comprises 250 trading days.

dt      = 1/250;
returns = tick2ret(Dataset.CAC);
sigma   = std(returns)*sqrt(250);
yields  = Dataset.EB3M;
yields  = 360*log(1 + yields);

Compare the resulting sample paths obtained from two risk-neutral historical simulation approaches, where the daily Euribor yields serve as a proxy for the risk-free rate of return.

The first approach specifies the risk-neutral return as the sample average of Euribor yields, and therefore assumes a constant (non-dynamic) risk-free return.

nPeriods = length(yields);  % Simulated observations
rng(5713,'twister')
obj    = gbm(mean(yields),diag(sigma),'StartState',100)

obj = 
   Class GBM: Generalized Geometric Brownian Motion
   ------------------------------------------------
     Dimensions: State = 1, Brownian = 1
   ------------------------------------------------
      StartTime: 0
     StartState: 100
    Correlation: 1
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
         Return: 0.0278117
          Sigma: 0.231906

[X1,T] = simulate(obj,nPeriods,'DeltaTime',dt);

In contrast, the second approach specifies the risk-neutral return as the historical time series of Euribor yields. It therefore assumes a dynamic, yet deterministic, rate of return; this example does not illustrate stochastic interest rates. To illustrate this dynamic effect, use the ts2func utility.

r = ts2func(yields,'Times',(0:nPeriods - 1)');

ts2func packages a specified time series array inside a callable function of time and state, and synchronizes it with an optional time vector. For instance:

r(0,100)

ans = 
0.0470

evaluates the function at (t = 0, ${\mathit{X}}_{\mathit{t}}^{}$ = 100) and returns the first observed Euribor yield.

However, you can also evaluate the resulting function at any intermediate time t and state ${\mathit{X}}_{\mathit{t}}^{}$:

r(7.5,200)

ans = 
0.0472

Furthermore, the following command produces the same result when called with time alone:

r(7.5)

ans = 
0.0472

The equivalence of these last two commands highlights some important features.

When you specify parameters as functions, they must evaluate properly when passed a scalar, real-valued sample time (t), and an NVars-by-1 state vector (${\mathit{X}}_{\mathit{t}}^{}$). They must also generate an array of appropriate dimensions, which in the first case is a scalar constant, and in the second case is a scalar, piecewise constant function of time alone. You are not required to use either time (t) or state (${\mathit{X}}_{\mathit{t}}^{}$). In the current example, the function evaluates properly when passed time followed by state, thereby satisfying the minimal requirements. The fact that it also evaluates correctly when passed only time simply indicates that the function does not require the state vector ${\mathit{X}}_{\mathit{t}}^{}$. The important point to make is that it works when you pass it $\left(\mathit{t},{\mathit{X}}_{\mathit{t}}^{}\right)$.

Furthermore, the ts2func function performs a zero-order-hold (ZOH) piecewise constant interpolation. The notion of piecewise constant parameters is pervasive throughout the SDE architecture, and is discussed in more detail in Optimizing Accuracy: About Solution Precision and Error.

Complete the comparison by performing the second simulation using the same initial random number state. Simulate paths using a dynamic, deterministic rate of return.

rng(5713,'twister')
obj = gbm(r, diag(sigma),'StartState',100)

obj = 
   Class GBM: Generalized Geometric Brownian Motion
   ------------------------------------------------
     Dimensions: State = 1, Brownian = 1
   ------------------------------------------------
      StartTime: 0
     StartState: 100
    Correlation: 1
          Drift: drift rate function F(t,X(t)) 
      Diffusion: diffusion rate function G(t,X(t)) 
     Simulation: simulation method/function simByEuler
         Return: function ts2func/vector2Function
          Sigma: 0.231906

X2  = simulate(obj,nPeriods,'DeltaTime',dt);

Plot the series of risk-free reference rates to compare the two simulation trials.

subplot(2,1,1)
plot(dates,100*yields)
datetick('x')
xlabel('Date')
ylabel('Annualized Yield (%)')
title('Risk Free Rate(3-Mo Euribor Continuously-Compounded)')
subplot(2,1,2)
plot(T,X1,'red',T,X2,'blue')
xlabel('Time (Years)')
ylabel('Index Level')
title('Constant vs. Dynamic Rate of Return: CAC 40')
legend({'Constant Interest Rates' 'Dynamic Interest Rates'},...
    'Location', 'Best')

The paths are close but not exact. The blue line in the last plot uses all the historical Euribor data, and illustrates a single trial of a historical simulation.

Price European Stock Options Using Monte Carlo Simulation

As discussed in Simulate Positive Interest Rates, all simulation and interpolation methods allow you to specify one or more functions of the form:

to evaluate at the end of every sample time.

The related example illustrates a simple, common end-of-period processing function to ensure nonnegative interest rates. This example illustrates a processing function that allows you to avoid simulation outputs altogether.

Consider pricing European stock options by Monte Carlo simulation within a Black-Scholes-Merton framework. Assume that the stock has the following characteristics:

• The stock currently trades at 100.

• The stock pays no dividends.

• The stock's volatility is 50% per annum.

• The option strike price is 95.

• The option expires in three months.

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

To solve this problem, model the evolution of the underlying stock by a univariate geometric Brownian motion (GBM) model with constant parameters:

Furthermore, assume that the stock price is simulated daily, and that each calendar month comprises 21 trading days.

strike   = 95;
rate     = 0.1;
sigma    = 0.5;
dt       = 1/252;
nPeriods = 63;
T        = nPeriods*dt;
obj = gbm(rate,sigma,'StartState',100);

The goal is to simulate independent paths of daily stock prices, and calculate the price of European options as the risk-neutral sample average of the discounted terminal option payoff at expiration 63 days from now. This example calculates option prices by two approaches:

• A Monte Carlo simulation that explicitly requests the simulated stock paths as an output. The output paths are then used to price the options.

• An end-of-period processing function, accessible by time and state, that records the terminal stock price of each sample path. This processing function is implemented as a nested function with access to shared information.

Before simulation, invoke the example file to access the end-of-period processing function.

nTrials = 10000; % Number of independent trials (that is, paths)
f = Example_BlackScholes(nPeriods,nTrials)

f = struct with fields:
    BlackScholes: @Example_BlackScholes/saveTerminalStockPrice
       CallPrice: @Example_BlackScholes/getCallPrice
        PutPrice: @Example_BlackScholes/getPutPrice

Simulate 10000 independent trials (sample paths). Request the simulated stock price paths as an output, and specify an end-of-period processing function.

rng(88161,'twister')
X = simBySolution(obj,nPeriods,'DeltaTime',dt,... 
    'nTrials',nTrials,'Processes',f.BlackScholes);

Calculate the option prices directly from the simulated stock price paths. Because these are European options, ignore all intermediate stock prices.

call = mean(exp(-rate*T)*max(squeeze(X(end,:,:)) - strike, 0))

call = 
13.9342

put  = mean(exp(-rate*T)*max(strike - squeeze(X(end,:,:)), 0))

put = 
6.4166

Price the options indirectly by invoking the nested functions.

f.CallPrice(strike,rate)

ans = 
13.9342

f.PutPrice(strike,rate)

ans = 
6.4166

Although steps 3 and 4 produce the same option prices, the latter approach works directly with the terminal stock prices of each sample path. Therefore, it is much more memory efficient. In this example, there is no compelling reason to request an output.