Efficient Algorithms for exponential martingales generated by SDE with jumps

Hello,

I am experimenting with jump processes, and in particular with stochastic diffrential equations of the form
dμ(t)=μ(t) λ (dν(t)−r dt )
where
dν(t)= increment of Poisson process
r dt= compensator
The result should be an exponential martinagale, i.e. a process constant on average.
I notice that naive implementation of the process quickly loose the martingale property probably owning to lack of statistics. Example

function mart!(du, u, p, t)
  du[1] = - p.λ * p.jump_rate * u[1]
end
function const_rate1(u,p,t)
  return p.jump_rate[1]
end  
function const_affect1!(integrator)
 v = copyto!(integrator.p.mart_tmp, integrator.u)
 integrator.u = v * (integrator.p.λ+1.0)
end  
mart_ode = ODEProblem(mart!, mu0, (0.0, tf), p)
const_jump = ConstantRateJump(const_rate1, const_affect1!)
mart_prob = JumpProblem(mart_ode, Direct(), const_jump)
ensemble_mart = EnsembleProblem(mart_prob)
sol_mart_prob = DifferentialEquations.solve(ensemble_mart, Tsit5(), EnsembleSerial(), trajectories=num_traj);

I would be very grateful to get advice even just reference to the numerical mathematics literature for learning how to efficiently integrate stochastic differential equations (driven by jump or point processes) whose solution enjoys the martingale property.
Many thanks!

It should be true if the ensembe is large enough? How large of an ensemble did you use?

I made my tests with ensembles of 10^3 and 10^4 trajectories. Indeed, increasing the statisics improves the results. Yet the results remain far from convergence but for a region for small values of the time (t) parameter and worsen for larger value of it. I made these tests because I am interested in more complicated examples of exponential martingales. I am wondering whether some importance sampling strategy can improve the situation (unless I am doing something very silly..). It seems to me that the problem is due to the growth of the variance of an exponential martingale.
Many thanks