# simByTransition

Simulate `CIR` sample paths with transition density

## Syntax

``[Paths,Times] = simByTransition(MDL,NPeriods)``
``[Paths,Times] = simByTransition(___,Name,Value)``

## Description

example

````[Paths,Times] = simByTransition(MDL,NPeriods)` simulates `NTrials` sample paths of `NVars` independent state variables driven by the Cox-Ingersoll-Ross (CIR) process sources of risk over `NPeriods` consecutive observation periods. `simByTransition` approximates a continuous-time CIR model using an approximation of the transition density function. ```

example

````[Paths,Times] = simByTransition(___,Name,Value)` specifies options using one or more name-value pair arguments in addition to the input arguments in the previous syntax.You can perform quasi-Monte Carlo simulations using the name-value arguments for `MonteCarloMethod` and `QuasiSequence`. For more information, see Quasi-Monte Carlo Simulation.```

## Examples

collapse all

Using the short rate, simulate the rate dynamics and term structures in the future using a CIR model. The CIR model is expressed as

`$dr\left(t\right)=\alpha \left(b-r\left(t\right)\right)dt+\sigma \sqrt{r\left(t\right)}dW\left(t\right)$`

The exponential affine form of the bond price is

`$B\left(t,T\right)={e}^{-A\left(t,T\right)r\left(t\right)+C\left(t,T\right)}$`

where

`$A\left(t,T\right)=\frac{2\left({e}^{\gamma \left(T-t\right)}-1\right)}{\left(\gamma +\alpha \right)\left({e}^{\gamma \left(T-t\right)}-1\right)+2\gamma }$`

`$B\left(t,T\right)=\frac{2\alpha b}{{\sigma }^{2}}\mathrm{log}\left(\frac{2\gamma {e}^{\left(\alpha +\gamma \right)\left(T-t\right)/2}}{\left(\gamma +\alpha \right)\left({e}^{\gamma \left(T-t\right)}-1\right)+2\gamma }\right)$`

and

`$\gamma =\sqrt{{\alpha }^{2}+2{\sigma }^{2}}$`

Define the parameters for the `cir` object.

```alpha = .1; b = .05; sigma = .05; r0 = .04;```

Define the function for bond prices.

```gamma = sqrt(alpha^2 + 2*sigma^2); A_func = @(t, T) ... 2*(exp(gamma*(T-t))-1)/((alpha+gamma)*(exp(gamma*(T-t))-1)+2*gamma); C_func = @(t, T) ... (2*alpha*b/sigma^2)*log(2*gamma*exp((alpha+gamma)*(T-t)/2)/((alpha+gamma)*(exp(gamma*(T-t))-1)+2*gamma)); P_func = @(t,T,r_t) exp(-A_func(t,T)*r_t+C_func(t,T));```

Create a `cir` object.

`obj = cir(alpha,b,sigma,'StartState',r0)`
```obj = Class CIR: Cox-Ingersoll-Ross ---------------------------------------- Dimensions: State = 1, Brownian = 1 ---------------------------------------- StartTime: 0 StartState: 0.04 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.05 Speed: 0.1 ```

Define the simulation parameters.

```nTrials = 100; nPeriods = 5; % Simulate future short over the next five years nSteps = 12; % Set intermediate steps to improve the accuracy```

Simulate the short rates. The returning path is a (NPeriods + 1)-by-`NVars`-by-`NTrials `three-dimensional time-series array. For this example, the size of the output is `6`-by-`1`-by-`100`.

```rng('default'); % Reproduce the same result rPaths = simByTransition(obj,nPeriods,'nTrials',nTrials,'nSteps',nSteps); size(rPaths)```
```ans = 1×3 6 1 100 ```
`rPathsExp = mean(rPaths,3);`

Determine the term structure over the next 30 years.

```maturity = 30; T = 1:maturity; futuresTimes = 1:nPeriods+1; % Preallocate simTermStruc simTermStructure = zeros(nPeriods+1,30); for i = futuresTimes for t = T bondPrice = P_func(i,i+t,rPathsExp(i)); simTermStructure(i,t) = -log(bondPrice)/t; end end plot(simTermStructure') legend('Current','1-year','2-year','3-year','4-year','5-year') title('Projected Term Structure for Next 5 Years') ylabel('Long Rate Maturity R(t,T)') xlabel('Time')```

The Cox-Ingersoll-Ross (CIR) short rate class derives directly from SDE with mean-reverting drift (`SDEMRD`): $d{X}_{t}=S\left(t\right)\left[L\left(t\right)-{X}_{t}\right]dt+D\left(t,{X}_{t}^{\frac{1}{2}}\right)V\left(t\right)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\left(0.1-{X}_{t}\right)dt+0.05{X}_{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 ```

Define the quasi-Monte Carlo simulation using the optional name-value arguments for `'MonteCarloMethod'` and `'QuasiSequence'`.

`[paths,time] = simByTransition(cir_obj,10,'ntrials',4096,'montecarlomethod','quasi','quasisequence','sobol');`

## Input Arguments

collapse all

Stochastic differential equation model, specified as a `cir` object. For more information on creating a `CIR` object, see `cir`.

Data Types: `object`

Number of simulation periods, specified as a positive scalar integer. The value of `NPeriods` determines the number of rows of the simulated output series.

Data Types: `double`

### Name-Value Arguments

Specify optional pairs of arguments as `Name1=Value1,...,NameN=ValueN`, where `Name` is the argument name and `Value` is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose `Name` in quotes.

Example: ```[Paths,Times] = simByTransition(CIR,NPeriods,'DeltaTimes',dt)```

Simulated trials (sample paths) of `NPeriods` observations each, specified as the comma-separated pair consisting of `'NTrials'` and a positive scalar integer.

Data Types: `double`

Positive time increments between observations, specified as the comma-separated pair consisting of `'DeltaTimes'` and a scalar or a `NPeriods`-by-`1` column vector.

`DeltaTime` represents the familiar dt found in stochastic differential equations, and determines the times at which the simulated paths of the output state variables are reported.

Data Types: `double`

Number of intermediate time steps within each time increment dt (defined as `DeltaTimes`), specified as the comma-separated pair consisting of `'NSteps'` and a positive scalar integer.

The `simByTransition` function partitions each time increment dt into `NSteps` subintervals of length dt/`NSteps`, and refines the simulation by evaluating the simulated state vector at `NSteps − 1` intermediate points. Although `simByTransition` does not report the output state vector at these intermediate points, the refinement improves accuracy by enabling the simulation to more closely approximate the underlying continuous-time process.

Data Types: `double`

Flag for storage and return method that indicates how the output array `Paths` is stored and returned, specified as the comma-separated pair consisting of `'StorePaths'` and a scalar logical flag with a value of `True` or `False`.

• If `StorePaths` is `True` (the default value) or is unspecified, then `simByTransition` returns `Paths` as a three-dimensional time series array.

• If `StorePaths` is `False` (logical `0`), then `simByTransition` returns the `Paths` output array as an empty matrix.

Data Types: `logical`

Monte Carlo method to simulate stochastic processes, specified as the comma-separated pair consisting of `'MonteCarloMethod'` and a string or character vector with one of the following values:

• `"standard"` — Monte Carlo using pseudo random numbers that has a convergence rate of O(N).

• `"quasi"` — Quasi-Monte Carlo rate of convergence is faster than standard Monte Carlo with errors approaching size of O(N-1).

• `"randomized-quasi"` — Quasi-random sequences, also called low-discrepancy sequences, are deterministic uniformly distributed sequences which are specifically designed to place sample points as uniformly as possible.

Data Types: `string` | `char`

Low discrepancy sequence to drive the stochastic processes, specified as the comma-separated pair consisting of `'QuasiSequence'` and a string or character vector with one of the following values:

• `"sobol"` — Quasi-random low-discrepancy sequences that use a base of two to form successively finer uniform partitions of the unit interval and then reorder the coordinates in each dimension

Note

If `MonteCarloMethod` option is not specified or specified as`"standard"`, `QuasiSequence` is ignored.

Data Types: `string` | `char`

Sequence of end-of-period processes or state vector adjustments, specified as the comma-separated pair consisting of `'Processes'` and a function or cell array of functions of the form

`${X}_{t}=P\left(t,{X}_{t}\right)$`

.

`simByTransition` applies processing functions at the end of each observation period. The processing functions accept the current observation time t and the current state vector Xt, and return a state vector that may adjust the input state.

If you specify more than one processing function, `simByTransition` invokes the functions in the order in which they appear in the cell array.

Data Types: `cell` | `function`

## Output Arguments

collapse all

Simulated paths of correlated state variables, returned as an ```(NPeriods + 1)```-by-`NVars`-by-`NTrials` three-dimensional time series array.

For a given trial, each row of `Paths` is the transpose of the state vector Xt at time t. When the input flag `StorePaths` = `False`, `simByTransition` returns `Paths` as an empty matrix.

Observation times associated with the simulated paths, returned as an `(NPeriods + 1)`-by-`1` column vector. Each element of `Times` is associated with the corresponding row of `Paths`.

collapse all

### Transition Density Simulation

The SDE has no solution such that r(t) = f(r(0),⋯).

In other words, the equation is not explicitly solvable. However, the transition density for the process is known.

The exact simulation for the distribution of r(t_1 ),⋯,r(t_n) is that of the process at times t_1,⋯,t_n for the same value of r(0). The transition density for this process is known and is expressed as

`$\begin{array}{l}r\left(t\right)=\frac{{\sigma }^{2}\left(1-{e}^{-\alpha \left(t-u\right)}}{4\alpha }{x}_{d}^{2}\left(\frac{4\alpha {e}^{-\alpha \left(t-u\right)}}{{\sigma }^{2}\left(1-{e}^{-\alpha \left(t-u\right)}\right)}r\left(u\right)\right),t>u\\ \text{where}\\ d\equiv \frac{4b\alpha }{{\sigma }^{2}}\end{array}$`

## Algorithms

Use the `simByTransition` function to simulate any vector-valued CIR process of the form

`$d{X}_{t}=S\left(t\right)\left[L\left(t\right)-{X}_{t}\right]dt+D\left(t,{X}_{t}^{\frac{1}{2}}\right)V\left(t\right)d{W}_{t}$`

where

• Xt is an `NVars`-by-`1` state vector of process variables.

• S is an `NVars`-by-`NVars` matrix of mean reversion speeds (the rate of mean reversion).

• L is an `NVars`-by-`1` vector of mean reversion levels (long-run mean or level).

• D is an `NVars`-by-`NVars` diagonal matrix, where each element along the main diagonal is the square root of the corresponding element of the state vector.

• V is an `NVars`-by-`NBrowns` instantaneous volatility rate matrix.

• dWt is an `NBrowns`-by-`1` Brownian motion vector.

## References

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

## Version History

Introduced in R2018b