| GARCH Toolbox™ | |
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:
 | (5-13) |
where:
μ is a diagonal matrix of index returns.
D is a diagonal matrix with Xt along the diagonal.
σ is a diagonal matrix of index standard deviation.
Create an SDE object to represent the equity market model.
load SDE_Data
SDE_Data =
Dates: [1359x1 double]
Canada: [1359x1 double]
France: [1359x1 double]
Germany: [1359x1 double]
Japan: [1359x1 double]
UK: [1359x1 double]
US: [1359x1 double]
Euribor3M: [1359x1 double]
prices = [SDE_Data.Canada SDE_Data.France SDE_Data.Germany ...
SDE_Data.Japan SDE_Data.UK SDE_Data.US];
Convert daily prices to returns:
returns = price2ret(prices);
Compute data statistics to input to simulation methods:
nVariables = size(returns, 2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns); covariance = cov(returns); t = 0; X = 100; X = X(ones(nVariables,1));
Create simple anonymous drift and diffusion functions accessible by (t, Xt):
F = @(t,X) diag(expReturn) * X; G = @(t,X) diag(X) * diag(sigma);
Use these functions to create an SDE object to represent the market model in Equation 5-13:
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 constructor requires additional information to determine the dimensionality of the model, because the functions passed to the SDE constructor are known only by their (t, Xt) interface. In other words, the SDE constructor 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 (t, Xt).
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. Note that the SDE engine has assigned the default simulation method, simByEuler, to the Simulation parameter.
Create an SDEDDO object to represent the market model in Equation 5-13:
Create drift and diffusion objects:
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 the SDEDDO constructor:
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 constructor requires two input objects that provide more information than the two functions from step 4 of Implementation 1: Using SDE Objects. 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 constructor, 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, F(t,Xt) and G(t,Xt), respectively. The additional parameters, A, B, Alpha, and Sigma, are arrays of appropriate dimension, indicating static (non-time-varying) parameters. In other words, A(t,Xt), B(t,Xt), Alpha(t,Xt), and Sigma(t,Xt) are constant functions of time and state.
Create SDELD, CEV, and GBM objects to represent the market model in Equation 5-13.
Create an SDELD object:
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 arrayCreate a CEV object:
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
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 class as Return. This is an intuitive name commonly associated with equity models.
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
randn('state', 100)
[S,T] = SDE.simulate(nPeriods, 'DeltaTime', dt);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 Xt at time t.
whos S
Name Size Bytes Class Attributes
S 250x6 12000 double
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')
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 Implementation 4: Using the Default Simulate Method:
randn('state', 100)
[S,T] = SDE.simByEuler(nPeriods, 'DeltaTime', dt);Simulate 10 trials with the same initial conditions:
randn('state', 100)
[S,T] = SDE.simulate(nPeriods, 'DeltaTime', dt, 'nTrials', 10);Now the output array S is an (NPERIODS + 1)-by-nVariables-by-nTrials time-series array:
whos S Name Size Bytes Class Attributes S 250x6x10 120000 double
whose first realization is identical to the single paths just plotted:
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')
Finally, consider simulation using GBM simulation methods. Separable GBM models have two specific simulation methods:
An overloaded Euler simulation method, designed for optimal performance
A method that provides an approximate solution of the underlying stochastic differential equation, designed for accuracy
To illustrate the performance benefit of the overloaded Euler approximation method, increase the number of trials to 10000:
randn('state', 100)
[X,T] = GBM.simulate(nPeriods, 'DeltaTime', dt, ...
'nTrials', 10000);The output X is a much larger time-series array:
whos X Name Size Bytes Class Attributes X 250x6x10000 120000000 double
With this sample size, you can examine the terminal distribution of, for example, Canada's TSX Composite to verify qualitatively the log-normal character of the data:
hist(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:
randn('state', 100)
[X,T] = GBM.simBySolution(nPeriods,...
'DeltaTime', dt, 'nTrials', 10);
subplot(2,1,1)
plot(T, S(:,:,1)), xlabel('Trading Day'),...
ylabel('Price')
title('First Path of Multi-Dimensional Market Model: ...
Euler Approximation')
subplot(2,1,2)
plot(T, X(:,:,1)), xlabel('Trading Day'),...
ylabel('Price')
title('First Path of Multi-Dimensional 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.
The following plot illustrates 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 Approximation 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 can be 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 an 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 piece-wise constant parameters.
When all model parameters are piece-wise 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 of Solutions.
All simulation methods require that you specify a time grid by specifying the number of periods (NPERIODS). You can also optionally specify a scalar or vector of strictly positive time increments (DeltaTime) and intermediate time steps (NSTEPS). These parameters, along with an initial sample time associated with the object (StartTime), uniquely determine the sequence of times at which the state vector is sampled. Thus, simulation methods allow you to traverse the time grid from beginning to end (that is, from left to right).
In contrast, interpolation methods allow you to traverse the time grid in any order, allowing both forward and backward movements in time. They allow you to specify a vector of interpolation times whose elements do not have to be unique.
Many references define the Brownian Bridge as a conditional simulation combined with a scheme for traversing the time grid, effectively merging two distinct algorithms. In contrast, the interpolation method offered here provides additional flexibility by intentionally separating the algorithms. In this method for moving about a time grid, you perform an initial Monte Carlo simulation to sample the state at the terminal time, and then successively sample intermediate states by stochastic interpolation. The first few samples determine the overall behavior of the paths, while later samples progressively refine the structure. Such algorithms are often called variance reduction techniques. This algorithm is particularly simple when the number of interpolation times is a power of 2. In this case, each interpolation falls midway between two known states, refining the interpolation using a method like bisection. This example highlights the flexibility of refined interpolation by implementing this power-of-two algorithm.
Load a historical data set of three-month Euribor rates (http://www.euribor.org/default.htm), observed daily, and corresponding trading dates spanning the time interval from February 7, 2001 through April 24, 2006:
clear, clf
load SDE_Data
plot(SDE_Data.Dates, 100 * SDE_Data.Euribor3M)
datetick('x'), xlabel('Date'), ylabel('Daily Yield (%)')
title('3-Month Euribor as a Daily Effective Yield')
Now fit a simple univariate Vasicek model to the daily equivalent yields of the three-month Euribor data:
Given initial conditions, the distribution of the short rate at some time T in the future is Gaussian with mean:
and variance:
To calibrate this simple short rate model, rewrite it in more familiar regression format:
where:
perform an ordinary linear regression where the model volatility is proportional to the standard error of the residuals:
yields = SDE_Data.Euribor3M; regressors = [ones(length(yields) - 1, 1) yields(1:end-1)]; [coefficients, intervals, residuals] = ... regress(diff(yields), regressors); dt = 1; % time increment = 1 day speed = -coefficients(2)/dt; level = -coefficients(1)/coefficients(2); sigma = std(residuals)/sqrt(dt);
Create an HWV object with an initial StartState set to the most recently observed short rate:
obj = hwv(speed, level, sigma, 'StartState', yields(end))
obj =
Class HWV: Hull-White/Vasicek
----------------------------------------
Dimensions: State = 1, Brownian = 1
----------------------------------------
StartTime: 0
StartState: 7.70408e-005
Correlation: 1
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
Sigma: 4.77637e-007
Level: 6.00424e-005
Speed: 0.00228854
Assume, for example, that you simulate the fitted model over 64 (26) trading days, using a refined Brownian bridge with the power-of-two algorithm instead of the usual beginning-to-end Monte Carlo simulation approach. Furthermore, assume that the initial time and state coincide with those of the last available observation of the historical data, and that the terminal state is the expected value of the Vasicek model 64 days into the future. In this case, you can assess the behavior of various paths that all share the same initial and terminal states, perhaps to support pricing path-dependent interest rate options over a three-month interval.
Create a vector of interpolation times to traverse the time grid by moving both forward and backward in time. Specifically, the first interpolation time is set to the most recent short rate observation time, the second interpolation time is set to the terminal time, and subsequent interpolation times successively sample intermediate states:
T = 64;
times = (1:T)';
t = NaN(length(times) + 1, 1);
t(1) = obj.StartTime;
t(2) = T;
delta = T;
jMax = 1;
iCount = 3;
for k = 1:log2(T)
i = delta / 2;
for j = 1:jMax
t(iCount) = times(i);
i = i + delta;
iCount = iCount + 1;
end
jMax = 2 * jMax;
delta = delta / 2;
endExamine the sequence of interpolation times generated by this algorithm:
stem(1:length(t), t, 'filled')
xlabel('Index'), ylabel('Interpolation Time (Days)')
title ('Sampling Scheme for the Power-of-Two Algorithm')
The first few samples are widely separated in time and determine the course structure of the paths. Later samples are closely spaced and progressively refine the detailed structure.
Now that you have generated the sequence of interpolation times, initialize a course time-series grid to begin the interpolation. The sampling process begins at the last observed time and state taken from the historical short rate series, and ends 64 days into the future at the expected value of the Vasicek model derived from the calibrated parameters:
average = obj.StartState * exp(-speed * T) + ... level * (1 - exp(-speed * T)); X = [obj.StartState ; average];
Generate 5 sample paths, setting the Refine input flag to TRUE to insert each new interpolated state into the time series grid as it becomes available. Perform interpolation on a trial-by-trial basis. Because the input time series X has five trials (where each page of the 3-dimensional time series represents an independent trial), the interpolated output series Y also has five pages:
nTrials = 5;
randn('state', 0)
Y = obj.interpolate(t, X(:,:,ones(nTrials,1)), ...
'Times',[obj.StartTime T], 'Refine', true);Plot the resulting sample paths. Because the interpolation times do not monotonically increase, sort the times and reorder the corresponding short rates:
[t,i] = sort(t);
Y = squeeze(Y);
Y = Y(i,:);
plot(t, 100 * Y), hold('on')
plot(t([1 end]), 100 * Y([1 end],1),'. black', ...
'MarkerSize', 20)
xlabel('Interpolation Time (Days into the Future)')
ylabel('Yield (%)'), hold('off')
title ('Euribor Yields Obtained by Brownian Bridge ...
Interpolation')
The short rates in this plot represent alternative sample paths that share the same initial and terminal values. They illustrate a special, though simplistic, case of a broader sampling technique known as stratified sampling. For a more sophisticated example of stratified sampling, see User-Specified Random Number Generation: Stratified Sampling.
Although this simple example simulated a univariate Vasicek interest rate model, it applies to problems of any dimensionality.
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 are 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 Implementing Multidimensional Equity Market Models). The market model to simulate is:
where μ is a diagonal matrix.
Load the SDE_Data data set:
load SDE_Data, SDE_Data
SDE_Data =
Dates: [1359x1 double]
Canada: [1359x1 double]
France: [1359x1 double]
Germany: [1359x1 double]
Japan: [1359x1 double]
UK: [1359x1 double]
US: [1359x1 double]
Euribor3M: [1359x1 double]
prices = [SDE_Data.Canada SDE_Data.France...
SDE_Data.Germany SDE_Data.Japan...
SDE_Data.UK SDE_Data.US];
Convert the daily prices to returns:
returns = price2ret(prices);
Specify Sigma and Correlation using the first technique:
Using the first technique, 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:
Using the second technique, 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:
randn('state', 0)
[X1,T] = GBM1.simByEuler(1000); % correlated Brownian motion
randn('state', 0)
[X2,T] = GBM2.simByEuler(1000); % standard Brownian motionWhen based on the same initial random number state, each technique generates identical asset price paths:
subplot(2,1,1), plot(T, X1)
title('Correlated Sample Paths ...
from Correlated Brownian Motion')
ylabel('Asset Price')
subplot(2,1,2), plot(T, X2)
title('Correlated Sample Paths ...
from Standard Brownian Motion')
xlabel('Trading Day'), ylabel('Asset Price')
As previously discussed, 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 generally 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 (t, Xt) interface. This reflects the fact that GBM models (and others) 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 (t, Xt) 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 (t, Xt) 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 (t, Xt) 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 of 3-month Euribor rates and closing index levels of France's CAC 40 spanning the time interval February 7, 2001 to April 24, 2006:
clear load SDE_Data
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 = price2ret(SDE_Data.France); sigma = std(returns) * sqrt(250); yields = SDE_Data.Euribor3M; 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
randn('state', 25)
obj = gbm(mean(yields), diag(sigma), 'StartState', 100)
[X1,T] = obj.simulate(nPeriods, 'DeltaTime', dt);
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.231875
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, X t = 100) and returns the first observed Euribor yield. However, you can also evaluate the resulting function at any intermediate time t and state Xt:
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 a NVARS-by-1 state vector (Xt). 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, piece-wise constant function of time alone.
You are not required to use either time (t) or state (Xt). 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 Xt. The important point to make is that it works when you pass it (t, Xt).
Furthermore, the ts2func function performs a zero-order-hold (ZOH) piece-wise constant interpolation. The notion of piece-wise constant parameters is pervasive throughout the SDE architecture, and is discussed in more detail in Optimizing Accuracy of Solutions.
Complete the comparison by performing the second simulation using the same initial random number state:
randn('state', 25)
obj = gbm(r, diag(sigma), 'StartState', 100)
X2 = obj.simulate(nPeriods, 'DeltaTime', dt);
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.231875
Plot the series of risk-free reference rates to compare the two simulation trials:
subplot(2,1,1)
plot(SDE_Data.Dates, 100 * yields)
datetick('x'), xlabel('Date'), ...
ylabel('Annualized Yield (%)')
title('Risk Free Rate ...
(3-Month 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.
All simulation and interpolation methods allow you to specify a sequence of functions, or background processes, to evaluate at the end of every sample time period. This period includes any intermediate time steps determined by the optional NSTEPS input, as discussed in Optimizing Accuracy of Solutions. These functions are specified as callable functions of time and state, and must return an updated state vector Xt:
You must specify multiple processing functions as a cell array of functions. These functions are invoked in the order in which they appear in the cell array.
Processing functions are not required to use time (t) or state (Xt). They are also not required to update or change the input state vector. In fact, simulation and interpolation methods have no knowledge of any implementation details, and in this respect, they only adhere to a published interface.
Such processing functions provide a powerful modeling tool that can solve a variety of problems. Such functions allow you to, for example, specify boundary conditions, accumulate statistics, plot graphs, and price path-dependent options.
Except for Brownian motion (BM) models, the individual components of the simulated state vector typically represent variables whose real-world counterparts are inherently positive quantities, such as asset prices or interest rates. However, by default, most of the simulation and interpolation methods provided here model the transition between successive sample times as a scaled (possibly multivariate) Gaussian draw. Consequently, when approximating a continuous-time process in discrete time, the state vector may not remain positive. The only exception is the simBySolution logarithmic transform of separable geometric Brownian motion models. Moreover, by default, none of the simulation and interpolation methods make adjustments to the state vector. Therefore, you are responsible for ensuring that all components of the state vector remain positive as appropriate.
Fortunately, specifying non-negative states ensures a simple end-of-period processing adjustment. Although this adjustment is widely applicable, it is revealing when applied to a univariate CIR square-root diffusion model:
Perhaps the primary appeal of univariate CIR models where:
is that the short rate remains positive. However, the positivity of short rates only holds for the underlying continuous-time model.
To illustrate the latter statement, simulate daily short rates of the CIR model over one calendar year (approximately 250 trading days):
randn('state', 10)
obj = cir(0.25, @(t,X) 0.1, 0.2, 'StartState', 0.02);
[X,T] = obj.simByEuler(250, 'DeltaTime', ...
1/250, 'nTrials', 5);Interest rates can become negative if the resulting paths are simulated in discrete time. Moreover, since CIR models incorporate a square root diffusion term, the short rates might even become complex:
[T(200:210) X(200:210,1,5)] ans = 0.7960 0.0023 0.8000 0.0023 0.8040 0.0022 0.8080 0.0010 0.8120 0.0005 0.8160 0.0003 0.8200 -0.0001 0.8240 -0.0000 - 0.0002i 0.8280 0.0002 - 0.0003i 0.8320 0.0005 - 0.0004i 0.8360 0.0007 - 0.0004i
Repeat the simulation, this time specifying a processing function that takes the absolute magnitude of the short rate at the end of each period. You can access the processing function by time and state (t, Xt), but it only uses the state vector of short rates Xt:
randn('state', 10)
[Y,T] = obj.simByEuler(250, 'DeltaTime', 1/250, ...
'nTrials', 5, 'Processes', @(t,X) abs(X));Graphically compare the magnitude of the unadjusted path (with negative and complex numbers!) to the corresponding path kept positive by using an end-of-period processing function over the time span of interest:
clf
plot(T, 100 * abs(X(:,1,5)), 'red', T, ...
100 * Y(:,1,5), 'blue')
axis([0.75 1 0 1.2])
xlabel('Time (Years)'), ylabel('Short Rate (%)')
title('Univariate CIR Short Rates')
legend({'Negative/Complex Rates' 'Positive Rates'}, ...
'Location', 'Best')
As discussed in Ensuring Positive State Variables, 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 previous example illustrated a simple, common end-of-period processing function to ensure non-negative 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. For more information, see the demo blackScholesExample.m.
Before simulation, invoke the example file to access the end-of-period processing function:
nTrials = 10000; % Number of independent trials (i.e., paths)
f = blackScholesExample(nPeriods, nTrials)
f =
BlackScholes: @blackScholesExample/saveTerminalStockPrice
CallPrice: @blackScholesExample/getCallPrice
PutPrice: @blackScholesExample/getPutPrice
Simulate 10000 independent trials (sample paths). Request the simulated stock price paths as an output, and specify an end-of-period processing function:
randn('state', 0)
X = obj.simBySolution(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))
put = mean(exp(-rate * T) * max(strike - ...
squeeze(X(end,:,:)), 0))
call =
13.9964
put =
6.2028Price the options indirectly by invoking the nested functions:
f.CallPrice(strike, rate)
f.PutPrice (strike, rate)
ans =
13.9964
ans =
6.2028
For reference, the theoretical call and put prices computed from the Black-Scholes option formulas are 13.6953 and 6.3497, respectively.
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.
Simulation methods allow you to specify a noise process directly, as a callable function of time and state:
Stratified sampling is a variance reduction technique that constrains a proportion of sample paths to specific subsets (or strata) of the sample space.
This example specifies a noise function to stratify the terminal value of a univariate equity price series. Starting from known initial conditions, the function first stratifies the terminal value of a standard Brownian motion, and then samples the process from beginning to end by drawing conditional Gaussian samples using a Brownian bridge.
The stratification process assumes that each path is associated with a single stratified terminal value such that the number of paths is equal to the number of strata. This technique is called proportional sampling. This example is similar to, yet more sophisticated than, the one discussed in Stochastic Interpolation and the Brownian Bridge. Since stratified sampling requires knowledge of the future, it also requires more sophisticated time synchronization; specifically, the function in this example requires knowledge of the entire sequence of sample times. For more information, see the demo stratifiedExample.m.
The function implements proportional sampling by partitioning the unit interval into bins of equal probability by first drawing a random number uniformly distributed in each bin. The inverse cumulative distribution function of a standard N(0,1) Gaussian distribution then transforms these stratified uniform draws. Finally, the resulting stratified Gaussian draws are scaled by the square root of the terminal time to stratify the terminal value of the Brownian motion.
The noise function does not return the actual Brownian paths, but rather the Gaussian draws Z(t,Xt) that drive it.
This example first stratifies the terminal value of a univariate, zero-drift, unit-variance-rate Brownian motion (BM) model:
Assume that 10 paths of the process are simulated daily over a three-month period. Also assume that each calendar month and year consist of 21 and 252 trading days, respectively:
randn('state', 10), rand('twister', 0)
dt = 1 / 252; % 1 day = 1/252 years
nPeriods = 63; % 3 months = 63 trading days
T = nPeriods * dt; % 3 months = 0.25 years
nPaths = 10; % # of simulated paths
obj = bm(0, 1, 'StartState', 0);
sampleTimes = cumsum([obj.StartTime ...
dt(ones(nPeriods,1))]);
z = stratifiedExample(nPaths, sampleTimes)
z = @stratifiedExample/stratifiedSampling
Simulate the standard Brownian paths by explicitly passing the stratified sampling function to the simulation method:
X = obj.simulate(nPeriods, 'DeltaTime', dt, ... 'nTrials', nPaths, 'Z', z);
For convenience, reorder the output sample paths by reordering the 3-dimensional output to a 2-dimensional equivalent array:
X = squeeze(X);
Verify the stratification:
Recreate the uniform draws with proportional sampling:
rand('twister', 0)
U = ((1:nPaths)' - 1 + rand(nPaths,1))/nPaths; Transform them to obtain the terminal values of standard Brownian motion:
WT = norminv(U) * sqrt(T); % Stratified Brownian motion.
Plot the terminal values and output paths on the same figure:
plot(sampleTimes, X), hold('on')
xlabel('Time (Years)'), ylabel('Brownian State')
title('Terminal Stratification: Standard Brownian Motion')
plot(T, WT, '. black', T, WT, 'o black')
hold('off')
The last value of each sample path (the last row of the output array X) coincides with the corresponding element of the stratified terminal value of the Brownian motion. This occurs because the simulated model and the noise generation function both represent the same standard Brownian motion.
However, you can use the same stratified sampling function to stratify the terminal price of constant-parameter geometric Brownian motion models. In fact, you can use the stratified sampling function to stratify the terminal value of any constant-parameter model driven by Brownian motion if the model's terminal value is a monotonic transformation of the terminal value of the Brownian motion.
To illustrate this, load the SDE_data data set and simulate risk-neutral sample paths of the FTSE 100 index using a geometric Brownian motion (GBM) model with constant parameters:
where the average Euribor yield represents the risk-free rate of return.
Assume that the relevant information derived from the daily data is annualized, and that each calendar year comprises 252 trading days:
returns = price2ret(SDE_Data.UK); sigma = std(returns) * sqrt(252); rate = SDE_Data.Euribor3M; rate = mean(360 * log(1 + rate));
Create the GBM model, assuming the FTSE 100 starts at 100:
obj = gbm(rate, sigma, 'StartState', 100);
Determine the sample time and simulate the price paths.
In what follows, NSTEPS specifies the number of intermediate time steps within each time increment DeltaTime. Each increment DeltaTime is partitioned into NSTEPS subintervals of length DeltaTime/nSteps each, refining the simulation by evaluating the simulated state vector at NSTEPS–1 intermediate points. This refinement improves accuracy by allowing the simulation to more closely approximate the underlying continuous-time process without storing the intermediate information:
nSteps = 1;
sampleTimes = cumsum([obj.StartTime ;
dt(ones(nPeriods * nSteps,1))/nSteps]);
z = stratifiedExample(nPaths, sampleTimes);
randn('state', 10), rand('twister', 0)
[Y, Times] = obj.simBySolution(nPeriods, 'nTrials, nPaths, ...
'DeltaTime', dt, 'nSteps', nSteps, 'Z', z);
Y = squeeze(Y); % Reorder to a 2-D array
plot(Times, Y)
xlabel('Time (Years)'), ylabel('Index Level')
title('FTSE 100 Terminal Stratification: ...
Geometric Brownian Motion')
Although the terminal value of the Brownian motion shown in the latter plot is normally distributed, and the terminal price in the previous plot is lognormally distributed, the corresponding paths of each graph are similar.
| Using SDE Objects to Create Models | Creating User-Specified Functions | |
| © 1984-2008- The MathWorks, Inc. - Site Help - Patents - Trademarks - Privacy Policy - Preventing Piracy - RSS |