I'm trying to learn how to numerically simulate SDE's of the form
$$
dX(t) = f(t,X)dt + g(t,X)dZ(t)
$$
Where Z(t) is a Lévy process with triplet $(a,0,\nu)$.
My question: what is currently considered the standard approach to numerically simulate equations like these?
I'm particularly interested in the case of stable processes, where
$$
\nu(dz) = \left[ c_{-}|z|^{-\alpha – 1}\mathbb{I}_{z<0} + c_{+}|z|^{-\alpha – 1}\mathbb{I}_{z>0} \right]dz
$$
From what I have read, I saw two possible "easier" techniques:
-
Euler scheme:
$$
X_{t+\delta t} – X_t = f(t,X_t) \delta t + g(t,X_t) (Z_{t+\delta t} – Z_t)
$$ -
Through Lévy-Itô decomposition, using truncation of small jumps + diffusion correction:
$$
X(t) = \left[f(t,X) + g(t,X)\left(a – \int_{\epsilon \leq |z| < 1} z \nu(dz) \right) \right] dt + g(t,X)\left( \int_{0 < |z| \leq \epsilon}z^2\nu(dz) \right) dW(t) + \int_{|z|>\epsilon}g(t,X_{-})z N_{\epsilon}(dz,dt)
$$
Where $N_{\epsilon}$ is the compound Poisson process truncation of jumps at an $\epsilon \in (0,1)$.
Afterwards, apply the Euler scheme to the previous truncated version of the SDE, simulating $N_{\epsilon}$ with activation $\lambda_{\epsilon} = \nu( (\epsilon, \infty))$ and jump sizes with distribution $U = \int_{-\infty}^{z} \nu(dz)/\lambda_{\epsilon}$.
I have already prepared code to use 2, because I already had code for simulating Brownian motion and compound Poisson process; however, I read in this doctoral thesis (Numerical Approximation of Stochastic Differential Equations Driven by Levy Motion with Infinitely Many Jumps) that the author was not able to prove this scheme has strong convergence for Lévy processes with infinite second moment (like stable processes); for me this is worrying. Additionally, the diffusive correction is only valid for some processes (not Gamma processes, for example) and it is inherently symmetric and cannot express possible asymmetry in the small jumps.
Regarding scheme 1 (some remark about its convergence properties on THE EULER SCHEME FOR LÉVY DRIVEN STOCHASTIC DIFFERENTIAL EQUATIONS: LIMIT THEOREMS), it seems easier to implement, but the integrating interface I use does not contemplate custom increment distributions, I would have to develop it all from scratch in a way that works with everything else.
Alternatively, I have found Simulation of non-Lipschitz stochastic differential equations driven by α-stable noise: a method based on deterministic homogenisation, but I am currently absolutely ignorant about the Marcus integral and I would rather simulate Itô integrals because I already have a lot of code for it.
Best Answer
The answer is effectively "it depends". There are multiple takes in the literature and it seems although the small-jump truncation schemes were often referred to as "easy to implement", this was a great euphemism that did not consider past very standard types of jump processes. A good example of this are tempered stable processes, for which the Lévy measure is $$\nu(dx) = \frac{c}{|x|^{\alpha +1}} \mathrm{e}^{-\theta x}dx$$ For the symmetric case; the general case is equivalent to the stable process in the question. The truncation yields a "big jump" activity (for each side of the jumps, separating the positive jumps from the negative ones) of the form $$ \lambda_{\epsilon} = \int_{\epsilon}^{\infty} \nu (dx) = \frac{c}{\alpha (\alpha -1)} \left[\frac{(\alpha -1)\mathrm{e}^{-\theta x}}{\epsilon^{\alpha}} - \frac{\theta \mathrm{e}^{-\theta x}}{\epsilon^{\alpha-1}} + \theta^\alpha \Gamma(2-\alpha, \theta \epsilon)\right]$$ Where $\Gamma(.,.)$ is the incomplete gamma function. To simulate the size of these big jumps, we need to numerically integrate and invert this function, as the methodology involves simulation by inverse transform. Additionally, there's the problem of the activity exploding when we approach zero, but that is already expected of this methodology. A more detailed description of this can be found in Kawai, Masuda (2011) - On simulation of tempered stable random variates, where the authors analyze different simulation methodologies for this type of process. The conclusion here is that an exact simulation scheme for the increments has better performance for this case, and it demystifies what one could find in many earlier texts indicating that this small jump truncation scheme technique is easy and straight forward.
Now, for the case mentioned in the question, on alpha stable processes, the convergence to the true solution is better if we use the exact simulation of the increments, as there's no need for losing precision by using approximate schemes when one can simulate the increments directly. So, from those two choices we should choose option 1.