# simBySolution

Simulate approximate solution of diagonal-drift Merton jump diffusion process

## Description

example

[Paths,Times,Z,N] = simBySolution(MDL,NPeriods) simulates NNTrials sample paths of NVars correlated state variables driven by NBrowns Brownian motion sources of risk and NJumps compound Poisson processes representing the arrivals of important events over NPeriods consecutive observation periods. The simulation approximates continuous-time Merton jump diffusion process by an approximation of the closed-form solution.

example

[Paths,Times,Z,N] = simBySolution(___,Name,Value) specifies options using one or more name-value pair arguments in addition to the input arguments in the previous syntax.

## Examples

collapse all

Simulate the approximate solution of diagonal-drift Merton process.

Create a merton object.

AssetPrice = 80;
Return = 0.03;
Sigma = 0.16;
JumpMean = 0.02;
JumpVol = 0.08;
JumpFreq = 2;

mertonObj = merton(Return,Sigma,JumpFreq,JumpMean,JumpVol,...
'startstat',AssetPrice)
mertonObj =
Class MERTON: Merton Jump Diffusion
----------------------------------------
Dimensions: State = 1, Brownian = 1
----------------------------------------
StartTime: 0
StartState: 80
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.16
Return: 0.03
JumpFreq: 2
JumpMean: 0.02
JumpVol: 0.08

Use simBySolution to simulate NTrials sample paths of NVARS correlated state variables driven by NBrowns Brownian motion sources of risk and NJumps compound Poisson processes representing the arrivals of important events over NPeriods consecutive observation periods. The function approximates continuous-time Merton jump diffusion process by an approximation of the closed-form solution.

nPeriods = 100;
[Paths,Times,Z,N] = simBySolution(mertonObj, nPeriods,'nTrials', 3)
Paths =
Paths(:,:,1) =

1.0e+03 *

0.0800
0.0662
0.1257
0.1863
0.2042
0.2210
0.2405
0.3143
0.4980
0.4753
0.4088
0.5627
0.6849
0.6662
0.7172
0.7710
0.6758
0.5528
0.4777
0.6314
0.7290
0.7265
0.6018
0.6630
0.5531
0.5919
0.5580
0.7209
0.8122
0.6494
0.8194
0.7434
0.6887
0.6873
0.7052
0.8532
0.5498
0.4686
0.5445
0.4291
0.5118
0.4138
0.4986
0.4331
0.4687
0.5235
0.4944
0.4616
0.3621
0.4860
0.4461
0.4268
0.4179
0.3913
0.5225
0.4346
0.3433
0.3635
0.3604
0.3736
0.3771
0.4883
0.4785
0.4859
0.5719
0.6593
0.7232
0.8269
0.7894
0.8895
0.9131
0.7396
0.9902
1.4258
1.1410
1.1657
1.2759
1.2797
1.2587
1.5073
1.5914
1.2676
1.5111
1.4698
1.5310
1.0471
1.3415
1.2142
1.3649
1.9905
1.9329
1.5042
1.7000
2.2315
2.6107
2.2992
2.6765
2.7024
1.6837
1.0520
1.1556

Paths(:,:,2) =

80.0000
67.0894
98.3231
108.1133
102.2668
116.5130
92.6337
94.7715
110.7864
125.7798
120.6730
116.9214
106.8356
118.3119
190.3385
228.3806
271.8072
272.0175
306.3696
249.6461
427.2599
310.1494
471.7915
370.6712
426.4875
393.6037
423.9768
436.6450
423.3666
415.2689
578.7237
448.8291
358.5539
314.4588
284.7537
345.2281
379.3241
432.3968
284.6978
428.3203
314.5781
326.2297
236.1605
178.9878
175.8927
177.5584
140.5670
124.3399
111.5921
114.6988
101.7877
72.8823
61.0876
54.7438
53.9104
44.3239
32.8282
35.8978
44.7213
37.6385
34.8707
33.4812
35.0828
37.3844
50.3077
49.7005
41.2006
58.0578
51.8254
42.3636
38.3241
40.1687
35.9465
44.4746
36.3203
31.4723
25.3097
23.4042
14.5024
11.9513
11.7996
13.2874
14.9033
14.9986
14.9639
18.8188
16.5700
17.8684
13.5567
13.5978
11.3215
10.6453
9.9437
10.9639
14.0077
16.5691
12.1943
10.7238
11.5439
9.3313
10.3501

Paths(:,:,3) =

80.0000
79.6896
69.0705
57.4353
54.6468
61.1361
78.0797
104.5536
107.1168
87.1463
54.5801
59.8430
67.0858
74.7163
65.0742
90.0205
70.0329
94.1883
88.2437
100.7302
127.2244
111.4070
81.0410
93.1479
72.5876
74.3940
71.8182
78.4764
90.1952
89.6539
70.3198
50.4493
58.2573
52.1928
67.7723
81.1286
112.6400
108.8060
103.0418
104.3689
120.8792
89.2307
66.3967
76.2541
57.1963
56.8041
40.4475
34.5959
45.2467
44.6159
52.2680
63.3114
69.8554
102.0669
76.8265
84.8615
62.4934
70.3915
54.4665
60.1859
68.3690
73.3205
87.8904
82.7777
94.8760
88.8936
103.9546
103.4198
99.0468
135.2132
117.9348
120.8927
126.9568
120.5084
119.4830
154.8170
165.2276
180.3558
150.8172
155.2828
138.6475
179.8007
158.8069
166.0540
229.0607
253.4962
240.1957
192.3787
225.7069
311.1060
353.6839
463.5303
515.0606
569.4017
488.1785
331.1247
392.7017
379.5234
238.3932
186.9090
209.5703

Times = 101×1

0
1
2
3
4
5
6
7
8
9
⋮

Z =
Z(:,:,1) =

-1.3077
3.5784
3.0349
0.7147
1.4897
0.6715
1.6302
0.7269
-0.7873
-1.0689
1.4384
1.3703
-0.2414
-0.8649
0.6277
-0.8637
-1.1135
-0.7697
1.1174
0.5525
0.0859
-1.0616
0.7481
-0.7648
0.4882
1.4193
1.5877
0.8351
-1.1658
0.7223
0.1873
-0.4390
-0.8880
0.3035
0.7394
-2.1384
-1.0722
1.4367
-1.2078
1.3790
-0.2725
0.7015
-0.8236
0.2820
1.1275
0.0229
-0.2857
-1.1564
0.9642
-0.0348
-0.1332
-0.2248
-0.8479
1.6555
-0.8655
-1.3320
0.3335
-0.1303
0.8620
-0.8487
1.0391
0.6601
-0.2176
0.0513
0.4669
0.1832
0.3071
0.2614
-0.1461
-0.8757
-1.1742
1.5301
1.6035
-1.5062
0.2761
0.3919
-0.7411
0.0125
1.2424
0.3503
-1.5651
0.0983
-0.0308
-0.3728
-2.2584
1.0001
-0.2781
0.4716
0.6524
1.0061
-0.9444
0.0000
0.5946
0.9298
-0.6516
-0.0245
0.8617
-2.4863
-2.3193
0.4115

Z(:,:,2) =

-0.4336
2.7694
0.7254
-0.2050
1.4090
-1.2075
0.4889
-0.3034
0.8884
-0.8095
0.3252
-1.7115
0.3192
-0.0301
1.0933
0.0774
-0.0068
0.3714
-1.0891
1.1006
-1.4916
2.3505
-0.1924
-1.4023
-0.1774
0.2916
-0.8045
-0.2437
-1.1480
2.5855
-0.0825
-1.7947
0.1001
-0.6003
1.7119
-0.8396
0.9610
-1.9609
2.9080
-1.0582
1.0984
-2.0518
-1.5771
0.0335
0.3502
-0.2620
-0.8314
-0.5336
0.5201
-0.7982
-0.7145
-0.5890
-1.1201
0.3075
-0.1765
-2.3299
0.3914
0.1837
-1.3617
-0.3349
-1.1176
-0.0679
-0.3031
0.8261
-0.2097
-1.0298
0.1352
-0.9415
-0.5320
-0.4838
-0.1922
-0.2490
1.2347
-0.4446
-0.2612
-1.2507
-0.5078
-3.0292
-1.0667
-0.0290
-0.0845
0.0414
0.2323
-0.2365
2.2294
-1.6642
0.4227
-1.2128
0.3271
-0.6509
-1.3218
-0.0549
0.3502
0.2398
1.1921
-1.9488
0.0012
0.5812
0.0799
0.6770

Z(:,:,3) =

0.3426
-1.3499
-0.0631
-0.1241
1.4172
0.7172
1.0347
0.2939
-1.1471
-2.9443
-0.7549
-0.1022
0.3129
-0.1649
1.1093
-1.2141
1.5326
-0.2256
0.0326
1.5442
-0.7423
-0.6156
0.8886
-1.4224
-0.1961
0.1978
0.6966
0.2157
0.1049
-0.6669
-1.9330
0.8404
-0.5445
0.4900
-0.1941
1.3546
0.1240
-0.1977
0.8252
-0.4686
-0.2779
-0.3538
0.5080
-1.3337
-0.2991
-1.7502
-0.9792
-2.0026
-0.0200
1.0187
1.3514
-0.2938
2.5260
-1.2571
0.7914
-1.4491
0.4517
-0.4762
0.4550
0.5528
1.2607
-0.1952
0.0230
1.5270
0.6252
0.9492
0.5152
-0.1623
1.6821
-0.7120
-0.2741
-1.0642
-0.2296
-0.1559
0.4434
-0.9480
-0.3206
-0.4570
0.9337
0.1825
1.6039
-0.7342
0.4264
2.0237
0.3376
-0.5900
-1.6702
0.0662
1.0826
0.2571
0.9248
0.9111
1.2503
-0.6904
-1.6118
1.0205
-0.0708
-2.1924
-0.9485
0.8577

N =
N(:,:,1) =

1
2
1
0
2
0
1
3
4
2
1
0
1
1
1
1
0
0
3
2
2
1
0
1
1
3
3
4
2
4
1
1
2
0
2
2
3
2
1
3
2
2
1
1
1
3
0
2
2
1
0
1
1
1
1
0
2
2
1
1
6
7
3
2
2
1
3
3
4
3
0
1
7
2
0
5
2
2
1
2
1
3
0
2
5
2
2
4
2
3
1
2
6
1
0
3
3
1
1
3

N(:,:,2) =

2
2
2
0
4
1
2
3
1
2
1
4
2
4
2
2
2
2
1
5
3
1
3
3
1
3
5
1
4
2
2
1
2
1
1
6
0
2
2
3
2
2
1
0
1
5
5
0
1
1
2
1
2
3
2
2
1
2
2
0
3
1
5
3
3
0
2
1
2
0
4
1
3
1
2
2
2
1
0
2
2
2
2
1
1
3
1
2
2
1
4
1
3
3
0
1
1
1
2
3

N(:,:,3) =

3
2
2
1
4
2
3
0
0
4
3
2
3
1
1
1
1
3
4
1
2
3
1
1
1
1
0
3
0
1
0
5
0
2
4
3
1
0
1
4
3
3
2
1
2
3
1
4
4
1
1
2
2
1
1
1
2
1
6
1
2
1
3
2
2
1
3
1
7
0
1
5
1
1
3
4
3
1
2
2
1
2
1
1
1
1
1
2
3
4
2
1
3
2
1
1
0
1
3
0

## Input Arguments

collapse all

Merton model, specified as a merton object. You can create a merton object using merton.

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 comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: [Paths,Times,Z,N] = simBySolution(merton,NPeriods,'DeltaTimes',dt,'NNTrials',10)

Simulated NTrials (sample paths) of NPeriods observations each, specified as the comma-separated pair consisting of 'NNTrials' 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 an NPeriods-by-1 column vector.

DeltaTimes 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 (specified as DeltaTimes), specified as the comma-separated pair consisting of 'NSteps' and a positive scalar integer.

The simBySolution 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 simBySolution does not report the output state vector at these intermediate points, the refinement improves accuracy by allowing the simulation to more closely approximate the underlying continuous-time process.

Data Types: double

Flag to use antithetic sampling to generate the Gaussian random variates that drive the Brownian motion vector (Wiener processes), specified as the comma-separated pair consisting of 'Antithetic' and a scalar numeric or logical 1 (true) or 0 (false).

When you specify true, simBySolution performs sampling such that all primary and antithetic paths are simulated and stored in successive matching pairs:

• Odd NTrials (1,3,5,...) correspond to the primary Gaussian paths.

• Even NTrials (2,4,6,...) are the matching antithetic paths of each pair derived by negating the Gaussian draws of the corresponding primary (odd) trial.

Note

If you specify an input noise process (see Z), simBySolution ignores the value of Antithetic.

Data Types: logical

Direct specification of the dependent random noise process for generating the Brownian motion vector (Wiener process) that drives the simulation, specified as the comma-separated pair consisting of 'Z' and a function or an (NPeriods * NSteps)-by-NBrowns-by-NNTrials three-dimensional array of dependent random variates.

The input argument Z allows you to directly specify the noise generation process. This process takes precedence over the Correlation parameter of the input merton object and the value of the Antithetic input flag.

Specifically, when Z is specified, Correlation is not explicitly used to generate the Gaussian variates that drive the Brownian motion. However, Correlation is still used in the expression that appears in the exponential term of the log[Xt] Euler scheme. Thus, you must specify Z as a correlated Gaussian noise process whose correlation structure is consistently captured by Correlation.

Note

If you specify Z as a function, it must return an NBrowns-by-1 column vector, and you must call it with two inputs:

• A real-valued scalar observation time t

• An NVars-by-1 state vector Xt

Data Types: double | function

Dependent random counting process for generating the number of jumps, specified as the comma-separated pair consisting of 'N' and a function or an (NPeriodsNSteps) -by-NJumps-by-NNTrials three-dimensional array of dependent random variates. If you specify a function, N must return an NJumps-by-1 column vector, and you must call it with two inputs: a real-valued scalar observation time t followed by an NVars-by-1 state vector Xt.

Data Types: double | function

Flag that indicates how the output array Paths is stored and returned, specified as the comma-separated pair consisting of 'StorePaths' and a scalar numeric or logical 1 (true) or 0 (false).

If StorePaths is true (the default value) or is unspecified, simBySolution returns Paths as a three-dimensional time series array.

If StorePaths is false (logical 0), simBySolution returns Paths as an empty matrix.

Data Types: logical

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)$

simBySolution applies processing functions at the end of each observation period. These functions must accept the current observation time t and the current state vector Xt, and return a state vector that can be an adjustment to the input state.

The end-of-period Processes argument allows you to terminate a given trial early. At the end of each time step, simBySolution tests the state vector Xt for an all-NaN condition. Thus, to signal an early termination of a given trial, all elements of the state vector Xt must be NaN. This test enables a user-defined Processes function to signal early termination of a trial, and offers significant performance benefits in some situations (for example, pricing down-and-out barrier options).

If you specify more than one processing function, simBySolution invokes the functions in the order in which they appear in the cell array. You can use this argument to specify boundary conditions, prevent negative prices, accumulate statistics, plot graphs, and more.

Data Types: cell | function

## Output Arguments

collapse all

Simulated paths of correlated state variables, returned as an (NPeriods + 1)-by-NVars-by-NNTrials 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 StorePaths is set to false, simBySolution 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.

Dependent random variates for generating the Brownian motion vector (Wiener processes) that drive the simulation, returned as a (NPeriods * NSteps)-by-NBrowns-by-NNTrials three-dimensional time-series array.

Dependent random variates for generating the jump counting process vector, returned as an (NPeriods ⨉ NSteps)-by-NJumps-by-NNTrials three-dimensional time-series array.

collapse all

### 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 that has the same expected value but a 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 other pairs, 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 NTrials.

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

The simBySolution function simulates the state vector Xt by an approximation of the closed-form solution of diagonal drift Merton jump diffusion models. Specifically, it applies a Euler approach to the transformed log[Xt] process (using Ito's formula). In general, this is not the exact solution to the Merton jump diffusion model because the probability distributions of the simulated and true state vectors are identical only for piecewise constant parameters.

This function simulates any vector-valued merton process of the form

$d{X}_{t}=B\left(t,{X}_{t}\right){X}_{t}dt+D\left(t,{X}_{t}\right)V\left(t,{x}_{t}\right)d{W}_{t}+Y\left(t,{X}_{t},{N}_{t}\right){X}_{t}d{N}_{t}$

Here:

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

• B(t,Xt) is an NVars-by-NVars matrix of generalized expected instantaneous rates of return.

• D(t,Xt) is an NVars-by-NVars diagonal matrix in which each element along the main diagonal is the corresponding element of the state vector.

• V(t,Xt) is an NVars-by-NVars matrix of instantaneous volatility rates.

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

• Y(t,Xt,Nt) is an NVars-by-NJumps matrix-valued jump size function.

• dNt is an NJumps-by-1 counting process vector.

## References

[1] Aït-Sahalia, Yacine. “Testing Continuous-Time Models of the Spot Interest Rate.” Review of Financial Studies 9, no. 2 ( Apr. 1996): 385–426.

[2] Aït-Sahalia, Yacine. “Transition Densities for Interest Rate and Other Nonlinear Diffusions.” The Journal of Finance 54, no. 4 (Aug. 1999): 1361–95.

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

[4] Hull, John C. Options, Futures and Other Derivatives. 7th ed, Prentice Hall, 2009.

[5] Johnson, Norman Lloyd, Samuel Kotz, and Narayanaswamy Balakrishnan. Continuous Univariate Distributions. 2nd ed. Wiley Series in Probability and Mathematical Statistics. New York: Wiley, 1995.

[6] Shreve, Steven E. Stochastic Calculus for Finance. New York: Springer-Verlag, 2004.

Introduced in R2020a