Numerical simulation of SDE's driven by Lévy processes (particularly stable processes)

Question

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:

1. Euler scheme: $$ X_{t+\delta t} - X_t = f(t,X_t) \delta t + g(t,X_t) (Z_{t+\delta t} - Z_t) $$

2. 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.