Simulating Jumps in EMT

I promised to show how distributions with discrete jumps can be simulated. The easiest way is to use a fine granularity dt of time. If we assume that n events occur during a time length 1, then we get the times of these events in the following function.

>function events (T,n,m=10000) ...
$  t=linspace(0,T,m); dt=T/m;
$  return t[nonzeros(random(m)<n*dt)];
 [0.1706,  0.2899,  0.3807,  0.4174,  0.8357,  0.995]

The trick is to assume an event, if R<n*dt, where R is a random number in [0,1].

Instead of this, we can use the fact that the waiting times are distributed with respect to an exponential distribution with lambda equal to n.

\(P(\text{Wait} = a) = \int\limits_0^a n e^{-nt} \, dt.\)

In statistics, this is proved by computing the limit of the probability that there is no event until t=T and then an event at time T.

\(\left( 1 – n \, dt \right)^{T/dt} \, n \, dt  \to n e^{-nT} \, dt\)

Mathematically, we should be more precise and not use the infinitesimal value dt on both sides. Instead we should compute the cumulative distribution F(a) of the waiting time and sum up these values for T from 0 to a. Then we have to find a careful argument that the limit is correct.

Now that we have the correct distribution of the waiting times, can we use it for a simulation. Obviously the problem is that we do not know how many events occur until time T. A loop is the easiest solution. But it is not efficient in EMT. So we make a guess, and update an insufficient or superfluous number of events.

>function events (T,n) ...
$  until t[-1]>T;
$  t=t|(t[-1]+cumsum(randexponential(1,floor(n*T)+1,1/n)));
$return t[1:max(nonzeros(t<1))]
 [0.380738,  0.806141,  0.851266,  0.895463,  0.96134]

The function randexponential() is already present in EMT. But it is very simple to write.

For the simulation of random data with jumps, the first method is far easier. Let us assume a development of the form

\(dS = \sigma_X X \sqrt{dt} + \sigma_Y Y J_\lambda\)

where X and Y are 0-1 normal distributions and J is a jump process. Then the easy way to simulate this is the following.

>T=3; m=1000; n=5; sx=0.02; sy=0.1;  ...
>t=linspace(0,T,m); dt=T/m; ...
>s=sx*sqrt(dt)*normal(m)+sy*normal(m)*(random(m)<n*dt);  ...
>s=cumsum(1|s); ...


Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.