Is it valid to simulate a shot noise Poisson process with discretization (R)

Question

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?

Answer 1

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