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!