Is it valid to simulate a shot noise Poisson process with discretization (R) I have a two-factor shot noise process as follows (sorry for the picture of text):

I want to simulate from this process. Is it valid to do so in a discretized way? The solution to these SDEs is a compound Poisson process, so I simulated it as follows:
shot_noise <- function(T, n, alpha, lambda, gamma, x0){
  dt = T/n
  dN = rpois(n, lambda*dt)
  x = c(x0)
  j = rep(0, n)
  for (i in 2:(n+1)){
    jumps = sum(rexp(dN[i-1], rate = gamma)*exp(-alpha * runif(dN[i-1], max = dt)))
    j[i-1] = jumps
    x[i] = x[i-1] - alpha*x[i-1] + jumps
  }
  return(x)
}

This would be the code to simulate x_+ and x_-, given the parameters. It seems correct to me, and the output looks like what I think it should. However, the other examples of simulating compound Poisson processes that I see tend not to discretize. I am doing so because this simulation will be used to test whether maximum likelihood estimation for this model is valid, and the observed data is discretized in that case ($\epsilon_t$ at each time-step).
Is the method I've used valid and correct? Is there a better way?
 A: Let us derive the simulation. It will be enough to consider one of the two
process $x^{+}_t$ and $x^{-}_t$, that we can denote by $x(t)$. It has the representation $$ x(t) = \sum_{\{k: \, T_k \leq t\}}
e^{- \alpha [t - T_k ]} \eta_k $$ where $T_k$ are the arrivals of the
Poisson Process. For $\Delta > 0$ we have $$ x(t + \Delta) =
\sum_{\{k: \, T_k \leq t\}} e^{- \alpha [t + \Delta - T_k ]} \eta_k
   + \sum_{\{k: \, t < T_k \leq t + \Delta\}} e^{- \alpha
[t + \Delta - T_k ]} \eta_k = e^{-\alpha \Delta} x(t) + \varepsilon . $$
In the second sum denoted as $\varepsilon$, the
arrivals $T_k$ and the jumps $\eta_k$ are independent of the history
$\{x(s);\, s \leq t \}$, and so is $\varepsilon$. So $x(t)$ is
a Markov process. Moreover, conditional on the number of
arrivals on $(t,\,t+ \Delta)$, say $J$, we know
that the r.vs $T_k - t$ have the distribution of the order statistics
of a sample of size $J$ of the uniform distribution on $(0,\,\Delta)$. So
\begin{equation}
   \varepsilon =  \sum_{j=1}^J e^{-\alpha U_j}
\eta_j' \tag{1}
\end{equation}
where the $\eta_j'$ are i.i.d. from the
distribution of $\eta_k$ and the $U_j$ are uniform on $(0,\,\Delta)$. The sum is understood as $0$ when $J =0$.
Therefore the simulation of $x(t + \Delta)$ conditional on $\{x(s);\,
s \leq t \}$ can be as follows


*

*Draw the number of jumps $J$ on $(t,\,t+ \Delta)$ from the Poisson distribution with mean
$\lambda \Delta$.

*Draw $J$ i.i.d r.vs $U_j$ from the uniform on $(0,\,\Delta)$ and $J$
i.i.d r.vs $\eta_j'$ from
the jump distribution.

*Take $x(t + \Delta) = e^{-\alpha \Delta} x(t) + \varepsilon$ where
$\varepsilon$ is obtained by (1).
If $t_i^\star$ is an increasing sequence of observation times which is
independent of the process $x(t)$, then we can simulate the sequence
$x(t_i^\star)$ by simulating $x(t_i^\star)$ conditional on the history
at $t = t_{i-1}^\star$ for $i > 0$, as in the code of the
question. The sequence $x(t_i^\star)$ is a Markov chain, and it is time-homogeneous if $t_i^\star-t_{i-1}^\star$ is constant.
Interestingly, it can be shown that when the jump distribution is
exponential with scale $1/ \gamma$, the stationary distribution of
$x(t)$ is gamma with shape $\lambda / \alpha$ and scale $1 / \gamma$.
This distribution can be used to draw the initial value
$x(t_0^\star)$, and we have an example of Markov chain with gamma margin
as in this question.
A classical reference for shot noise is section 5.6 in D.R Cox and
V. Isham Point Processes, Chapman & Hall (1980).
